Meta-analysis vs. random replication: which is better?
Init
library(kirkegaard)## Loading required package: tidyverse## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.4 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Loading required package: weights## Loading required package: Hmisc## Loading required package: lattice## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, units## Loading required package: assertthat##
## Attaching package: 'assertthat'## The following object is masked from 'package:tibble':
##
## has_name## Loading required package: magrittr##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extract## Loading required package: psych##
## Attaching package: 'psych'## The following object is masked from 'package:Hmisc':
##
## describe## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## Loading required package: metafor## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack##
## Loading the 'metafor' package (version 3.0-2). For an
## introduction to the package please type: help(metafor)## Loading required package: rlang##
## Attaching package: 'rlang'## The following object is masked from 'package:magrittr':
##
## set_names## The following object is masked from 'package:assertthat':
##
## has_name## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
## splice##
## Attaching package: 'kirkegaard'## The following object is masked from 'package:rlang':
##
## is_logical## The following object is masked from 'package:psych':
##
## rescale## The following object is masked from 'package:assertthat':
##
## are_equal## The following objects are masked from 'package:purrr':
##
## is_logical, is_numeric## The following object is masked from 'package:base':
##
## +load_packages(
metafor,
readxl
)
theme_set(theme_bw())Data
Data from: https://osf.io/dhp63/ Paper: https://psycnet.apa.org/record/2019-19962-005
#read the values
d = read_excel("MoneyPrimingMetaAnalysis_Dataset_FINAL_correction2.xlsx",
skip = 1,
sheet = 2) %>%
df_legalize_names()## New names:
## * Notes -> Notes...16
## * Notes -> Notes...40
## * Notes -> Notes...44
## * `Effect size (D) (Reported)` -> `Effect size (D) (Reported)...54`
## * `Direction of Effect` -> `Direction of Effect...63`
## * ...#recode
d %<>% mutate(
prereg = case_when(
Preregistered==0 ~ FALSE,
Preregistered==1 ~ TRUE)
)
#subsets
d_unreg = d %>% filter(!prereg)
d_reg = d %>% filter(prereg)Analysis
#meta-analysis unregistered studies
meta_unreg = d_unreg %$% rma(yi = d, vi = var_of_d)
meta_unreg##
## Random-Effects Model (k = 236; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.1407 (SE = 0.0176)
## tau (square root of estimated tau^2 value): 0.3751
## I^2 (total heterogeneity / total variability): 80.95%
## H^2 (total variability / sampling variability): 5.25
##
## Test for Heterogeneity:
## Q(df = 235) = 1005.6633, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.3148 0.0287 10.9846 <.0001 0.2586 0.3709 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#registered
meta_reg = d_reg %$% rma(yi = d, vi = var_of_d)
meta_reg##
## Random-Effects Model (k = 51; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0020 (SE = 0.0035)
## tau (square root of estimated tau^2 value): 0.0451
## I^2 (total heterogeneity / total variability): 9.48%
## H^2 (total variability / sampling variability): 1.10
##
## Test for Heterogeneity:
## Q(df = 50) = 59.2982, p-val = 0.1727
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0171 0.0215 0.7946 0.4269 -0.0250 0.0592
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#how often is registered more accurate than meta of unregistered?
d_reg %<>% mutate(
d_error = d - meta_reg$b[1],
d_abs_error = abs(d_error),
)
#trivial results by comparing to own mean, but we see the SD is .19
d_reg$d_error %>% describe2()#how large is the error from unreg?
meta_unreg$b[1] - meta_reg$b[1]## [1] 0.2976858#plot comparison
GG_denhist(d_reg, "d_error", vline = NULL) +
geom_vline(data = tibble(d_error = meta_unreg$b[1] - meta_reg$b[1]), aes(xintercept = d_error), linetype = "dashed", color = "blue", size = 1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
GG_save("figs/compare.png")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.#which fraction are more wrong, directionally?
(d_reg$d_error > (meta_unreg$b[1] - meta_reg$b[1])) %>% table2()#absolute?
(d_reg$d_abs_error > abs(meta_unreg$b[1] - meta_reg$b[1])) %>% table2()#with moderator
meta_mod = d %$% rma(yi = d, vi = var_of_d, mods = ~ prereg)
meta_mod##
## Mixed-Effects Model (k = 287; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1073 (SE = 0.0127)
## tau (square root of estimated tau^2 value): 0.3276
## I^2 (residual heterogeneity / unaccounted variability): 78.46%
## H^2 (unaccounted variability / sampling variability): 4.64
## R^2 (amount of heterogeneity accounted for): 10.06%
##
## Test for Residual Heterogeneity:
## QE(df = 285) = 1064.9615, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 22.3629, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3087 0.0260 11.8894 <.0001 0.2578 0.3596 ***
## preregTRUE -0.2752 0.0582 -4.7289 <.0001 -0.3893 -0.1612 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_mod %>% plot()
#density plot
d %>% GG_denhist("d", group = "prereg") +
scale_fill_discrete("Registered?") +
xlab("Effect size (d)")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
GG_save("figs/dist.png")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.#boxplot
d %>%
ggplot(aes(prereg, d, fill = prereg)) +
geom_boxplot()
#histograms
d %>%
ggplot(aes(d, fill = prereg)) +
geom_histogram() +
facet_wrap("prereg", nrow = 2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.