Random effect predictions from gamm model error: cannot evaluate groups for desired levels on 'newdata'
Random effect predictions from gamm model error: cannot evaluate groups for desired levels on 'newdata'
我正在尝试使用 newdata 参数从 gamm 模型（来自 mgcv 包）生成预测。我想对模型的 lme 部分进行预测，以便预测包括 运行dom 效应。但是，我认为，由于模型系数的命名方式，我 运行 遇到了问题。
我的问题是，newdata 参数应该如何构造/命名以允许预测。谢谢。
阿姆维
mod <- gamm(outcome ~ s(time) + predvar, data=d,
random=list(groupvar=~1),
correlation = corARMA(form=~1|groupvar, p = 1))
# okay
pred <- predict(mod$lme)
# Not okay
pred <- predict(mod$lme, newdata=d)
产生错误
Error in predict.lme(mod$lme, newdata = d) :
cannot evaluate groups for desired levels on 'newdata'
如果我 运行 nlme 中的模型没有样条项，newdata 执行没有问题
mod2 <- lme(outcome ~ time + predvar, data=d,
random=list(groupvar=~1),
correlation = corARMA(form=~1|groupvar, p = 1))
# okay
pred2 <- predict(mod2, newdata=d)
d <- structure(list(time = c(0, 1, 2, 3, 4, 5, 6, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16, 17, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14), outcome = c(-1.85, -1.57, -1.38, -1.22, -1.27, -1.63,
-2.07, -1.36, -0.33, 0.08, 0.3, 0.44, 0.78, 1.03, 1.13, 1.14,
1.05, 0.94, 0.73, 0.51, 0.08, 0.01, 0.42, 0.59, 0.71, 0.79, 0.87,
0.75, 0.6, 0.38, 0.01, -0.63), predvar = c(-1.83, -1.77, -1.7,
-1.84, -1.84, -1.72, -1.69, 0.01, -0.07, 0.16, -0.04, 0.04, 0.25,
0.19, 0.17, 0.22, 0.34, 0.54, 0.7, 0.81, 0.92, 1.12, 0.58, 0.63,
0.63, 0.68, 0.62, 0.56, 0.61, 0.73, 0.92, 1.07), groupvar = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("a",
"b", "c"), class = "factor")), .Names = c("time", "outcome",
"predvar", "groupvar"), row.names = c(NA, -32L), class = "data.frame")
信息：我没有将 运行dom 效果指定为样条曲线 (s(. , bs="re")) 因为我的 RE 比上面的例子更复杂。
如果您需要随机效应，对新数据进行预测的一种方法是对模型的 gam 部分进行预测，然后添加随机效应。
使用上面的例子，
library(mgcv)
mod <- gamm(outcome ~ s(time) + predvar, data=d,
random=list(groupvar=~1),
correlation = corARMA(form=~1|groupvar, p = 1))
# For comparison: predict with RE: we cant use the newdata arg here
pred <- predict(mod$lme)
# Extract the random effects from the model and match with the relevant observation
re <- coef(mod$lme)[ncol(coef(mod$lme))]
pred_ref <- re[[1]][match(d$groupvar, gsub(".*/", "", rownames(re)) )]
# Predict on gam part of model and adjust for RE
pred2 <- as.vector(predict(mod$gam, data=d) - pred_ref)
# Compare
all.equal(pred, pred2, check.attributes = F, use.names = F)
我正在尝试使用 newdata 参数从 gamm 模型（来自 mgcv 包）生成预测。我想对模型的 lme 部分进行预测，以便预测包括 运行dom 效应。但是，我认为，由于模型系数的命名方式，我 运行 遇到了问题。
我的问题是，newdata 参数应该如何构造/命名以允许预测。谢谢。
阿姆维
mod <- gamm(outcome ~ s(time) + predvar, data=d,
random=list(groupvar=~1),
correlation = corARMA(form=~1|groupvar, p = 1))
# okay
pred <- predict(mod$lme)
# Not okay
pred <- predict(mod$lme, newdata=d)
产生错误
Error in predict.lme(mod$lme, newdata = d) : cannot evaluate groups for desired levels on 'newdata'
如果我 运行 nlme 中的模型没有样条项，newdata 执行没有问题
mod2 <- lme(outcome ~ time + predvar, data=d,
random=list(groupvar=~1),
correlation = corARMA(form=~1|groupvar, p = 1))
# okay
pred2 <- predict(mod2, newdata=d)
d <- structure(list(time = c(0, 1, 2, 3, 4, 5, 6, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16, 17, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14), outcome = c(-1.85, -1.57, -1.38, -1.22, -1.27, -1.63,
-2.07, -1.36, -0.33, 0.08, 0.3, 0.44, 0.78, 1.03, 1.13, 1.14,
1.05, 0.94, 0.73, 0.51, 0.08, 0.01, 0.42, 0.59, 0.71, 0.79, 0.87,
0.75, 0.6, 0.38, 0.01, -0.63), predvar = c(-1.83, -1.77, -1.7,
-1.84, -1.84, -1.72, -1.69, 0.01, -0.07, 0.16, -0.04, 0.04, 0.25,
0.19, 0.17, 0.22, 0.34, 0.54, 0.7, 0.81, 0.92, 1.12, 0.58, 0.63,
0.63, 0.68, 0.62, 0.56, 0.61, 0.73, 0.92, 1.07), groupvar = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("a",
"b", "c"), class = "factor")), .Names = c("time", "outcome",
"predvar", "groupvar"), row.names = c(NA, -32L), class = "data.frame")
信息：我没有将 运行dom 效果指定为样条曲线 (s(. , bs="re")) 因为我的 RE 比上面的例子更复杂。
如果您需要随机效应，对新数据进行预测的一种方法是对模型的 gam 部分进行预测，然后添加随机效应。
使用上面的例子，
library(mgcv)
mod <- gamm(outcome ~ s(time) + predvar, data=d,
random=list(groupvar=~1),
correlation = corARMA(form=~1|groupvar, p = 1))
# For comparison: predict with RE: we cant use the newdata arg here
pred <- predict(mod$lme)
# Extract the random effects from the model and match with the relevant observation
re <- coef(mod$lme)[ncol(coef(mod$lme))]
pred_ref <- re[[1]][match(d$groupvar, gsub(".*/", "", rownames(re)) )]
# Predict on gam part of model and adjust for RE
pred2 <- as.vector(predict(mod$gam, data=d) - pred_ref)
# Compare
all.equal(pred, pred2, check.attributes = F, use.names = F)