引导元分析模型时双峰分布的原因

机器算法验证 r 混合模式 引导程序 荟萃分析 双峰
2022-03-24 00:22:17

我帮助一位同事使用@Wolfgang 编写的metafor R 包框架引导一个元分析混合效应模型。

有趣且令人担忧的是，对于模型的一个系数，我在引导时得到双峰分布（参见下图的右下角面板）。

我想主要原因之一可能是在引导时，假设一半模型收敛于本地解决方案，另一半收敛于另一个。我尝试按照 metafor 文档中的建议调整收敛算法 - Convergence Problems with the rma() function。此外，我尝试了其他收敛算法bobyqa，如rma.mvnewuoa函数的帮助文档中所建议的，但得到了相同的双峰响应。

我还尝试按照如何解释自举相关性的多模态分布中的建议从有问题的组中消除一些潜在的异常值，但无济于事。

我找不到重现此问题的方法，因此我将数据上传到GitHub 存储库（下面代码部分中的链接也应加载到您的环境中，以测试案例所需的所有内容）。我在 Linux 集群上将引导程序作为数组作业运行（以防万一，shell 脚本是job.sh，它在每个 CPU 上执行运行下述模型的 R 脚本bootstrap.r ）。单次运行需要 2-3 分钟。请注意，自举 100 次也足以检测双峰响应。以下是 1000 次迭代的示例。我熟悉 R 和其他方法，但对元分析不太熟悉。

如果双峰分布正常（尽管可能是由于收敛问题），我将不胜感激，如果不是，那么可以做些什么呢？（除了我已经尝试过的）

下面 - 比较自举（红线）和单个完整模型运行（蓝线）的系数。直方图描绘了每个系数的自举分布。对引导数据进行抽样是通过从两个固定效应形成的每个组/组合中选择替换来完成的。他们的原始样本量是：

table(dt$f1, dt$f2)
#>       
#>        f2_1 f2_2 f2_3
#>   f1_1  177  174   41
#>   f1_2  359  363  107

library(data.table)
library(ggplot2)
library(metafor)
#> Loading required package: Matrix
#> Loading 'metafor' package (version 2.0-0). For an overview 
#> and introduction to the package please type: help(metafor).

load(url("https://github.com/valentinitnelav/test/raw/master/bimodal_distrib_boot/coef_boot_dt_1010.rda"))
load(url("https://github.com/valentinitnelav/test/raw/master/bimodal_distrib_boot/rmamv_model.rda"))
load(url("https://github.com/valentinitnelav/test/raw/master/bimodal_distrib_boot/data.rda"))

coef_dt <- data.frame(estim = rmamv_model[["beta"]])
coef_dt$coef_name <- rownames(coef_dt)
coef_dt <- rbind(coef_dt,
                 coef_boot_dt[, .(estim = mean(coef)), by = coef_name])
coef_dt[, gr := rep(c("estim_model", "estim_boot"), each = 6)]

ggplot(data = coef_boot_dt,
       aes(x = coef,
           group = coef_name)) +
  geom_histogram(bins = 100) +
  geom_vline(aes(xintercept = estim,
                 group = gr,
                 color = gr),
             lwd = 1,
             data = coef_dt) +
  facet_wrap(vars(coef_name), ncol = 2)

由reprex 包（v0.2.1）于 2019 年 5 月 2 日创建

模型是这样的：

rmamv_model <- rma.mv(y ~ f2:f1 - 1,
                  V = var_y,
                  random = list(~ 1|r1,
                                ~ 1|r2),
                  R = list(r2 = cor_mat),
                  data = dt,
                  method = "REML",
                  # Tune the convergence algorithm / optimizer
                  control = list(optimizer = "nlminb",
                                 iter.max = 1000,
                                 step.min = 0.4,
                                 step.max = 0.5))

R会话信息：

devtools::session_info()
#> - Session info ----------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.5.2 (2018-12-20)
#>  os       Windows 7 x64 SP 1          
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  English_United States.1252  
#>  ctype    English_United States.1252               
#>  date     2019-05-02                  
#> 
#> - Packages --------------------------------------------------------------
#>  package     * version date       lib source        
#>  assertthat    0.2.1   2019-03-21 [1] CRAN (R 3.5.2)
#>  backports     1.1.3   2018-12-14 [1] CRAN (R 3.5.2)
#>  callr         3.2.0   2019-03-15 [1] CRAN (R 3.5.3)
#>  cli           1.1.0   2019-03-19 [1] CRAN (R 3.5.3)
#>  colorspace    1.4-1   2019-03-18 [1] CRAN (R 3.5.3)
#>  crayon        1.3.4   2017-09-16 [1] CRAN (R 3.5.1)
#>  curl          3.3     2019-01-10 [1] CRAN (R 3.5.2)
#>  data.table  * 1.12.0  2019-01-13 [1] CRAN (R 3.5.2)
#>  desc          1.2.0   2018-05-01 [1] CRAN (R 3.5.1)
#>  devtools      2.0.1   2018-10-26 [1] CRAN (R 3.5.1)
#>  digest        0.6.18  2018-10-10 [1] CRAN (R 3.5.1)
#>  dplyr         0.8.0.1 2019-02-15 [1] CRAN (R 3.5.2)
#>  evaluate      0.13    2019-02-12 [1] CRAN (R 3.5.2)
#>  fs            1.2.7   2019-03-19 [1] CRAN (R 3.5.3)
#>  ggplot2     * 3.1.0   2018-10-25 [1] CRAN (R 3.5.1)
#>  glue          1.3.1   2019-03-12 [1] CRAN (R 3.5.2)
#>  gtable        0.2.0   2016-02-26 [1] CRAN (R 3.5.1)
#>  highr         0.8     2019-03-20 [1] CRAN (R 3.5.3)
#>  htmltools     0.3.6   2017-04-28 [1] CRAN (R 3.5.1)
#>  httr          1.4.0   2018-12-11 [1] CRAN (R 3.5.2)
#>  knitr         1.22    2019-03-08 [1] CRAN (R 3.5.2)
#>  labeling      0.3     2014-08-23 [1] CRAN (R 3.5.0)
#>  lattice       0.20-38 2018-11-04 [2] CRAN (R 3.5.2)
#>  lazyeval      0.2.2   2019-03-15 [1] CRAN (R 3.5.3)
#>  magrittr      1.5     2014-11-22 [1] CRAN (R 3.5.1)
#>  Matrix      * 1.2-15  2018-11-01 [2] CRAN (R 3.5.2)
#>  memoise       1.1.0   2017-04-21 [1] CRAN (R 3.5.1)
#>  metafor     * 2.0-0   2017-06-22 [1] CRAN (R 3.5.2)
#>  mime          0.6     2018-10-05 [1] CRAN (R 3.5.1)
#>  munsell       0.5.0   2018-06-12 [1] CRAN (R 3.5.1)
#>  nlme          3.1-137 2018-04-07 [2] CRAN (R 3.5.2)
#>  pillar        1.3.1   2018-12-15 [1] CRAN (R 3.5.2)
#>  pkgbuild      1.0.3   2019-03-20 [1] CRAN (R 3.5.3)
#>  pkgconfig     2.0.2   2018-08-16 [1] CRAN (R 3.5.1)
#>  pkgload       1.0.2   2018-10-29 [1] CRAN (R 3.5.1)
#>  plyr          1.8.4   2016-06-08 [1] CRAN (R 3.5.1)
#>  prettyunits   1.0.2   2015-07-13 [1] CRAN (R 3.5.1)
#>  processx      3.3.0   2019-03-10 [1] CRAN (R 3.5.3)
#>  ps            1.3.0   2018-12-21 [1] CRAN (R 3.5.2)
#>  purrr         0.3.2   2019-03-15 [1] CRAN (R 3.5.3)
#>  R6            2.4.0   2019-02-14 [1] CRAN (R 3.5.2)
#>  Rcpp          1.0.1   2019-03-17 [1] CRAN (R 3.5.3)
#>  remotes       2.0.2   2018-10-30 [1] CRAN (R 3.5.1)
#>  rlang         0.3.4   2019-04-07 [1] CRAN (R 3.5.3)
#>  rmarkdown     1.12    2019-03-14 [1] CRAN (R 3.5.3)
#>  rprojroot     1.3-2   2018-01-03 [1] CRAN (R 3.5.1)
#>  scales        1.0.0   2018-08-09 [1] CRAN (R 3.5.1)
#>  sessioninfo   1.1.1   2018-11-05 [1] CRAN (R 3.5.1)
#>  stringi       1.4.3   2019-03-12 [1] CRAN (R 3.5.2)
#>  stringr       1.4.0   2019-02-10 [1] CRAN (R 3.5.1)
#>  tibble        2.1.1   2019-03-16 [1] CRAN (R 3.5.3)
#>  tidyselect    0.2.5   2018-10-11 [1] CRAN (R 3.5.1)
#>  usethis       1.4.0   2018-08-14 [1] CRAN (R 3.5.1)
#>  withr         2.1.2   2018-03-15 [1] CRAN (R 3.5.1)
#>  xfun          0.5     2019-02-20 [1] CRAN (R 3.5.2)
#>  xml2          1.2.0   2018-01-24 [1] CRAN (R 3.5.1)
#>  yaml          2.2.0   2018-07-25 [1] CRAN (R 3.5.1)
2个回答

感谢提供数据和代码。我重新拟合了您正在使用的模型，并且第二个方差分量（为其cor_mat指定）漂移到一个非常大的值，这很奇怪。然而，分析这个方差分量（用profile(rmamv_model, sigma2=2)）表示没有问题，所以我不认为这是一个收敛问题。相反，我认为问题的出现是因为该模型不包括估计水平的随机效应（基本上每个元分析模型都应该包括）。所以，我建议适合：

dt$id <- 1:nrow(dt)

res <- rma.mv(y ~ f2:f1 - 1,
              V = var_y,
              random = list(~ 1|r1,
                            ~ 1|r2, 
                            ~ 1|id),
              R = list(r2 = cor_mat),
              data = dt,
              method = "REML")

结果看起来更合理。我怀疑这也可能解决最后一个系数的双峰自举分布问题。

如果无法访问可重现的示例，则很难对这种引导行为给出明确的答案。假设确实没有异常值，我怀疑我们观察到了斯坦因现象的一个温和案例，特别是作为混合效应方法表明我们在数据中存在一些聚类。

如上所述，我建议继续从f2f2_3:f1f1_2交互的“不寻常”值中查看一些运行，其中存在非常不同的值，并调查这两个随机子样本的边际分布。例如在某些情况下，f2f2_3:f1f1_2远低于1而估计的模型建议的值接近 2.4. 边际分布是否相似？是否存在重叠不足的情况？也许“简单”的引导是不合适的，我们需要根据某些因素对手头的样本进行分层。

其它你可能感兴趣的问题
上一篇变压器网络中位置余弦编码背后的直觉是什么？ 下一篇我应该期望亚当工作到什么程度？