如何从metafor元回归结果中提取特定效应修饰符的b值与se值

问题背景 #

用R的metafor包对9个因变量（存储在metrics4中）执行元回归分析，代码如下：

output5_MR = map(metrics4,
                 #magrittr::extract(!. %in% c("Soil NPK availability", "Nutrient use efficiency")),
                 function(i) metadata1 %>%
                   dplyr::filter(measurement_n==i) %>%
                   rma.mv(lnrr, v, random = ~ 1 | publication_title / unique_id, mods = ~ duration_exp + temp_group + soil_texture + country + Biochar_app_rate, 
method = "REML", data=.)) 

注：原代码的mods参数未包含Biochar_app_rate，需添加后才能提取该变量的结果。

运行模型后，希望提取特定效应修饰符（如Biochar_app_rate）的估计值b和标准误se，但之前的代码提取了所有效应修饰符的结果，无法针对性获取目标变量的值。

解决方案 #

metafor的rma.mv模型结果对象中，b和se都是命名向量，变量名对应模型中的效应修饰符名称。因此可以通过变量名索引来提取特定变量的结果，结合purrr的map系列函数实现批量处理：

1. 提取特定变量的系数（b） #

# 提取Biochar_app_rate的b值，不存在则返回NA
output5_MR_b <- map_dbl(output5_MR, 
                        function(x) as.numeric(x$b["Biochar_app_rate"]),
                        .default = NA_real_)

2. 提取特定变量的标准误（se） #

# 提取Biochar_app_rate的se值，不存在则返回NA
output5_MR_se <- map_dbl(output5_MR, 
                         function(x) as.numeric(x$se["Biochar_app_rate"]),
                         .default = NA_real_)

3. 整理成数据框（更直观） #

如果希望把结果整合到一个数据框中，方便查看，可以用map_df：

library(tibble)

output_MR_specific <- map_df(output5_MR, 
                             function(x) {
                               tibble(
                                 metric = names(x$model$y), # 获取当前模型对应的因变量
                                 b_biochar = as.numeric(x$b["Biochar_app_rate"]),
                                 se_biochar = as.numeric(x$se["Biochar_app_rate"])
                               )
                             },
                             .default = tibble(metric = NA, b_biochar = NA, se_biochar = NA))

注意事项 #

• 如果某个因变量的子数据集中，Biochar_app_rate没有变异（比如所有值相同），rma.mv会自动从模型中移除该变量，此时提取会返回NA，添加.default = NA_real_可以避免报错。
• 原示例数据集中没有Biochar_app_rate变量，测试时需先添加该变量：

# 给示例数据集添加Biochar_app_rate变量
metadata1$Biochar_app_rate <- c(10, 20, 15, 5, 25, 18, 12, 8, 22, 30)

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