在R中如何基于混合GAM模型仅对固定效应因子正确预测？

问题背景 #

我使用library(mgcv)为包含5个存在嵌套关系的数据集拟合GAM模型，其中自变量为x，响应变量为y。分析目标是针对固定因子（示例中factor4和factor5为固定因子，其余为随机因子）的不同水平，预测y作为x的函数。希望纳入随机因子完善模型，但仅基于固定效应层面进行预测（类似re_form = NA的效果）。

数据集 #

data_exp <- read.csv(text = "
, factor1, factor2, factor3, factor4, factor5, x, y
1, A, 1, A15, DE, RS1, 1, 20.6
2, B, 1, A15, DE, RS1, 1.5, 26.2
3, C, 1, A15, DE, RS1, 1.5, 24.2
4, D, 1, A15, DE, RS1, 1.8, 31
5, E, 1, A15, DE, RS1, 2.2, 30.7
6, F, 1, A15, DE, RS1, 2.2, 28.4
7, A, 2, A15, DE, RS1, 2.2, 31
8, B, 2, A15, DE, RS1, 3.5, 36.3
9, C, 2, A15, DE, RS1, 3.6, 36.2
10, D, 2, A15, DE, RS1, 3.7, 34.7
11, E, 2, A15, DE, RS1, 5.4, 41
12, F, 2, A15, DE, RS1, 6.4, 39.6
13, A, 1, A16, DE, RS1, 6.7, 37.7
14, B, 1, A16, DE, RS1, 6.7, 34.5
15, C, 1, A16, DE, RS1, 6.7, 41.2
16, D, 1, A16, DE, RS1, 6.7, 38.9
17, E, 1, A16, DE, RS1, 6.7, 38.9
18, F, 1, A16, DE, RS1, 7.8, 41.3
19, A, 2, A16, DE, RS1, 8, 41.1
20, B, 2, A16, DE, RS1, 8.5, 40.8
21, C, 2, A16, DE, RS1, 8.7, 41.5
22, D, 2, A16, DE, RS1, 9.5, 41.3
23, E, 2, A16, DE, RS1, 10, 39.9
24, F, 2, A16, DE, RS1, 10, 35.8
25, A, 1, A15, DE, RS2, 1, 24
26, B, 1, A15, DE, RS2, 1.5, 31.1
27, C, 1, A15, DE, RS2, 1.8, 36.3
28, D, 1, A15, DE, RS2, 1.8, 30.9
29, E, 1, A15, DE, RS2, 1.8, 33
30, F, 1, A15, DE, RS2, 2.5, 37.5
31, A, 2, A15, DE, RS2, 2.7, 39.5
32, B, 2, A15, DE, RS2, 3.5, 40.2
33, C, 2, A15, DE, RS2, 3.6, 44.9
34, D, 2, A15, DE, RS2, 3.7, 40.9
35, E, 2, A15, DE, RS2, 3.7, 42
36, F, 2, A15, DE, RS2, 3.7, 43.2
37, A, 1, A16, DE, RS2, 4.5, 41.9
38, B, 1, A16, DE, RS2, 4.8, 41
39, C, 1, A16, DE, RS2, 5, 46.8
40, D, 1, A16, DE, RS2, 5, 38.6
41, E, 1, A16, DE, RS2, 5, 45
42, F, 1, A16, DE, RS2, 5.5, 43.8
43, A, 2, A16, DE, RS2, 6.2, 46.2
44, B, 2, A16, DE, RS2, 6.2, 45
45, C, 2, A16, DE, RS2, 6.2, 48.3
46, D, 2, A16, DE, RS2, 7.5, 42.2
47, E, 2, A16, DE, RS2, 7.5, 47.8
48, F, 2, A16, DE, RS2, 7.5, 48.5
49, A, 1, A15, PE, RS1, 1, 27.2
50, B, 1, A15, PE, RS1, 1.5, 34.1
51, C, 1, A15, PE, RS1, 2, 39.1
52, D, 1, A15, PE, RS1, 2.5, 39.1
53, E, 1, A15, PE, RS1, 2.5, 41.8
54, F, 1, A15, PE, RS1, 2.5, 41.5
55, A, 2, A15, PE, RS1, 2.5, 41.6
56, B, 2, A15, PE, RS1, 3.5, 47
57, C, 2, A15, PE, RS1, 4, 43.3
58, D, 2, A15, PE, RS1, 4.1, 48.4
59, E, 2, A15, PE, RS1, 4.2, 44.5
60, F, 2, A15, PE, RS1, 4.5, 44.1
61, A, 1, A16, PE, RS1, 4.5, 44.9
62, B, 1, A16, PE, RS1, 4.5, 46.9
63, C, 1, A16, PE, RS1, 4.5, 46.7
64, D, 1, A16, PE, RS1, 5, 46.2
65, E, 1, A16, PE, RS1, 5, 44.9
66, F, 1, A16, PE, RS1, 5.2, 45.7
67, A, 2, A16, PE, RS1, 5.3, 49.5
68, B, 2, A16, PE, RS1, 5.4, 47.5
69, C, 2, A16, PE, RS1, 5.5, 46
70, D, 2, A16, PE, RS1, 5.6, 45.4
71, E, 2, A16, PE, RS1, 5.8, 46.9
72, F, 2, A16, PE, RS1, 6.5, 44.3
73, A, 1, A15, PE, RS2, 1, 32.5
74, B, 1, A15, PE, RS2, 1, 30.9
75, C, 1, A15, PE, RS2, 1.8, 42.3
76, D, 1, A15, PE, RS2, 1.9, 42.4
77, E, 1, A15, PE, RS2, 2, 45.2
78, F, 1, A15, PE, RS2, 2.3, 45
79, A, 2, A15, PE, RS2, 2.3, 41.1
80, B, 2, A15, PE, RS2, 2.3, 43.2
81, C, 2, A15, PE, RS2, 2.3, 47.4
82, D, 2, A15, PE, RS2, 3.2, 50.5
83, E, 2, A15, PE, RS2, 3.5, 48.7
84, F, 2, A15, PE, RS2, 3.6, 47.8
85, A, 1, A16, PE, RS2, 3.6, 48.5
86, B, 1, A16, PE, RS2, 3.6, 47.7
87, C, 1, A16, PE, RS2, 3.7, 48
88, D, 1, A16, PE, RS2, 3.8, 48.4
89, E, 1, A16, PE, RS2, 3.9, 50.6
90, F, 1, A16, PE, RS2, 4, 48.3
91, A, 2, A16, PE, RS2, 4, 48.2
92, B, 2, A16, PE, RS2, 4, 47.8
93, C, 2, A16, PE, RS2, 4, 51.5
94, D, 2, A16, PE, RS2, 4.5, 49.2
95, E, 2, A16, PE, RS2, 4.6, 51.6
96, F, 2, A16, PE, RS2, 4.8, 51.7") %>%
  mutate(factor1 = as.factor(factor1),
         factor2 = as.factor(factor2),
         factor3 = as.factor(factor3),
         factor4 = as.factor(factor4),
         factor5 = as.factor(factor5))

尝试的代码 #

模型拟合 #

library(mgcv)
bam_ex.2 <- bam (y ~
                   factor5 +
                   factor4 + 
                   s(x, by = interaction(factor5, factor4), k= 3) + # 为两个因子的每个水平创建一个平滑项
                   s(factor1, by = factor4, bs = "re") + # factor1嵌套于factor4，为随机因子
                   s(factor2, by = factor3, bs = "re") + # factor2嵌套于factor3，为随机因子       
                   s(factor3, bs = "re"), # factor3为随机因子
                 data = data_exp,
                 discrete = TRUE)

预测与绘图 #

pred_ex.2 <- data_exp %>% 
  mutate(y.p = predict(bam_ex.2, 
                       re_form = NA,
                       type="response"))

library(ggplot2)
ggplot(pred_ex.2, aes(x = x, y = y.p, color = factor4))+
  geom_line(aes(linetype= factor5), size=1) +
  scale_linetype_manual(values=c("solid", "dotted", "dashed"))

遇到的问题 #

绘制的预测结果图出现不连续线条，检查后发现：当factor4、factor5和x的水平完全相同时，y.p的预测值却不一致，这和预期的4条连续平滑线条（factor4×factor5共4种组合）不符。

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