如何在使用predict.gam进行gamm模型预测时排除随机斜率效应？

我试图理解当使用exclude从预测函数mgcv::predict.gam()中排除主题随机斜率项时实际发生了什么。我发现无论是I exclude还是include，预测值都是相同的。有人能帮我理解为什么吗？
一些其他上下文：我以前问过一个关于排除随机效应here的问题。在我的最后一个问题中，我被建议使用exclude来排除某些随机效应，如果使用mgcv::predict.gam。在这里，我还读到，总是建议在测试/预测数据中包含随机效应协变量，然后在预测函数中排除这些。那篇文章提到，测试数据中的随机效应协变量可以包含任意值，因为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)"可以/不能被排除？目前，因为pr和pr1是相同的，我不知道这些是被排除还是包含在预测函数中。任何指导将非常感谢。
根据上面的观点，如果我在预测数据中将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

我刚刚意识到exclude对传递给它的向量中的空格敏感。因此，如果我排除了多个随机效应，它需要完全像它在summary(fit)中出现的那样出现，而不是像它在我的原始模型规范中出现的那样。否则，不排除这些影响。
例如：

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) 

> pr1
         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

小问题，但谜团解开了！顺便说一下，这也证实了排除所有随机效应项给出的预测值与将零向量传递给测试数据中的Reference_fact协变量相同（这防止使用任何随机效应;因此仅使用总体平均值）。

如何在使用predict.gam进行gamm模型预测时排除随机斜率效应？

1条答案

相关问题

热门标签

最新问答