使用predict.gam排除随机斜率效应时预测值无变化的问题咨询

我在尝试用exclude参数从mgcv::predict.gam()中排除个体随机斜率项时，发现无论是否排除这些项，预测值完全一致，想请教背后的原因。

此前我了解到，使用mgcv::predict.gam时可通过exclude排除特定随机效应；另有资料指出，预测数据需包含随机效应协变量（可设任意值，predict.gam会忽略它们）。但我发现这一逻辑不适用于建模为随机斜率项的协变量——协变量的值会通过随机斜率项影响预测值。

示例代码 #

# create training data
set.seed(123)
x <- seq(1,10,1)
resampled <- sample(x,100,replace=T)

dat <- tibble(
  Reference_fact = as.factor(resampled),
  Country_fact = as.factor(rep(stri_rand_strings(10,2),10)),
  exp_A = runif(100,0,500),
  exp_B = runif(100,0,20),
  resp = exp_A^2 + 0.1*exp_B * 0.3*exp_A*exp_B + rnorm(100, 0, 10)
)
# define model specification
model_spec <- c("resp ~ 
                 s(exp_A) +
                 s(exp_B) +
                 s(Reference_fact, bs = 're') + 
                 s(Country_fact, bs = 're') +
                 s(exp_A, Reference_fact,  bs = 're') +
                 s(exp_B, Reference_fact, bs = 're')"
)

fit <- mgcv::gam(formula(str_replace_all(model_spec, "[\r\n]", "")),
                 method = 'REML',
                 family = 'gaussian',
                 data=dat)

set.seed(123)
dat_pred <- tibble(
  Reference_fact = as.factor(rep(seq(1,10,1),5)),
  Country_fact = as.factor(rep(stri_rand_strings(10,2),5)), 
  exp_A = runif(50,0,500),
  exp_B = runif(50,0,20)
)

pr <- mgcv::predict.gam(
  fit, 
  dat_pred, 
  exclude = c("s(Reference_fact)"),
  include=c(
              "s(exp_A, Reference_fact)",
              "s(exp_B, Reference_fact)"
              ),
  keep_prediction_data = FALSE, 
  newdata.guaranteed = TRUE, 
  se.fit = FALSE) 

pr1 <- mgcv::predict.gam(
  fit, 
  dat_pred, 
  exclude = c("s(Reference_fact)",
    "s(exp_A, Reference_fact)", 
    "s(exp_B, Reference_fact)"
  ),
  keep_prediction_data = FALSE, 
  newdata.guaranteed = TRUE, 
  se.fit = FALSE) 

预测结果对比 #

我预期pr和pr1的输出不同，但实际二者完全一致：

> pr
          1           2           3           4           5 
200516.2393 121796.3625 104512.3974 247897.1752 108595.6892 
          6           7           8           9          10 
126261.5061  74434.1023  89602.9046  20708.9620   5000.3689 
         11          12          13          14          15 
232978.0968 205927.3882 120337.7975 160699.1463   -467.0361 
         16          17          18          19          20 
 57676.7601 147143.7946  13544.9029  27155.7492  13026.8152 
         21          22          23          24          25 
  4670.3186  44073.4473  42882.7365  34971.4159   5406.6882 
         26          27          28          29          30 
  4348.2832  14815.5571  54557.0164  17806.6748 185469.5308 
         31          32          33          34          35 
  1041.5902  49111.2511 160814.6226   5908.1544  79419.0126 
         36          37          38          39          40 
 12468.5101   5869.1082 143593.8867 201719.6693  34835.5460 
         41          42          43          44          45 
113740.5525   1948.7769  36794.7991  21055.7098 168001.1933 
         46          47          48          49          50 
 50446.2849 165451.7506 168251.5211 159329.6230  48757.5280

> pr1
          1           2           3           4           5 
200516.2393 121796.3625 104512.3974 247897.1752 108595.6892 
          6           7           8           9          10 
126261.5061  74434.1023  89602.9046  20708.9620   5000.3689 
         11          12          13          14          15 
232978.0968 205927.3882 120337.7975 160699.1463   -467.0361 
         16          17          18          19          20 
 57676.7601 147143.7946  13544.9029  27155.7492  13026.8152 
         21          22          23          24          25 
  4670.3186  44073.4473  42882.7365  34971.4159   5406.6882 
         26          27          28          29          30 
  4348.2832  14815.5571  54557.0164  17806.6748 185469.5308 
         31          32          33          34          35 
  1041.5902  49111.2511 160814.6226   5908.1544  79419.0126 
         36          37          38          39          40 
 12468.5101   5869.1082 143593.8867 201719.6693  34835.5460 
         41          42          43          44          45 
113740.5525   1948.7769  36794.7991  21055.7098 168001.1933 
         46          47          48          49          50 
 50446.2849 165451.7506 168251.5211 159329.6230  48757.5280  

核心疑问 #

为什么s(exp_A, Reference_fact)、s(exp_B, Reference_fact)这类随机斜率项无法被exclude参数排除？目前因两次预测结果完全相同，我无法确认这些项是否被纳入预测，恳请指导。

补充问题：Reference_fact设为任意值的情况 #

若在预测数据中将Reference_fact设为任意值0，即使在预测时排除所有随机效应项，预测值也会不同。我的理解是：Reference_fact=0在训练数据中不存在，没有对应的随机斜率系数可用，因此仅使用总体均值进行预测，这个理解是否正确？

对应代码及输出 #

dat_pred_arbitrary <- tibble(
  Reference_fact = 0,
  Country_fact = as.factor(rep(stri_rand_strings(10,2),5)), 
  exp_A = runif(50,0,500),
  exp_B = runif(50,0,20)
)

pa <- mgcv::predict.gam(
  fit, 
  dat_pred_arbitrary, 
  exclude = c("s(Reference_fact)",
    "s(exp_A, Reference_fact)",
    "s(exp_B, Reference_fact)"
  ),
  keep_prediction_data = FALSE, 
  newdata.guaranteed = TRUE, 
  se.fit = FALSE) 

> pa
         1          2          3          4          5 
200488.956 121788.106 104520.843 247897.029 108629.512 
         6          7          8          9         10 
126259.139  74472.436  89594.384  20708.678   4998.486 
        11         12         13         14         15 
232969.281 205918.624 120342.779 160626.310   -459.715 
        16         17         18         19         20 
 57672.117 147243.220  13532.432  27155.012  13023.850 
        21         22         23         24         25 
  4665.593  44064.876  42886.827  34910.752   5429.489 
        26         27         28         29         30 
  4346.282  14894.535  54555.745  17806.292 185460.863 
        31         32         33         34         35 
  1019.899  49106.885 160820.450   5820.028  79453.377 
        36         37         38         39         40 
 12459.032   5961.430 143585.435 201719.366  34833.051 
        41         42         43         44         45 
113706.734   1944.823  36795.537  20968.192 168052.472 
        46         47         48         49         50 
 50444.755 165507.188 168238.245 159329.166  48750.677

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