predict.rma() 到 metafor 中的复杂新数据（多项式和因子水平）
predict.rma() onto complex new data in metafor (polynomials and factor levels)
我有一个混合效应元分析模型 (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()
但是假设我真的希望这些线条非常平滑（点之间没有直线）。上面 link 中的示例展示了如何为单个连续调节器执行此操作，方法是创建一个包含大量 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 帮助文件中的示例，以及上面 link 中的示例）。
任何关于如何使这项工作的建议将不胜感激，提前致谢！
这将完成工作：
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()
这里有一个使用 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)
我有一个混合效应元分析模型 (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()
但是假设我真的希望这些线条非常平滑（点之间没有直线）。上面 link 中的示例展示了如何为单个连续调节器执行此操作，方法是创建一个包含大量 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 帮助文件中的示例，以及上面 link 中的示例）。
任何关于如何使这项工作的建议将不胜感激，提前致谢！
这将完成工作：
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()
这里有一个使用 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)