2) Moderator analysis of PN and VOE using meta-analysis
1 Document Overview:
In this document, using condition-level data, we:
- Find moderators of PN and VOE in experimental conditions. This includes:
finding the moderators for PN and VOE separately
for significant moderators, estimating both effects within the same additive and interaction models
comparing additive and interaction models to see which better fits the data
- Find moderators of PN and VOE in both experimental and control conditions. This includes:
finding the moderators for PN and VOE separately
for exp_vs_control, estimating both effects within the same additive and interaction models
comparing additive and interaction models to see which better fits the data
2 Imports
options(scipen = 1, digits = 3)
library(pacman)
pacman::p_load(plyr,
tidyverse,
metafor,
lme4,
lmerTest,
cowplot,
easystats,
sjPlot,
kableExtra,
easystats,
MASS,
performance,
effects,
patchwork,
viridis,
janitor)
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.6.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] janitor_2.2.0 viridis_0.6.4 viridisLite_0.4.2
## [4] patchwork_1.1.3 effects_4.2-2 carData_3.0-5
## [7] MASS_7.3-60 kableExtra_1.3.4 sjPlot_2.8.15
## [10] see_0.8.0 report_0.5.7 parameters_0.21.2
## [13] performance_0.10.8 modelbased_0.8.6 insight_0.19.6
## [16] effectsize_0.8.6 datawizard_0.9.0 correlation_0.8.4
## [19] bayestestR_0.13.1 easystats_0.6.0 cowplot_1.1.1
## [22] lmerTest_3.1-3 lme4_1.1-34 metafor_4.4-0
## [25] numDeriv_2016.8-1.1 metadat_1.2-0 Matrix_1.6-1.1
## [28] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.0
## [31] dplyr_1.1.3 purrr_1.0.2 readr_2.1.4
## [34] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4
## [37] tidyverse_2.0.0 plyr_1.8.9 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] DBI_1.1.3 gridExtra_2.3 sandwich_3.0-2 rlang_1.1.1
## [5] magrittr_2.0.3 multcomp_1.4-25 snakecase_0.11.1 compiler_4.3.2
## [9] systemfonts_1.0.5 vctrs_0.6.4 rvest_1.0.3 pkgconfig_2.0.3
## [13] fastmap_1.1.1 backports_1.4.1 utf8_1.2.4 rmarkdown_2.25
## [17] tzdb_0.4.0 nloptr_2.0.3 xfun_0.40 cachem_1.0.8
## [21] jsonlite_1.8.7 sjmisc_2.8.9 ggeffects_1.3.2 broom_1.0.5
## [25] R6_2.5.1 bslib_0.5.1 stringi_1.7.12 boot_1.3-28.1
## [29] jquerylib_0.1.4 estimability_1.4.1 Rcpp_1.0.11 knitr_1.45
## [33] modelr_0.1.11 zoo_1.8-12 splines_4.3.2 nnet_7.3-19
## [37] timechange_0.2.0 tidyselect_1.2.0 rstudioapi_0.15.0 yaml_2.3.7
## [41] codetools_0.2-19 sjlabelled_1.2.0 lattice_0.21-9 withr_2.5.2
## [45] evaluate_0.22 survival_3.5-7 survey_4.2-1 xml2_1.3.5
## [49] pillar_1.9.0 generics_0.1.3 mathjaxr_1.6-0 hms_1.1.3
## [53] munsell_0.5.0 scales_1.2.1 minqa_1.2.6 glue_1.6.2
## [57] emmeans_1.8.9 tools_4.3.2 webshot_0.5.5 mvtnorm_1.2-3
## [61] grid_4.3.2 mitools_2.4 colorspace_2.1-0 nlme_3.1-163
## [65] cli_3.6.1 fansi_1.0.5 svglite_2.1.2 sjstats_0.18.2
## [69] gtable_0.3.4 sass_0.4.7 digest_0.6.33 TH.data_1.1-2
## [73] htmltools_0.5.6.1 lifecycle_1.0.3 httr_1.4.7i2_from_q_df <- function(Q, df) {
return ((Q-df)/Q*100)
}3 Read Meta-analysis data
meta_analysis_data <- read_csv('processed_data/meta-analysis_data.csv') %>%
filter(study_ID != "sanal-hayes_2022") %>% #comment out to replicate results with Sanal-Hayes (2022) found in SM
mutate(age_z = scale(mean_age_1))## Rows: 109 Columns: 66
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (34): study_ID, long_cite, short_cite, doi, peer_reviewed, coder, expt_c...
## dbl (18): expt_num, n_1, mean_age_1, sd_age_1, x_1, x_2, x_3, SD_1, SD_2, SD...
## lgl (14): group_name_1, group_name_2, n_2, mean_age_2, t, F, r, d, d_var, co...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.4 Moderator analysis for experimental conditions
meta_analysis_data_exp <- meta_analysis_data %>%
filter(exp_or_control == "experimental")4.1 Modeling each effect separately
4.1.1 PN Effect
#Prepare data and contrasts
meta_analysis_data_exp$equal_per_nov <- as.factor(meta_analysis_data_exp$equal_per_nov)
meta_analysis_data_exp$exposure_phase <- as.factor(meta_analysis_data_exp$exposure_phase)
meta_analysis_data_exp$domain <- as.factor(meta_analysis_data_exp$domain)
meta_analysis_data_exp$stim_loop <- as.factor(meta_analysis_data_exp$stim_loop)
meta_analysis_data_exp$equal_per_nov <- relevel(meta_analysis_data_exp$equal_per_nov, "yes")
meta_analysis_data_exp$exposure_phase <- relevel(meta_analysis_data_exp$exposure_phase, "habituation")
meta_analysis_data_exp$domain <- relevel(meta_analysis_data_exp$domain, "psychology")
meta_analysis_data_exp$stim_loop <- relevel(meta_analysis_data_exp$stim_loop, "yes")
contrasts(meta_analysis_data_exp$equal_per_nov) = contr.sum(2)
contrasts(meta_analysis_data_exp$exposure_phase) = contr.sum(2)
contrasts(meta_analysis_data_exp$domain) = contr.sum(2)
contrasts(meta_analysis_data_exp$stim_loop) = contr.sum(2)#PN model
smd_pn <- metafor::escalc(measure = "SMD", m1i = x_2, sd1i = SD_2, n1i = n_1, m2i = x_1, sd2i = SD_1, n2i = n_1, data = meta_analysis_data_exp)
rma_pn <- rma.mv(yi, vi, mods = ~equal_per_nov + exposure_phase + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = smd_pn)## Warning: 1 row with NAs omitted from model fitting.rma_pn##
## Multivariate Meta-Analysis Model (k = 75; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1394 0.3733 75 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 69) = 176.4956, p-val < .0001
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 36.3222, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.2998 0.0787 3.8094 0.0001 0.1456 0.4541 ***
## equal_per_nov1 0.0177 0.0738 0.2402 0.8101 -0.1269 0.1623
## exposure_phase1 0.3403 0.0593 5.7407 <.0001 0.2241 0.4565 ***
## age_z 0.0867 0.0617 1.4051 0.1600 -0.0342 0.2075
## domain1 -0.0615 0.0927 -0.6631 0.5073 -0.2432 0.1202
## stim_loop1 0.1029 0.0768 1.3405 0.1801 -0.0476 0.2533
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_from_q_df(rma_pn$QE, rma_pn$QEdf)## [1] 60.9#PN Table
rma_pn_results <-data.frame(b = round(rma_pn$b, 3), SE = round(rma_pn$se, 3), z = round(rma_pn$zval, 3), p = ifelse(rma_pn$pval <= .001, "<.001", round(rma_pn$pval, 3)), CI = paste(round(rma_pn$ci.lb, 3), round(rma_pn$ci.ub, 3), sep = " - "))
rma_pn_results1 <- rma_pn_results
rownames(rma_pn_results) <- c("Intercept", "Equal PN: yes", "Exposure Phase: habituation", "Age (z-score)", "Domain: psychology", "Stimulus Loop: yes")
colnames(rma_pn_results) <- c("Estimate", "Standard Error", "z", "p", "95% CI")
rma_pn_results%>%
kable("html",align = 'clc',table.attr="id=rmapntable")%>%
kable_styling("striped")| Estimate | Standard Error | z | p | 95% CI | |
|---|---|---|---|---|---|
| Intercept | 0.300 | 0.079 | 3.809 | <.001 | 0.146 - 0.454 |
| Equal PN: yes | 0.018 | 0.074 | 0.240 | 0.81 | -0.127 - 0.162 |
| Exposure Phase: habituation | 0.340 | 0.059 | 5.741 | <.001 | 0.224 - 0.457 |
| Age (z-score) | 0.087 | 0.062 | 1.405 | 0.16 | -0.034 - 0.208 |
| Domain: psychology | -0.061 | 0.093 | -0.663 | 0.507 | -0.243 - 0.12 |
| Stimulus Loop: yes | 0.103 | 0.077 | 1.340 | 0.18 | -0.048 - 0.253 |
#PN Visualize coefficients
PN_mod_plot <- plot_model(rma_pn, terms = c("Overall", "equal_per_nov1", "exposure_phase1", "age_z", "domain1", "stim_loop1"), ci.lvl = 0.95, show.values = TRUE, show.p = TRUE, colors = "bw") ## Warning: 1 row with NAs omitted from model fitting.4.1.2 VOE Effect
#VOE model
smd_voe <- metafor::escalc(measure = "SMD", m1i = x_3, sd1i = SD_3, n1i = n_1, m2i = x_2, sd2i = SD_2, n2i = n_1, data = meta_analysis_data_exp)
rma_voe <- rma.mv(yi, vi, mods = ~equal_per_nov + exposure_phase + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = smd_voe)## Warning: 1 row with NAs omitted from model fitting.rma_voe##
## Multivariate Meta-Analysis Model (k = 75; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0406 0.2015 75 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 69) = 107.2781, p-val = 0.0022
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 7.3429, p-val = 0.1964
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.2621 0.0577 4.5387 <.0001 0.1489 0.3752 ***
## equal_per_nov1 -0.0301 0.0536 -0.5614 0.5745 -0.1351 0.0750
## exposure_phase1 -0.0563 0.0445 -1.2654 0.2057 -0.1434 0.0309
## age_z -0.1062 0.0467 -2.2735 0.0230 -0.1977 -0.0146 *
## domain1 0.0613 0.0668 0.9186 0.3583 -0.0695 0.1922
## stim_loop1 0.0025 0.0549 0.0454 0.9638 -0.1051 0.1101
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_from_q_df(rma_voe$QE, rma_voe$QEdf)## [1] 35.7#VOE Table
rma_voe_results <-data.frame(b = round(rma_voe$b, 3), SE = round(rma_voe$se, 3), z = round(rma_voe$zval, 3), p = ifelse(rma_voe$pval <= .001, "<.001", round(rma_voe$pval, 3)), CI = paste(round(rma_voe$ci.lb, 3), round(rma_voe$ci.ub, 3), sep = " - "))
rma_voe_results1 <- rma_voe_results
rownames(rma_voe_results) <- c("Intercept", "Equal PN: yes", "Exposure Phase: habituation", "Age (z-score)", "Domain: psychology", "Stimulus Loop: yes")
colnames(rma_voe_results) <- c("Estimate", "Standard Error", "z", "p", "95% CI")
rma_voe_results%>%
#kbl() %>%
kable("html",align = 'clc',table.attr="id=rmavoetable")%>%
kable_styling("striped")| Estimate | Standard Error | z | p | 95% CI | |
|---|---|---|---|---|---|
| Intercept | 0.262 | 0.058 | 4.539 | <.001 | 0.149 - 0.375 |
| Equal PN: yes | -0.030 | 0.054 | -0.561 | 0.575 | -0.135 - 0.075 |
| Exposure Phase: habituation | -0.056 | 0.044 | -1.265 | 0.206 | -0.143 - 0.031 |
| Age (z-score) | -0.106 | 0.047 | -2.273 | 0.023 | -0.198 - -0.015 |
| Domain: psychology | 0.061 | 0.067 | 0.919 | 0.358 | -0.07 - 0.192 |
| Stimulus Loop: yes | 0.002 | 0.055 | 0.045 | 0.964 | -0.105 - 0.11 |
#VOE Visualize coefficients
VOE_mod_plot <- plot_model(rma_voe, terms = c("Overall", "equal_per_nov1", "exposure_phase1", "age_z", "domain1", "stim_loop1"), ci.lvl = 0.95, show.values = TRUE, show.p = TRUE, colors = "bw")## Warning: 1 row with NAs omitted from model fitting.4.1.3 Visualize PN and VOE effects
PN_mod_plot + VOE_mod_plot
Modeling each effect separately, we found that PN and VOE were moderated by different predictors. For PN, the only significant moderator was whether the experiment used habituation or familiarization. Exposing infants to the training stimulus until the infant habituated, rather than familiarizing them for a fixed number of trials, evoked a greater PN effect (estimate = 0.34, SE = 0.059, [0.224 - 0.457], z = 5.741, p = <.001). By contrast, habituation versus familiarization did not moderate the size of the VOE effect (estimate = -0.056, SE = 0.044, [-0.143 - 0.031], z = -1.265, p = 0.206).
For the VOE effect, infant age was the only significant moderator. Studies that tested older infants, relative to younger infants, reported a smaller VOE effect (estimate = -0.106, SE = 0.047, [-0.198 - -0.015], z =. -2.273, p = 0.023). Infant age did not predict the size of the PN effect (estimate = 0.087, SE = 0.062, [-0.034 - 0.208], z = 1.405, p = 0.16). This analysis also revealed some potentially interesting negative results. We found that the VOE effect is comparable in size across the domains of intuitive psychology and physics (estimate = 0.061, SE = 0.067, [-0.07 - 0.192], z = 0.919, p = 0.358), and regardless of whether it was made to compete against an expected test event that was matched in novelty, or was more perceptually novel, than the unexpected test event (estimate = -0.03, SE = 0.054, [-0.135 - 0.075], z = -0.561, p = 0.575).
4.1.4 Boxplot Exposure Phase
PN_boxplot_exposure_phase <- ggplot(smd_pn, aes(x = "PN", y = yi)) +
geom_boxplot() +
geom_point(aes(size = n_1, colour = exposure_phase), alpha = .2) +
ylab("SMD")+
xlab("PN") +
ylim(-2, 2) +
scale_size(range = c(1, 3)) +
scale_x_discrete(labels = NULL) +
facet_wrap(~exposure_phase) +
ggtitle("PN") +
theme(plot.title = element_text(hjust = 0.5, size = 22), axis.title = element_text(size = 16), axis.text = element_text(size = 12))VOE_boxplot_exposure_phase <- ggplot(smd_voe, aes(x = "VOE", y = yi)) +
geom_boxplot() +
geom_point(aes(size = n_1, colour = exposure_phase), alpha = .2) +
ylab("SMD")+
xlab("VOE") +
ylim(-2, 2) +
scale_size(range = c(1, 3)) +
scale_x_discrete(labels = NULL) +
facet_wrap(~exposure_phase) +
ggtitle("VOE") +
theme(plot.title = element_text(hjust = 0.5, size = 22), axis.title = element_text(size = 16), axis.text = element_text(size = 12))PN_boxplot_exposure_phase + VOE_boxplot_exposure_phase## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
4.1.5 Scatterplot Age
PN_scatter_age <- ggplot(data = smd_pn, aes(mean_age_1, yi)) +
geom_point(aes(size = n_1, colour = mean_age_1)) +
scale_color_viridis() +
scale_size(range = c(1, 3)) +
geom_smooth(method = "lm") +
geom_errorbar(width = .1, alpha = 1, aes(ymin = yi - vi, ymax = yi + vi)) + #error bars are sampling variances here
ggtitle("PN") +
ylab("SMD") +
xlab("Mean age") +
coord_cartesian(ylim = c(-2,2)) +
theme(plot.title = element_text(hjust = 0.5, size = 22), axis.title = element_text(size = 16), axis.text = element_text(size = 12))VOE_scatter_age <- ggplot(data = smd_voe, aes(mean_age_1, yi)) +
geom_point(aes(size = n_1, colour = mean_age_1)) +
scale_color_viridis() +
scale_size(range = c(1, 3)) +
geom_smooth(method = "lm") +
geom_errorbar(width = .1, alpha = 1, aes(ymin = yi - vi, ymax = yi + vi)) + #error bars are sampling variances here
ggtitle("VOE") +
ylab("SMD") +
xlab("Mean age") +
coord_cartesian(ylim = c(-2,2)) +
theme(plot.title = element_text(hjust = 0.5, size = 22), axis.title = element_text(size = 16), axis.text = element_text(size = 12))PN_scatter_age + VOE_scatter_age## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Removed 1 rows containing missing values (`geom_point()`).
4.2 Estimating both effects within the same model to address statistical non-independence
4.2.1 Exposure Phase
#Prepare data and contrasts
meta_analysis_data_long <- meta_analysis_data_exp %>%
pivot_longer(cols = c(x_1, x_2, x_3, SD_1, SD_2, SD_3), names_to = c(".value", "trial_type"), names_sep = "_") %>%
mutate(trial_type = ifelse(trial_type == 1, "last_train", ifelse(trial_type == 2, "expected", "unexpected"))) %>%
mutate(variance = SD^2/n_1)
meta_analysis_data_long$trial_type <- as.factor(meta_analysis_data_long$trial_type)
contrasts(meta_analysis_data_long$exposure_phase) = contr.treatment(2)
contrasts(meta_analysis_data_long$trial_type) = contr.treatment(3)#Additive model
results_exposure_phase_add <- rma.mv(x, variance, mods = ~trial_type + exposure_phase + equal_per_nov + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = meta_analysis_data_long)## Warning: 3 rows with NAs omitted from model fitting.results_exposure_phase_add##
## Multivariate Meta-Analysis Model (k = 225; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 28.3510 5.3246 75 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 217) = 1891.5818, p-val < .0001
##
## Test of Moderators (coefficients 2:8):
## QM(df = 7) = 350.8220, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 13.0155 1.1860 10.9742 <.0001 10.6910 15.3401 ***
## trial_type2 -1.7214 0.2226 -7.7333 <.0001 -2.1577 -1.2852 ***
## trial_type3 2.0872 0.2500 8.3488 <.0001 1.5972 2.5771 ***
## exposure_phase2 5.9509 1.3259 4.4883 <.0001 3.3523 8.5496 ***
## equal_per_nov1 1.1040 0.8680 1.2719 0.2034 -0.5973 2.8052
## age_z -3.5009 0.6864 -5.1007 <.0001 -4.8462 -2.1557 ***
## domain1 0.1437 1.0758 0.1336 0.8937 -1.9649 2.2523
## stim_loop1 4.0095 0.9151 4.3816 <.0001 2.2160 5.8030 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Interaction model
results_exposure_phase_inter <- rma.mv(x, variance, mods = ~trial_type * exposure_phase + equal_per_nov + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = meta_analysis_data_long)## Warning: 3 rows with NAs omitted from model fitting.results_exposure_phase_inter##
## Multivariate Meta-Analysis Model (k = 225; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 27.9988 5.2914 75 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 215) = 1723.1646, p-val < .0001
##
## Test of Moderators (coefficients 2:10):
## QM(df = 9) = 487.7732, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt 13.9831 1.1868 11.7825 <.0001 11.6571
## trial_type2 -3.4310 0.2906 -11.8067 <.0001 -4.0006
## trial_type3 2.1713 0.3576 6.0719 <.0001 1.4704
## exposure_phase2 4.1744 1.3436 3.1068 0.0019 1.5410
## equal_per_nov1 1.1013 0.8629 1.2763 0.2018 -0.5899
## age_z -3.4638 0.6823 -5.0765 <.0001 -4.8012
## domain1 0.0907 1.0696 0.0848 0.9325 -2.0058
## stim_loop1 4.0188 0.9097 4.4178 <.0001 2.2359
## trial_type2:exposure_phase2 4.6433 0.4574 10.1510 <.0001 3.7467
## trial_type3:exposure_phase2 -0.2584 0.5002 -0.5166 0.6054 -1.2387
## ci.ub
## intrcpt 16.3091 ***
## trial_type2 -2.8614 ***
## trial_type3 2.8722 ***
## exposure_phase2 6.8079 **
## equal_per_nov1 2.7925
## age_z -2.1265 ***
## domain1 2.1871
## stim_loop1 5.8018 ***
## trial_type2:exposure_phase2 5.5398 ***
## trial_type3:exposure_phase2 0.7219
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Does the interaction model better explain the data
anova_test_exposure_phase <- anova(results_exposure_phase_inter, results_exposure_phase_add, refit = TRUE)## Warning: 3 rows with NAs omitted from model fitting.
## Warning: 3 rows with NAs omitted from model fitting.anova_test_exposure_phase##
## df AIC BIC AICc logLik LRT pval QE
## Full 11 1550.4591 1588.0362 1551.6986 -764.2296 1723.1646
## Reduced 9 1682.3438 1713.0887 1683.1810 -832.1719 135.8847 <.0001 1891.5818#Visualize data
ggplot(data = meta_analysis_data_long, aes(x = factor(trial_type, level=c('last_train', 'expected', 'unexpected')), y = x)) +
geom_boxplot() +
xlab("Trial Type")+
ylab("Looking Time") +
geom_point(aes(size = n_1, colour = exposure_phase), alpha = .2) +
facet_wrap(~exposure_phase) +
theme(plot.title = element_text(hjust = 0.5, size = 22), axis.title = element_text(size = 18), axis.text = element_text(size = 16), strip.text = element_text(size = 16), legend.text = element_text(size = 15), legend.title = element_text(size = 16))## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).## Warning: Removed 3 rows containing missing values (`geom_point()`).
We tested whether the two significant predictors that emerged impacted the two effects differently when VOE and PN were estimated in the same model, including an interaction between trial type and exposure phase, and the other moderators as fixed additive effects. For both effects, we found a significant interaction between each predictor and trial type (exposure phase X trial type: Likelihood Ratio Test, p <.001). Contrasts extracted from these models confirmed the findings that familiarization predicted a smaller PN effect than habituation (estimate = -4.643, SE = 0.457, [-5.54, -3.747], z = -10.151, p <.001) but did not affect the VOE effect (estimate = -0.258, SE = 0.5, [-1.239, 0.722], z = -0.517, p = 0.605).
4.2.2 Age
contrasts(meta_analysis_data_long$exposure_phase) = contr.sum(2)
#Additive model
results_age_add <- rma.mv(x, variance, mods = ~trial_type + age_z + exposure_phase + equal_per_nov + domain + stim_loop, random = ~1|paper_expt_info, data = meta_analysis_data_long)## Warning: 3 rows with NAs omitted from model fitting.results_age_add##
## Multivariate Meta-Analysis Model (k = 225; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 28.3510 5.3246 75 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 217) = 1891.5818, p-val < .0001
##
## Test of Moderators (coefficients 2:8):
## QM(df = 7) = 350.8220, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 15.9910 0.9048 17.6733 <.0001 14.2176 17.7644 ***
## trial_type2 -1.7214 0.2226 -7.7333 <.0001 -2.1577 -1.2852 ***
## trial_type3 2.0872 0.2500 8.3488 <.0001 1.5972 2.5771 ***
## age_z -3.5009 0.6864 -5.1007 <.0001 -4.8462 -2.1557 ***
## exposure_phase1 -2.9755 0.6629 -4.4883 <.0001 -4.2748 -1.6761 ***
## equal_per_nov1 1.1040 0.8680 1.2719 0.2034 -0.5973 2.8052
## domain1 0.1437 1.0758 0.1336 0.8937 -1.9649 2.2523
## stim_loop1 4.0095 0.9151 4.3816 <.0001 2.2160 5.8030 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Interaction model
results_age_inter <- rma.mv(x, variance, mods = ~trial_type * age_z + exposure_phase + equal_per_nov + domain + stim_loop, random = ~1|paper_expt_info, data = meta_analysis_data_long)## Warning: 3 rows with NAs omitted from model fitting.results_age_inter##
## Multivariate Meta-Analysis Model (k = 225; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 27.5792 5.2516 75 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 215) = 1850.4608, p-val < .0001
##
## Test of Moderators (coefficients 2:10):
## QM(df = 9) = 372.3040, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 15.6245 0.8988 17.3833 <.0001 13.8628 17.3861 ***
## trial_type2 -1.4666 0.2643 -5.5495 <.0001 -1.9846 -0.9487 ***
## trial_type3 2.8896 0.3105 9.3059 <.0001 2.2810 3.4982 ***
## age_z -2.9752 0.6969 -4.2693 <.0001 -4.3411 -1.6094 ***
## exposure_phase1 -2.9862 0.6544 -4.5632 <.0001 -4.2689 -1.7036 ***
## equal_per_nov1 1.1075 0.8567 1.2927 0.1961 -0.5717 2.7867
## domain1 0.2349 1.0623 0.2211 0.8250 -1.8472 2.3170
## stim_loop1 3.9877 0.9032 4.4152 <.0001 2.2175 5.7579 ***
## trial_type2:age_z -0.4893 0.2612 -1.8736 0.0610 -1.0012 0.0226 .
## trial_type3:age_z -1.3885 0.3192 -4.3502 <.0001 -2.0140 -0.7629 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Does the interaction model better explain the data?
anova_test_age <- anova(results_age_inter, results_age_add, refit = TRUE)## Warning: 3 rows with NAs omitted from model fitting.
## Warning: 3 rows with NAs omitted from model fitting.anova_test_age##
## df AIC BIC AICc logLik LRT pval QE
## Full 11 1667.2435 1704.8206 1668.4829 -822.6218 1850.4608
## Reduced 9 1682.3438 1713.0887 1683.1810 -832.1719 19.1003 <.0001 1891.5818#Visualize data
ggplot(data = meta_analysis_data_long, aes(x = mean_age_1, y = x)) +
xlab("Mean Infant Age (days)")+
ylab("Looking Time") +
geom_point(aes(size = n_1, colour = mean_age_1), alpha = .3) +
facet_wrap(~factor(trial_type, level=c('last_train', 'expected', 'unexpected'))) +
geom_errorbar(data = meta_analysis_data_long, width = .1, alpha = 0.1, inherit.aes = FALSE, aes(x = mean_age_1, ymin = x - SD/sqrt(n_1)*1.96, ymax = x + SD/sqrt(n_1)*1.96))+
scale_color_viridis() +
geom_smooth(method = "lm") +
theme(plot.title = element_text(hjust = 0.5, size = 22), axis.title = element_text(size = 18), axis.text = element_text(size = 16), strip.text = element_text(size = 16), legend.text = element_text(size = 15), legend.title = element_text(size = 16))## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 3 rows containing missing values (`geom_point()`).
We tested whether the two significant predictors that emerged impacted the two effects differently when VOE and PN were estimated in the same model, including an interaction between trial type and age, and the other moderators as fixed additive effects. For both effects, we found a significant interaction between each predictor and trial type (age X trial type: Likelihood Ratio Test, p <.001). Contrasts extracted from these models confirmed the findings that infant age predicted the size of the VOE effect (estimate = -1.388, SE = 0.319, [-2.014, -0.763], z = -4.35, p <.001) but not the PN effect (estimate = 0.489, SE = 0.261, [-0.023, 1.001], z = 1.874, p = 0.061). The only difference we observed when modeling both effects together was that there was a marginal effect of age on the PN effect, such that older infants showed a bigger effect than younger infants (the opposite of the pattern observed for the VOE effect).
5 Experimental or control moderator analysis
#Prepare data and contrasts
meta_analysis_data$equal_per_nov <- as.factor(meta_analysis_data$equal_per_nov)
meta_analysis_data$exposure_phase <- as.factor(meta_analysis_data$exposure_phase)
meta_analysis_data$domain <- as.factor(meta_analysis_data$domain)
meta_analysis_data$stim_loop <- as.factor(meta_analysis_data$stim_loop)
meta_analysis_data$exp_or_control <- as.factor(meta_analysis_data$exp_or_control)
meta_analysis_data$equal_per_nov <- relevel(meta_analysis_data$equal_per_nov, "yes")
meta_analysis_data$exposure_phase <- relevel(meta_analysis_data$exposure_phase, "habituation")
meta_analysis_data$domain <- relevel(meta_analysis_data$domain, "psychology")
meta_analysis_data$stim_loop <- relevel(meta_analysis_data$stim_loop, "yes")
meta_analysis_data$exp_or_control <- relevel(meta_analysis_data$exp_or_control, "experimental")
contrasts(meta_analysis_data$equal_per_nov) = contr.sum(2)
contrasts(meta_analysis_data$exposure_phase) = contr.sum(2)
contrasts(meta_analysis_data$domain) = contr.sum(2)
contrasts(meta_analysis_data$stim_loop) = contr.sum(2)
contrasts(meta_analysis_data$exp_or_control) = contr.sum(2)5.1 Modeling each effect separately
5.1.1 PN Effect
smd_pn_exp_control <- metafor::escalc(measure = "SMD", m1i = x_2, sd1i = SD_2, n1i = n_1, m2i = x_1, sd2i = SD_1, n2i = n_1, data = meta_analysis_data)
rma_pn_exp_control <- rma.mv(yi, vi, mods = ~equal_per_nov + exposure_phase + exp_or_control + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = smd_pn_exp_control)## Warning: 1 row with NAs omitted from model fitting.rma_pn_exp_control##
## Multivariate Meta-Analysis Model (k = 106; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1548 0.3935 106 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 99) = 259.6836, p-val < .0001
##
## Test of Moderators (coefficients 2:7):
## QM(df = 6) = 24.3909, p-val = 0.0004
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3693 0.0828 4.4623 <.0001 0.2071 0.5315 ***
## equal_per_nov1 0.0228 0.0692 0.3291 0.7421 -0.1129 0.1585
## exposure_phase1 0.2472 0.0531 4.6585 <.0001 0.1432 0.3513 ***
## exp_or_control1 -0.0513 0.0578 -0.8869 0.3751 -0.1646 0.0620
## age_z 0.0334 0.0544 0.6148 0.5387 -0.0732 0.1401
## domain1 -0.0839 0.0912 -0.9195 0.3578 -0.2627 0.0949
## stim_loop1 0.0718 0.0704 1.0189 0.3082 -0.0663 0.2098
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.1.2 VOE effect
smd_voe_exp_control <- metafor::escalc(measure = "SMD", m1i = x_3, sd1i = SD_3, n1i = n_1, m2i = x_2, sd2i = SD_2, n2i = n_1, data = meta_analysis_data)
rma_voe_exp_control <- rma.mv(yi, vi, mods = ~equal_per_nov + exposure_phase + exp_or_control + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = smd_voe_exp_control)## Warning: 1 row with NAs omitted from model fitting.rma_voe_exp_control##
## Multivariate Meta-Analysis Model (k = 106; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0325 0.1802 106 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 99) = 138.6782, p-val = 0.0053
##
## Test of Moderators (coefficients 2:7):
## QM(df = 6) = 15.4989, p-val = 0.0167
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.1005 0.0575 1.7478 0.0805 -0.0122 0.2132 .
## equal_per_nov1 -0.0220 0.0471 -0.4667 0.6407 -0.1144 0.0704
## exposure_phase1 -0.0247 0.0378 -0.6551 0.5124 -0.0988 0.0493
## exp_or_control1 0.1612 0.0414 3.8916 <.0001 0.0800 0.2424 ***
## age_z -0.0338 0.0388 -0.8722 0.3831 -0.1098 0.0422
## domain1 0.0314 0.0614 0.5115 0.6090 -0.0889 0.1518
## stim_loop1 -0.0024 0.0473 -0.0515 0.9589 -0.0952 0.0903
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We found that whether the condition was experimental or control did not predict the size of the PN effect (estimate = 0.051, SE = 0.058, [-0.062, 0.165], z = 0.887, p = 0.375), but did predict the size of the VOE, such that experimental conditions have a larger VOE effect (estimate = 0.161, SE = 0.041, [0.08, 0.242], z = 3.892, p = <.001).
5.2 Estimating both effects within the same model to address statistical non-independence: trial type and exp vs control
#Prepare data and contrasts
meta_analysis_data_long_exp_control <- meta_analysis_data %>%
pivot_longer(cols = c(x_1, x_2, x_3, SD_1, SD_2, SD_3), names_to = c(".value", "trial_type"), names_sep = "_") %>%
mutate(trial_type = ifelse(trial_type == 1, "last_train", ifelse(trial_type == 2, "expected", "unexpected"))) %>%
mutate(variance = SD^2/n_1)
meta_analysis_data_long_exp_control$trial_type <- as.factor(meta_analysis_data_long_exp_control$trial_type)
meta_analysis_data_long_exp_control$exp_or_control <- as.factor(meta_analysis_data_long_exp_control$exp_or_control)
meta_analysis_data_long_exp_control$exp_or_control <- relevel(meta_analysis_data_long_exp_control$exp_or_control, "control")
contrasts(meta_analysis_data_long_exp_control$exp_or_control) = contr.treatment(2)
contrasts(meta_analysis_data_long_exp_control$trial_type) = contr.treatment(3)#Additive model
results_expcontrol_add <- rma.mv(x, variance, mods = ~trial_type + exp_or_control + equal_per_nov + exposure_phase + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = meta_analysis_data_long_exp_control)## Warning: 3 rows with NAs omitted from model fitting.results_expcontrol_add##
## Multivariate Meta-Analysis Model (k = 318; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 28.4731 5.3360 106 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 309) = 2740.3076, p-val < .0001
##
## Test of Moderators (coefficients 2:9):
## QM(df = 8) = 412.2206, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 17.5281 1.2896 13.5918 <.0001 15.0005 20.0556 ***
## trial_type2 -1.7431 0.1954 -8.9185 <.0001 -2.1261 -1.3600 ***
## trial_type3 1.6258 0.2142 7.5909 <.0001 1.2060 2.0455 ***
## exp_or_control2 -0.4563 1.2331 -0.3700 0.7114 -2.8730 1.9605
## equal_per_nov1 1.3878 0.7784 1.7830 0.0746 -0.1378 2.9134 .
## exposure_phase1 -2.9173 0.5685 -5.1317 <.0001 -4.0315 -1.8031 ***
## age_z -3.2281 0.5856 -5.5125 <.0001 -4.3759 -2.0804 ***
## domain1 -1.2881 1.0225 -1.2597 0.2078 -3.2922 0.7160
## stim_loop1 3.8054 0.8021 4.7444 <.0001 2.2334 5.3774 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Interaction model
results_expcontrol_inter <- rma.mv(x, variance, mods = ~trial_type * exp_or_control + equal_per_nov + exposure_phase + age_z + domain + stim_loop, random = ~1|paper_expt_info, data = meta_analysis_data_long_exp_control)## Warning: 3 rows with NAs omitted from model fitting.results_expcontrol_inter##
## Multivariate Meta-Analysis Model (k = 318; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 28.5250 5.3409 106 no paper_expt_info
##
## Test for Residual Heterogeneity:
## QE(df = 307) = 2714.8088, p-val < .0001
##
## Test of Moderators (coefficients 2:11):
## QM(df = 10) = 427.2133, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt 17.9143 1.3077 13.6991 <.0001 15.3513
## trial_type2 -1.7459 0.4089 -4.2700 <.0001 -2.5473
## trial_type3 0.4018 0.4158 0.9663 0.3339 -0.4132
## exp_or_control2 -0.9626 1.2652 -0.7608 0.4468 -3.4425
## equal_per_nov1 1.4124 0.7791 1.8130 0.0698 -0.1145
## exposure_phase1 -2.9233 0.5690 -5.1378 <.0001 -4.0385
## age_z -3.2239 0.5861 -5.5001 <.0001 -4.3727
## domain1 -1.2878 1.0234 -1.2584 0.2082 -3.2936
## stim_loop1 3.7785 0.8028 4.7066 <.0001 2.2050
## trial_type2:exp_or_control2 0.0236 0.4655 0.0507 0.9596 -0.8888
## trial_type3:exp_or_control2 1.6800 0.4852 3.4627 0.0005 0.7291
## ci.ub
## intrcpt 20.4774 ***
## trial_type2 -0.9445 ***
## trial_type3 1.2167
## exp_or_control2 1.5172
## equal_per_nov1 2.9394 .
## exposure_phase1 -1.8081 ***
## age_z -2.0751 ***
## domain1 0.7180
## stim_loop1 5.3520 ***
## trial_type2:exp_or_control2 0.9359
## trial_type3:exp_or_control2 2.6309 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Does the interaction model better explain the data
anova_test_exp_control <- anova(results_expcontrol_inter, results_expcontrol_add, refit = TRUE)## Warning: 3 rows with NAs omitted from model fitting.
## Warning: 3 rows with NAs omitted from model fitting.anova_test_exp_control##
## df AIC BIC AICc logLik LRT pval QE
## Full 12 2245.9349 2291.0795 2246.9579 -1110.9675 2714.8088
## Reduced 10 2257.1523 2294.7728 2257.8689 -1118.5761 15.2174 0.0005 2740.3076We compared the fit of the model with this interaction to a simpler model including all predictors as main effects, using the anova() function. The model that contained an interaction between trial_type and exp_or_control better explained the data (Likelihood Ratio Test, p <.001), and also provided a better compromise between parsimony and fit (AIC = 2245.935 vs 2257.152, BIC = 2291.08 vs 2294.773).
5.2.1 Visualizing data
ggplot(data = meta_analysis_data_long_exp_control, aes(x = factor(trial_type, level=c('last_train', 'expected', 'unexpected')), y = x)) +
xlab("Trial Type") +
ylab("Looking Time") +
geom_boxplot() +
geom_point(aes(size = n_1, colour = exp_or_control), alpha = .2) +
#scale_size(c(1,3)) +
facet_wrap(~exp_or_control)## Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).## Warning: Removed 3 rows containing missing values (`geom_point()`).