R中约束性GAM实现：构建整体递增且满足阈值限制的模型

需求说明 #

• 模型整体需呈现递增趋势，无需严格单调
• 约束条件：当变量B超过15000时，变量A不能为0；当B小于0时，A可以为0
• 当前问题：其中一组XY配对的模型在B取最大值时，A的拟合值出现下降
• 要求：需通用函数适配多组XY配对，不能为每组单独设置权重

示例数据与当前代码进度 #

步骤1：加载包并读取数据 #

library(mgcv)
library(dplyr)

df <- structure(list(X = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2), Y = c(1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2), grp = c("0", "1", "2", "3", "4", "5", "6", 
"7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", 
"18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", 
"29", "30", "31", "32", "33", "34", "35", "0", "1", "2", "3", 
"4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", 
"16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", 
"27", "28", "29", "30", "31", "32", "33", "34", "35"), A = c(1.13659381866455, 
1.13659477233887, 1.13658058643341, 1.13659477233887, 1.13656425476074, 
1.13656532764435, 1.13658249378204, 1.13656723499298, 0.709016442298889, 
0.708269774913788, 0.707276940345764, 0.705555200576782, 0.703566908836365, 
0.705580234527588, 0.711442351341248, 0.709806680679321, 0.706768035888672, 
0.710436880588531, 0.711900770664215, 0.627014517784119, 0.624246716499329, 
0.637109100818634, 0.6302330493927, 0.712375164031982, 0.632026135921478, 
0.707980036735535, 1.02401030063629, 1.24682903289795, 1.16107773780823, 
1.72039294242859, 1.65212118625641, 1.77960109710693, 2.14993286132812, 
2.56561064720154, 2.75121307373047, 3.24887132644653, 0, 0.471771329641342, 
0.711109101772308, 0.79891711473465, 0.931526362895966, 0.968687295913696, 
0.995432734489441, 1.0633533000946, 1.16641199588776, 1.34169447422028, 
1.49520266056061, 1.63908088207245, 1.7665741443634, 1.88498985767365, 
1.99606990814209, 2.09789752960205, 2.19571089744568, 2.28780961036682, 
2.37474989891052, 2.44801378250122, 2.52833318710327, 2.60511207580566, 
2.67927670478821, 2.7647979259491, 2.82285761833191, 2.97204279899597, 
3.10502099990845, 3.23239874839783, 3.33626580238342, 3.47603249549866, 
3.57603311538696, 3.68949151039124, 3.95258760452271, 3.93370676040649, 
4.04461479187012, 4.92624521255493), B = c(300, 400, 530, 600, 
730, 770, 800, 880, 1000, 1250, 1500, 1750, 2000, 2250, 2500, 
2750, 3000, 3250, 3500, 3750, 4000, 4250, 4500, 4750, 5000, 5500, 
6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500, 10000, 15000, 
300, 400, 530, 600, 730, 770, 800, 880, 1000, 1250, 1500, 1750, 
2000, 2250, 2500, 2750, 3000, 3250, 3500, 3750, 4000, 4250, 4500, 
4750, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500, 
10000, 15000)), row.names = c(NA, -72L), class = c("tbl_df", 
"tbl", "data.frame"))

步骤2：为每个XY配对构建GAM模型 #

new_dat <- df %>% 
  summarise(gam_model=list(gam(B~s(A))), .by = c(X,Y))

步骤3：对A=0时进行逐行GAM预测 #

pred_df <- data.frame(A= 0)
gam_predict <- function(m) predict(m,pred_df)
new_dat$gam_pred_0 <- sapply(new_dat$gam_model,gam_predict)

内容的提问来源于stack exchange，提问作者user11057680