在gam中你到底是怎么预测的？具有可重现性的例子

Question

当你安装了一个可能包含随机效应的模型时，你如何在mgcv::gam中预测呢？

这个站点上的另一个使用“排除”技巧的线程不适合我(https://stats.stackexchange.com/questions/131106/predicting-with-random-effects-in-mgcv-gam)

ya <- rnorm(100, 0, 1)
yb <- rnorm(100,0,1.5)
yc <- rnorm(100, 0, 2)
yd <- rnorm(100, 0, 2.5)

yy <- c(ya,yb,yc,yd) #so, now we've got data from 4 different groups. 
xx <- c(rep("a", 100), rep("b",100), rep("c",100),rep("d",100)) #groups
zz <- rnorm(400,0,1) #some other covariate

model <- gam(yy ~ zz + s(xx, bs = "re")) #the model

predictdata <- data.frame( zz = 5 )   #new data
predict(model, newdata = predictdata, exclude = "s(xx)")   #prediction

这就产生了错误

Error in model.frame.default(ff, data = newdata, na.action = na.act) : 
  variable lengths differ (found for 'xx')
In addition: Warning messages:
1: In predict.gam(model, newdata = predictdata, exclude = "s(xx)") :
  not all required variables have been supplied in  newdata!

2: 'newdata' had 1 row but variables found have 400 rows 

我的mgcv包裹是最新的。

编辑：

如果将预测数据更改为

predictdata <- data.frame(zz = 5, xx = "f")

然后上面写着

Error in predict.gam(model, newdata = predictdata, exclude = "s(xx)") : 
  f not in original fit

Answer 1

我对您的示例进行了实验，而且“排除”语句似乎确实有效，尽管您必须在新数据值中指定用于拟合模型的原始数据集中包含的随机效果。然而，这让我有点不安。另一个警告是，“排除”似乎不适用于由组单独估算的方差结构模型(我用另一个数据集尝试过)，即类似s(xx，s="re"，by=group)的模型。您可能希望发布问题或将问题转移到交叉验证，以便其他统计学家/分析人员能够看到它，也许可以提供一个更好的答案。

下面是我的密码。注意，我改变了组a和d的平均值，但是总体平均值应该在零左右。

ya <- rnorm(100, 1, 1)
yb <- rnorm(100, 0,1.5)
yc <- rnorm(100, 0, 2)
yd <- rnorm(100, -1, 2.5)

yy <- c(ya,yb,yc,yd) #so, now we've got data from 4 different groups. 
xx <- c(rep("a", 100), rep("b",100), rep("c",100),rep("d",100)) #groups
zz <- rnorm(400,0,1) #some other covariate

some.data= data.frame(yy,xx,zz)
model <- gam(yy ~ zz + s(xx, bs = "re"),data=some.data) #the model


# the intercept is the overall mean when zz is zero
summary(model)

 predictdata <- data.frame(zz = c(0,0,0,0), xx =c("a","b","c","d"))  #new data

#excluding random effects. Estimate should be the same for all and should be the intercept  
predict(model, newdata = predictdata, exclude = "s(xx)") 

#including random effects. Estimates should differ by group with 'a' larger and 'd' smaller
predict(model, newdata = predictdata)