J'aide un collègue à amorcer un modèle d'effets mixtes de méta-analyse en utilisant le framework de package metafor R créé par @Wolfgang.
Fait intéressant et inquiétant, pour l'un des coefficients du modèle, j'obtiens une distribution bimodale lors du bootstrap (voir le panneau en bas à droite de la figure ci-dessous).
Je suppose que l'une des principales causes pourrait être le fait que lors du démarrage, disons que la moitié des modèles convergent dans une solution locale et l'autre moitié dans une autre. J'ai essayé de régler l'algorithme de convergence comme suggéré dans cette documentation metafor - Problèmes de convergence avec la fonction rma () . De plus, j'ai essayé d'autres algorithmes de convergence comme bobyqaet newuoacomme suggéré dans la documentation d' aide de la fonction rma.mv , mais j'ai obtenu la même réponse bimodale.
J'ai également essayé d'éliminer certaines des valeurs aberrantes potentielles du groupe problématique, comme suggéré dans Comment interpréter la distribution multimodale de la corrélation amorcée , mais en vain.
Je ne pouvais pas trouver un moyen de reproduire cela, j'ai donc téléchargé des données sur un référentiel GitHub (également les liens dans la section de code ci-dessous devraient charger dans votre environnement tout ce qui est nécessaire pour tester le cas). J'exécute le bootstrapping sur un cluster Linux en tant que travail de tableau (juste au cas où le script shell est job.sh , qui exécute sur chaque CPU le script R bootstrap.r qui exécute le modèle décrit ci-dessous). Un seul passage prend 2-3 minutes. Notez que l'amorçage 100 fois est également suffisant pour détecter la réponse bimodale. Voici un exemple pour 1000 itérations. Je connais R et d'autres méthodes mais pas tant que ça avec la méta-analyse.
J'apprécierais de l'aide pour comprendre si la distribution bimodale est correcte (bien que cela puisse être dû à des problèmes de convergence) et sinon, que peut-on faire à ce sujet? (en plus de ce que j'ai déjà essayé)
Ci-dessous - comparaison des coefficients du bootstrap (lignes rouges) et d'une seule exécution complète du modèle (lignes bleues). Les histogrammes illustrent les distributions bootstrapées pour chaque coefficient. L'échantillonnage des données pour le bootstrap a été effectué en sélectionnant avec remplacement de chaque groupe / combinaison formé par les deux effets fixes. Leurs tailles d'échantillon brut sont:
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)
Créé le 2019-05-02 par le package reprex (v0.2.1)
Le modèle va comme ceci:
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))
Informations sur la session 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)
la source
Sans avoir accès à un exemple reproductible, il est extrêmement difficile de donner une réponse définitive à ce comportement d'amorçage. En supposant qu'il n'y a effectivement pas de valeurs aberrantes, je soupçonne que nous observons un cas bénin du phénomène de Stein, d' autant plus qu'une méthodologie à effets mixtes nous suggère un certain regroupement dans nos données.
Cela dit, je suggère d'aller de l'avant et d'examiner certaines des séries à partir des valeurs d'1 2.4
f2f2_3:f1f1_2interaction «inhabituelles» , où il existe des valeurs très différentes, et d'étudier la distribution marginale de ces deux sous-échantillons aléatoires. Par exemple, dans certains cas,f2f2_3:f1f1_2est bien inférieur à tandis que le modèle estimé suggère des valeurs proches de . La distribution marginale est-elle similaire? Y a-t-il un cas de chevauchement insuffisant? Peut-être que le bootstrap "simple" est inapproprié et nous devons stratifier l'échantillon à portée de main en fonction de certains facteurs.la source