Saya membantu seorang rekan untuk mem - bootstrap model efek-campuran meta-analisis menggunakan kerangka paket metafor R yang ditulis oleh @Wolfgang.
Menarik dan mengkhawatirkan, untuk salah satu koefisien model saya mendapatkan distribusi bimodal ketika bootstrap (lihat panel kanan bawah pada gambar di bawah).
Saya kira salah satu penyebab utama bisa menjadi kenyataan bahwa ketika bootstrap, katakanlah setengah dari model berkumpul di solusi lokal dan setengah lainnya di yang lain. Saya mencoba untuk menyetel algoritma konvergensi seperti yang disarankan dalam dokumentasi metafor ini - Masalah Konvergensi dengan fungsi rma () . Juga, saya mencoba algoritma konvergensi lain seperti bobyqa
dan newuoa
seperti yang disarankan dalam dokumentasi bantuan fungsi rma.mv , tetapi mendapat respons bimodal yang sama.
Saya juga mencoba untuk menghilangkan beberapa outlier potensial dari kelompok bermasalah seperti yang disarankan dalam Bagaimana menafsirkan distribusi multimodal dari korelasi bootstrap , tetapi tidak berhasil.
Saya tidak dapat menemukan cara untuk mereproduksi ini jadi saya mengunggah data pada repositori GitHub (juga tautan di bagian kode di bawah ini akan memuat di lingkungan Anda semua yang diperlukan untuk menguji kasus ini). Saya menjalankan bootstrap pada Linux cluster sebagai pekerjaan array (untuk berjaga-jaga, skrip shell adalah job.sh , yang mengeksekusi pada setiap CPU script R bootstrap.r yang menjalankan model yang dijelaskan di bawah). Lari tunggal membutuhkan waktu 2-3 menit. Perhatikan bahwa bootstrap 100 kali juga cukup untuk mendeteksi respons bimodal. Di bawah ini adalah contoh untuk 1000 iterasi. Saya akrab dengan R dan metode lain tetapi tidak banyak dengan meta-analisis.
Saya akan sangat menghargai bantuan dengan pemahaman jika distribusi bimodal baik-baik saja (meskipun mungkin karena masalah konvergensi) dan jika tidak, lalu apa yang dapat dilakukan? (Selain apa yang sudah saya coba)
Bawah - membandingkan koefisien dari bootstrap (garis merah) dan dari menjalankan model penuh tunggal (garis biru). Histogram menggambarkan distribusi bootstrap untuk setiap koefisien. Pengambilan sampel data untuk bootstrap dilakukan dengan memilih dengan penggantian dari setiap kelompok / kombinasi yang dibentuk oleh dua efek tetap. Ukuran sampel mentah mereka adalah:
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)
Dibuat pada 2019-05-02 oleh paket reprex (v0.2.1)
Modelnya seperti ini:
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))
Info sesi 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)
sumber
Tanpa memiliki akses ke contoh yang dapat direproduksi sangat sulit untuk memberikan jawaban yang pasti untuk perilaku bootstrap ini. Dengan asumsi bahwa memang tidak ada outlier, saya menduga bahwa kami mengamati kasus ringan dari fenomena Stein terutama sebagai metodologi efek campuran menunjukkan kita beberapa pengelompokan dalam data kami.
Setelah mengatakan hal di atas, saya akan menyarankan untuk terus maju dan melihat beberapa proses dari nilai-nilai1 2.4
f2f2_3:f1f1_2
interaksi "tidak biasa" , di mana ada nilai yang sangat berbeda, dan menyelidiki distribusi marjinal dari dua subsampel acak ini. Misalnya dalam beberapa kasus,f2f2_3:f1f1_2
jauh di bawah sedangkan model perkiraan menyarankan nilai mendekati . Apakah distribusi marginalnya serupa? Apakah ada kasus tumpang tindih yang tidak memadai? Mungkin bootstrap "sederhana" tidak pantas dan kita perlu stratifikasi sampel yang ada sehubungan dengan beberapa faktor.sumber