predict.rma()到metafor中的复杂新数据(多项式和因子水平)

zaq34kh6 于 2023-10-13 发布在其他

关注(0)|答案(2)|浏览(136)

我有一个混合效应荟萃分析模型（rma.mv），包含一个连续多项式调节器和一个分类调节器。我想预测新的数据用于绘图目的。
作为一个可复制的例子，我们可以复制可用的框架：http://www.metafor-project.org/doku.php/tips:non_linear_meta_regression

dat <- structure(list(yi = c(0.99, 0.54, -0.01, 1.29, 0.66, -0.12, 1.18,
-0.23, 0.03, 0.73, 1.27, 0.38, -0.19, 0.01, 0.31, 1.41, 1.32, 1.22, 1.24,
-0.05, 1.17, -0.17, 1.29, -0.07, 0.04, 1.03, -0.16, 1.25, 0.27, 0.27, 0.07,
0.02, 0.7, 1.64, 1.66, 1.4, 0.76, 0.8, 1.91, 1.27, 0.62, -0.29, 0.17, 1.05,
-0.34, -0.21, 1.24, 0.2, 0.07, 0.21, 0.95, 1.71, -0.11, 0.17, 0.24, 0.78,
1.04, 0.2, 0.93, 1, 0.77, 0.47, 1.04, 0.22, 1.42, 1.24, 0.15, -0.53, 0.73,
0.98, 1.43, 0.35, 0.64, -0.09, 1.06, 0.36, 0.65, 1.05, 0.97, 1.28), vi =
c(0.018, 0.042, 0.031, 0.022, 0.016, 0.013, 0.066, 0.043, 0.092, 0.009,
0.018, 0.034, 0.005, 0.005, 0.015, 0.155, 0.004, 0.124, 0.048, 0.006, 0.134,
0.022, 0.004, 0.043, 0.071, 0.01, 0.006, 0.128, 0.1, 0.156, 0.058, 0.044,
0.098, 0.154, 0.117, 0.013, 0.055, 0.034, 0.152, 0.022, 0.134, 0.038, 0.119,
0.145, 0.037, 0.123, 0.124, 0.081, 0.005, 0.026, 0.018, 0.039, 0.062, 0.012,
0.132, 0.02, 0.138, 0.065, 0.005, 0.013, 0.101, 0.051, 0.011, 0.018, 0.012,
0.059, 0.111, 0.073, 0.047, 0.01, 0.007, 0.055, 0.019, 0.104, 0.056, 0.006,
0.094, 0.009, 0.008, 0.02 ), xi = c(9.4, 6.3, 1.9, 14.5, 8.4, 1.8, 11.3,
4.8, 0.7, 8.5, 15, 11.5, 4.5, 4.3, 4.3, 14.7, 11.4, 13.4, 11.5, 0.1, 12.3,
1.6, 14.6, 5.4, 2.8, 8.5, 2.9, 10.1, 0.2, 6.1, 4, 5.1, 12.4, 10.1, 13.3,
12.4, 7.6, 12.6, 12, 15.5, 4.9, 0.2, 6.4, 9.4, 1.7, 0.5, 8.4, 0.3, 4.3, 1.7,
15.2, 13.5, 6.4, 3.8, 8.2, 11.3, 11.9, 7.1, 9, 9.9, 7.8, 5.5, 9.9, 2.6,
15.5, 15.3, 0.2, 3.2, 10.1, 15, 10.3, 0, 8.8, 3.6, 15, 6.1, 3.4, 10.2, 10.1,
13.7)), class = "data.frame", row.names = c(NA, -80L))

然后添加一个因子：

dat$fac<-rep(letters[1:3],length=length(dat$xi))

然后运行一个三次多项式的模型（好吧，它不是rma.mv，因为我找不到一个方便的混合效应例子，但rma应该是一样的）：

res.cub<-rma(yi,vi,mods=~poly(xi,3,raw=T)*fac,data=dat)

我可以毫无问题地预测现有的调节值，如果我将这些附加到原始的矩阵中，我可以为不同的因素绘制不同的线条：

mypreds<-predict(res.cub,addx=TRUE)
dat$pred<-mypreds$pred    
ggplot(data=dat,aes(x=xi,y=pred,colour=fac,group=fac)) +
 geom_line() +
 geom_point() +
 theme_bw()

但是，假设我真的希望这些线是超级平滑的（点之间没有直线）。上面链接中的例子展示了如何为一个连续的缓和剂做这件事，通过创建一个有很多xi值的向量并预测这些值。然而，我不能正确地设置预测网格来为因子的每个水平做这件事。
我知道因子需要在predict.rma中指定为虚拟变量，所以我认为下面生成了我需要的新数据矩阵：

xs<-seq(-1, 16, length=50)
newmods1<-as.matrix(expand.grid(unname(poly(xs, degree=3, raw=TRUE)),unname(rep(c(0,1))),
                                unname(rep(c(0,1))),unname(rep(c(0,1)))))
lotsofpreds<-predict(res.cub, newmods=newmods1)

但我得到了“错误：在模型中找不到变量'Var 1'。”
我想我的问题是与新矩阵中变量的名称有关，但我找不到如何处理这个问题的例子，因为大多数其他例子似乎都没有在newmods参数中使用命名变量（例如：在predict.rma帮助文件中的例子，以及上面链接中的例子）。
任何关于如何让这个工作的建议将不胜感激，提前感谢！

来源：https://stackoverflow.com/questions/63554740/predict-rma-onto-complex-new-data-in-metafor-polynomials-and-factor-levels

2条答案

按热度按时间

30byixjq1#

这里有一个更通用的解决方案，使用model.matrix（），以防其他人从rma模型中进行奇怪而精彩的预测：

#model (from original post):
res.cub<-rma(yi,vi,mods=~poly(xi,3,raw=T)*fac,data=dat)
#copy and save model formula:
newform <- (~poly(xi,3,raw=T)*fac) 
#create new x values wanted:
n  <- 50
xs <- seq(-1, 16, length=n) 
#generate x values for all levels of factors using same term names in same order as in model formula
newgrid <- data.frame(expand.grid(xi=xs,fac=levels(as.factor(dat$fac)))) 
#create the new model matrix and remove the intercept
predgrid<-model.matrix(newform,data=newgrid)[,-1] 
#predict onto the new model matrix
mypreds <- predict(res.cub, newmods=predgrid) 
#attach predictions to variables for plotting
newgrid$pred<-mypreds$pred 
#plot
ggplot(data=newgrid, aes(x=xi, y=pred, colour=fac, group=fac)) + 
  geom_line() +
  geom_point() +
  theme_bw()

作为奖励，如果一个人正在做调整后的估计或边际均值（http://www.metafor-project.org/doku.php/tips:computing_adjusted_effects），并希望在某些变量的水平上求平均值，那么可以在预测新矩阵之前插入以下代码片段：

#turn matrix temporarily into a dataframe
predgrid2<-as.data.frame(predgrid)
#run a for loop to replace all columns containing 'fac' with the column mean
for (colname in colnames(predgrid2)) {
  if (grepl('fac', colname)){
    predgrid2[[colname]] <- colMeans(predgrid2[colname])
  }
}
#remove duplicates (combinations of x with different fac levels no longer exist)
predgrid2<-predgrid2[!duplicated(predgrid2),]
#turn back into a matrix so predict.rma can read it
predgrid<-as.matrix(predgrid2)

展开查看全部

赞(0）回复(0）举报 2023-10-13

a5g8bdjr2#

这将完成工作：

n  <- 50
xs <- seq(-1, 16, length=n)
xs <- cbind(xs, xs^2, xs^3)
grp.a <- cbind(xs, 0,0, xs*0, xs*0)
grp.b <- cbind(xs, 1,0, xs*1, xs*0)
grp.c <- cbind(xs, 0,1, xs*0, xs*1)
newmods <- rbind(grp.a, grp.b, grp.c)
mypreds <- predict(res.cub, newmods=newmods)
mypreds <- data.frame(x = newmods[,1], pred = mypreds$pred, fac = rep(c("a","b","c"), each=n))
ggplot(data=mypreds, aes(x=x, y=pred, colour=fac, group=fac)) +
 geom_line() +
 geom_point() +
 theme_bw()

展开查看全部

赞(0）回复(0）举报 2023-10-13

我来回答

predict.rma()到metafor中的复杂新数据(多项式和因子水平)

2条答案

相关问题

热门标签

最新问答