1

我有一个混合效应元分析模型 (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)

但我得到“错误：在模型中找不到变量‘Var1’。”

我想我的问题与新矩阵中变量的名称有关，但我找不到如何处理这个问题的示例，因为大多数其他示例似乎没有在 newmods 参数中使用命名变量（例如predict.rma 帮助文件中的示例，以及上面链接中的示例）。

任何有关如何使这项工作的建议将不胜感激，在此先感谢！

4

2 回答 2

1

这将完成工作：

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()
于 2020-08-24T15:08:25.563 回答
1

这是一个使用 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)
于 2020-10-23T02:58:25.613 回答