如何用metafor的predict.rma()指定部分调节变量时对grade做元分析？

使用metafor包进行多变量元分析时，希望通过predict.rma()预测新调节变量的效应，但不想为调节变量grade指定具体值，而是让函数自动对grade的所有水平进行整合分析。

使用示例数据集dat.bangertdrowns2004拟合的模型代码如下：

# 加载metafor包
library(metafor)
dat <- dat.bangertdrowns2004

# 拟合模型
res <- rma.mv(yi, vi, mods = ~ grade + length + pers - 1, random = ~ 1 | ni,  data=dat)
res

预测时指定length=4:10和pers=1，代码如下：

# 预测
gridpred <- expand.grid(length=4:10, pers=1)
pred.res <- predict.rma(object = res , newmods = as.matrix(gridpred))

运行后报错：

Error in predict.rma(object = res, newmods = as.matrix(gridpred)) : 
  Argument 'newmods' does not specify values for this variable: grade

请问如何实现predict.rma()在上述条件下对grade进行整合分析？

要实现不用指定grade具体值、同时整合它的效应来预测，核心就是计算边际均值——也就是固定length和pers的取值后，把grade所有水平的预测值整合起来求平均。这里有两种可行的办法：

方法1：手动生成全网格预测后求平均 #

1. 先把grade的所有水平都放进预测网格里，搭配你指定的length和pers
2. 对每个组合做预测
3. 按length和pers分组，算出grade各水平的平均预测值和置信区间

代码示例：

# 构建包含所有grade水平的完整预测网格
grid_full <- expand.grid(grade = unique(dat$grade), length = 4:10, pers = 1)

# 用model.matrix生成符合模型要求的newmods格式
newmods_mat <- model.matrix(~ grade + length + pers - 1, data = grid_full)

# 执行预测
pred_full <- predict.rma(res, newmods = newmods_mat)

# 把预测结果和网格信息合并成数据框
pred_full_df <- cbind(grid_full, pred_full)

# 按length和pers分组，计算平均预测值、标准误和置信区间
library(dplyr)
pred_avg <- pred_full_df %>%
  group_by(length, pers) %>%
  summarize(
    yi_avg = mean(pred),
    # 同时考虑grade内部的抽样变异和grade之间的效应变异
    se_avg = sqrt(mean(se^2) + var(pred)/n()),
    ci_lb = yi_avg - 1.96*se_avg,
    ci_ub = yi_avg + 1.96*se_avg
  )

print(pred_avg)

方法2：把grade改成随机效应 #

如果你觉得grade的效应是随机变化的（比如你想关注所有grade的平均效应，而非某几个特定年级的差异），可以把grade从固定效应移到随机效应部分。这样预测时不用指定grade，模型会自动整合它的随机变异。

重构后的模型代码：

# 将grade设为随机效应，只把length和pers作为固定效应
res_random_grade <- rma.mv(yi, vi, mods = ~ length + pers - 1, 
                           random = ~ 1 | ni + grade,  # 添加grade为随机效应层
                           data = dat)

# 直接用指定的length和pers做预测
gridpred <- expand.grid(length=4:10, pers=1)
newmods_mat <- model.matrix(~ length + pers - 1, data = gridpred)
pred_res <- predict.rma(res_random_grade, newmods = newmods_mat)

print(pred_res)

怎么选这两种方法？ #

• 要是你需要保留grade各水平的固定差异（比如明确区分小学、中学、大学的效应），就用方法1算边际均值
• 要是你关注的是所有grade的平均效应，把grade当成随机变异的因素，就用方法2重构模型

内容的提问来源于stack exchange，提问作者luquintanilla