Chapter 5 Method 3: Two-Stage Individual Participant Meta-Analysis
5.1 Analytic Plan
In the present study, we estimate associations between the Big Five personality traits and crystallized cognitive abilities using a two-stage individual participant data meta-analysis (also known as a coordinated data analysis when multiple analysts are used). The procedure is as follows:
5.1.1 1. Sample-Specific Statistical Modeling
Models were run separately for each sample, personality characteristic, outcome, covariate, and moderator combination. To do so, we wrote a series of functions in the R programming language (see online materials) to (1) set up and run the model and extract model coefficients, (2) extract simple-effects predictions for moderator models (i.e. predicted values across levels of the moderator values), and (3) extract effect size metrics for the later meta-analysis. The basic form of the model is as follows:
\(Y_{ij}=b_0+b_1\ast predictor_{ij}+\epsilon_{ij}\) \(\epsilon_{ij}\sim\mathcal{N}(0, \sigma^2)\),
where \(b_1\) represents the effect of personality predicting the outcome separately in each sample.
The modeling function used the base R lm() function to run the models and the tidy() function from the broom package to extract model coefficients and confidence intervals (CI). Inferences will be made based on the 95% confidence intervals. Effect sizes and their standard errors were extracted from the model summary(). Simple-effects predictions were calculated by providing the full range of personality levels (0-10) and average levels of the covariates from the data used to estimate the model as the “newdata” argument in the base R predict() function.
5.1.2 2. Results Pooling Using Meta-Analysis
Once all the models were run, we next combined the effects across samples for each personality characteristic (5), outcome (1), covariate (6; none, each alone, all covariates), and moderator (3) combination and estimated a meta-analysis model of the target effect. We chose to conduct a random effects meta-analysis, which assumes that all effect sizes are randomly drawn from a population of effect sizes. To do so, we constructed three helper functions to (1) set up and run the meta-analytic models, (2) extract the meta-analytic estimates, and (3) extract heterogeneity estimates. The functions and more detail on them can be found in the online materials. The meta-analytic model is as follows:
\(T_i=\mu+\zeta_i+\epsilon_i\)
\(\epsilon_i\sim \mathcal{N}\left(0,\sigma^2\right)\ \) \(\zeta_i\sim \mathcal{N}\left(0,\tau^2\right)\ \) \(Cov\left(\zeta_i,\epsilon_i\right)=0\),
where \(T_i\) is the sample-specific effect of sample \(i\), \(\mu\) is the overall meta-analytic estimate, \(\zeta_i\) is true sample variability of sample \(i\) from the overall estimate, and \(\epsilon_i\) is sampling error.
Random effects meta-analyses were estimated using the metafor package in R. Meta-analytic estimates and confidence intervals were extracted using the coef() function. Inferences were based on the 95% confidence intervals (CI). Heterogeneity estimates were directly extracted from the model. Simple effects predictions of participant-level moderators were calculated using the weighted average sample-level predictions across levels of the moderators.
5.2 Step 1: Combine Data
We again need to combine data. However, rather than combining data across studies, for the two-stage approach, we’ll be combining data within studies in order to run separate analyses for each before combining via meta-analytic tools.
loadRData <- function(fileName, type){
#loads an RData file, and returns it
path <- sprintf("%s/data/clean/%s_cleaned.RData", local_path, fileName)
print(path)
load(path)
get(ls()[grepl(type, ls())])
}
ipd3_meta_data <- tibble(
study = studies[!studies %in% c("CNLSY", "SLS")]
, data = map(str_to_lower(study), ~loadRData(., "combined"))
) %>% mutate(
data = map(data, ~(.) %>%
ungroup() %>%
mutate(SID = as.character(SID)))
, study = mapvalues(study, studies, studies_long)
) %>%
unnest(data) %>%
mutate(age = ifelse(is.na(age), p_year - yearBrth, age))## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/base-i_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/eas_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/gsoep_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/hilda_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/hrs_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/lasa_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/map_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/mars_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/octo-twin_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/ros_cleaned.RData"
## [1] "~/Documents/projects/data synthesis/crystallized/data/clean/satsa_cleaned.RData"## # A tibble: 119,597 × 13
## study Trait p_value p_year SID Outcome o_value o_year education gender SRhealth yearBrth age
## <chr> <chr> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BASE E 3.83 1990 10004 crystallized 2.17 2000 17 0 4 1918 72
## 2 BASE E 3.83 1990 10010 crystallized 3.08 1995 11 0 2 1911 79
## 3 BASE E 2.5 1990 10033 crystallized 0.909 1997 8 1 5 1910 80
## 4 BASE E 3 1990 10034 crystallized 6.54 1995 10 0 5 1914 76
## 5 BASE E 2.5 1990 10111 crystallized 6.54 1995 13 1 4 1901 89
## 6 BASE E 3.17 1990 10115 crystallized 10 1995 11 0 2 1918 72
## 7 BASE E 3 1990 10116 crystallized 3.85 1995 10 0 4 1917 73
## 8 BASE E 3 1990 10145 crystallized 6.36 1997 10 1 5 1897 93
## 9 BASE E 3.33 1990 10175 crystallized 2.31 1995 8 1 3 1911 79
## 10 BASE E 3.33 1990 10188 crystallized 5.77 2004 10 1 2 1903 87
## # ℹ 119,587 more rows5.2.1 Study-Level Moderators
save_fun <- function(d, trait, outcome){
save(d, file = sprintf("%s/data/two_stage/meta_data/%s_%s.RData", local_path, trait, outcome))
}
url <- "https://github.com/emoriebeck/data-synthesis-tutorial/raw/main/codebooks/crystallized_tables.xlsx"
destfile <- "tables.xlsx"
curl::curl_download(url, destfile)
ipd_metaMod_data <- readxl::read_xlsx(destfile, sheet = "Table 4") %>%
select(-Category, -Construct, -category) %>%
pivot_longer(cols = c("BASE-I":"SATSA")
, names_to = "study"
, values_to = "value") %>%
pivot_wider(names_from = "name"
, values_from = "value") %>%
mutate(continent = relevel(factor(continent), ref = "North America")
, country = relevel(factor(country), ref = "United States")
, scale = relevel(factor(scale), ref = "NEO-FFI")) %>%
right_join(
ipd3_meta_data %>%
select(-p_value, -o_value, -(education:yearBrth))
) %>%
group_by(study, Trait, Outcome) %>%
mutate(baseAge = mean(age, na.rm = T) - 60, # center at age 60
predInt = mean(o_year - p_year) - 5, # center at 5 years
baseYear = ifelse(study %in% c("MARS", "MAP", "ROS"), mean(p_year), mean(p_year) - 2000)) %>% # center at 2000
ungroup() %>%
select(Trait, Outcome, study, continent, country, scale, baseAge, baseYear, predInt) %>%
distinct() %>%
group_by(Trait, Outcome) %>%
nest() %>%
ungroup()
ipd_metaMod_data## # A tibble: 5 × 3
## Trait Outcome data
## <chr> <chr> <list>
## 1 A crystallized <tibble [6 × 7]>
## 2 C crystallized <tibble [7 × 7]>
## 3 E crystallized <tibble [9 × 7]>
## 4 N crystallized <tibble [11 × 7]>
## 5 O crystallized <tibble [7 × 7]>## # A tibble: 5 × 4
## Trait Outcome data `pmap(list(data, Trait, Outcome), save_fun)`
## <chr> <chr> <list> <list>
## 1 A crystallized <tibble [6 × 7]> <NULL>
## 2 C crystallized <tibble [7 × 7]> <NULL>
## 3 E crystallized <tibble [9 × 7]> <NULL>
## 4 N crystallized <tibble [11 × 7]> <NULL>
## 5 O crystallized <tibble [7 × 7]> <NULL>readxl::read_xlsx(destfile, sheet = "Table 4") %>%
select(-Category, -Construct, -category) %>%
pivot_longer(cols = c("BASE-I":"SATSA")
, names_to = "study"
, values_to = "value") %>%
pivot_wider(names_from = "name"
, values_from = "value") %>%
right_join(
ipd3_meta_data %>%
select(-p_value, -o_value, -(education:yearBrth))
) %>%
group_by(study, Trait, Outcome) %>%
mutate(baseAge = mean(age, na.rm = T), # center at age 60
baseAgeSD = sd(age, na.rm = T), # center at age 60
predInt = mean(o_year - p_year), # center at 5 years
predIntSD = sd(o_year - p_year), # center at 5 years
baseYear = ifelse(study %in% c("MARS", "MAP", "ROS"), median(p_year), median(p_year)),
baseYearSD = ifelse(study %in% c("MARS", "MAP", "ROS"), sd(p_year), sd(p_year)),
cogYear = ifelse(study %in% c("MARS", "MAP", "ROS"), median(o_year), median(o_year)),
cogYearSD = ifelse(study %in% c("MARS", "MAP", "ROS"), sd(o_year), sd(o_year))) %>% # center at 2000
ungroup() %>%
select(Trait, Outcome, study, continent, country, scale, baseAge, baseAgeSD, baseYear, baseYearSD, cogYear, cogYearSD, predInt, predIntSD) %>%
distinct() %>%
group_by(study, Outcome) %>%
summarize_at(vars(baseAge:predIntSD), mean) %>%
kable(., "html", digits = 2) %>%
kable_classic(full_width = F, html_font = "Times New Roman") | study | Outcome | baseAge | baseAgeSD | baseYear | baseYearSD | cogYear | cogYearSD | predInt | predIntSD |
|---|---|---|---|---|---|---|---|---|---|
| BASE | crystallized | 78.23 | 6.66 | 1990 | 0.00 | 1997.00 | 3.51 | 8.56 | 3.51 |
| EAS | crystallized | 79.47 | 5.36 | 2011 | 3.18 | 2014.00 | 3.09 | 2.15 | 2.41 |
| GSOEP | crystallized | 49.83 | 15.83 | 2005 | 0.00 | 2012.00 | 0.00 | 7.00 | 0.00 |
| HILDA | crystallized | 44.51 | 16.92 | 2005 | 0.00 | 2012.00 | 0.00 | 7.00 | 0.00 |
| HRS | crystallized | 71.71 | 6.97 | 2006 | 0.00 | 2010.00 | 0.00 | 4.00 | 0.00 |
| LASA | crystallized | 61.46 | 15.71 | 1992 | 1.19 | 1995.00 | 9.09 | 8.39 | 9.45 |
| MAP | crystallized | 79.45 | 7.32 | 0 | 0.00 | 6.33 | 4.53 | 6.75 | 4.53 |
| MARS | crystallized | 73.60 | 6.37 | 0 | 0.00 | 6.00 | 3.86 | 6.46 | 3.86 |
| OCTO-Twin | crystallized | 82.99 | 2.66 | 1991 | 0.00 | 1997.00 | 2.82 | 6.21 | 2.82 |
| ROS | crystallized | 75.87 | 7.38 | 0 | 0.03 | 8.00 | 6.42 | 9.53 | 6.42 |
| SATSA | crystallized | 54.77 | 9.84 | 1984 | 0.00 | 1999.00 | 0.00 | 15.00 | 0.00 |
5.2.2 Harmonize Data
ipd3_meta_data <- ipd3_meta_data %>%
group_by(study, Trait, Outcome) %>%
mutate_at(vars(p_value, o_value, SRhealth),
~((. - min(., na.rm = T))/(max(., na.rm = T) - min(., na.rm = T))*10)) %>%
mutate(gender = factor(gender, levels = c(0,1), labels = c("Male", "Female")),
education = education - 12,
age = age - mean(age, na.rm = T)) %>%
ungroup()
ipd3_meta_data## # A tibble: 119,597 × 13
## study Trait p_value p_year SID Outcome o_value o_year education gender SRhealth yearBrth age
## <chr> <chr> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <fct> <dbl> <dbl> <dbl>
## 1 BASE E 7 1990 10004 crystallized 2.17 2000 5 Male 7.5 1918 -6.23
## 2 BASE E 7 1990 10010 crystallized 3.08 1995 -1 Male 2.5 1911 0.766
## 3 BASE E 3 1990 10033 crystallized 0.909 1997 -4 Female 10 1910 1.77
## 4 BASE E 4.5 1990 10034 crystallized 6.54 1995 -2 Male 10 1914 -2.23
## 5 BASE E 3 1990 10111 crystallized 6.54 1995 1 Female 7.5 1901 10.8
## 6 BASE E 5 1990 10115 crystallized 10 1995 -1 Male 2.5 1918 -6.23
## 7 BASE E 4.5 1990 10116 crystallized 3.85 1995 -2 Male 7.5 1917 -5.23
## 8 BASE E 4.5 1990 10145 crystallized 6.36 1997 -2 Female 10 1897 14.8
## 9 BASE E 5.5 1990 10175 crystallized 2.31 1995 -4 Female 5 1911 0.766
## 10 BASE E 5.5 1990 10188 crystallized 5.77 2004 -2 Female 2.5 1903 8.77
## # ℹ 119,587 more rows5.2.3 Save Data Files
save_fun <- function(d, trait, outcome, study){
save(d, file = sprintf("%s/data/two_stage/%s_%s_%s.RData", local_path, trait, outcome, study))
}
ipd3_meta_data %>%
group_by(study, Trait, Outcome) %>%
nest() %>%
ungroup() %>%
mutate(data = pmap(list(data, Trait, Outcome, study), save_fun))## # A tibble: 40 × 4
## study Trait Outcome data
## <chr> <chr> <chr> <list>
## 1 BASE E crystallized <NULL>
## 2 BASE N crystallized <NULL>
## 3 BASE O crystallized <NULL>
## 4 EAS A crystallized <NULL>
## 5 EAS C crystallized <NULL>
## 6 EAS E crystallized <NULL>
## 7 EAS N crystallized <NULL>
## 8 EAS O crystallized <NULL>
## 9 GSOEP A crystallized <NULL>
## 10 GSOEP C crystallized <NULL>
## # ℹ 30 more rows5.3 Step 2: Run Models for Each Study
5.3.1 Functions
5.3.1.1 Model Function
ipd3_study_mod_fun <- function(trait, outcome, type, mod, study, cov){
## load the data
load(sprintf("%s/data/two_stage/%s_%s_%s.RData", local_path, trait, outcome, study))
## compiled Bayesian model to speed up processing and avoid crashing
if(type == "Bayesian") load(sprintf("%s/results/3_ipd_meta/bayes_sample_mod.RData", local_path))
## model formula
if (cov == "all") cv <- c("age", "gender", "education")
if (!cov %in% c("all", "none")) cv <- cov
rhs <- "p_value"
rhs <- if(cov != "none") c(rhs, cv) else rhs
if(mod != "none"){rhs <- c(rhs, paste("p_value", mod, sep = "*"))}
rhs <- paste(rhs, collapse = " + ")
f <- paste("o_value ~ ", rhs, collapse = "")
## run the models & save
m <- if(type == "Frequentist"){do.call("lm", list(formula = f, data = quote(d)))} else {update(m, formula = f, newdata = d, warmup = 1000, iter = 2000, cores = 1)}
save(m, file = sprintf("%s/results/3_ipd_meta/%s/studyModels/%s_%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov, study))
## extract model terms and confidence intervals & save
rx <- tidy(m, conf.int = T) %>%
select(term, estimate, conf.low, conf.high)
save(rx, file = sprintf("%s/results/3_ipd_meta/%s/studySummary/%s_%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov, study))
## calculate effect sizes for random effects meta analysis
es <- ipd3_es_fun(m, type, mod)
save(es, file = sprintf("%s/results/3_ipd_meta/%s/studyEffects/%s_%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov, study))
## calculate simple effects
if(mod != "none"){
pred.rx <- ipd3_study_simpeff_fun(m, mod, type)
save(pred.rx, file = sprintf("%s/results/3_ipd_meta/%s/studyPredicted/%s_%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov, study))
}
## clean up the local function environment
rm(list = c("d", "f", "m", "fx", "es", "rhs"))
gc()
}5.3.1.2 Effect Size Function
ipd3_es_fun <- function(m, type, mod){
## extract model features needed for meta-analysis
ts <- insight::clean_parameters(m)$Cleaned_Parameter
ts <- if(mod == "none") "p_value" else ts[grepl("p_value.", ts)]
ts <- str_replace(ts, "[.]", ":")
## standardize the model
# ms <- standardize(m)
## get standardized model coefficients and standard errors
## for bayesian models this is the sd of the posterior estimates
es <- if(type == "Frequentist"){
summary(m)$coef[c("(Intercept)", ts), c("Estimate", "Std. Error")] %>% as.data.frame() %>% setNames(c("Estimate", "SEI"))
} else {
fixef(m)[c("Intercept", ts), c("Estimate", "Est.Error")] %>% as.data.frame() %>% setNames(c("Estimate", "SEI"))
}
## format to standardized format
es <- es %>%
rownames_to_column("term") %>%
# t() %>%
# as.data.frame() %>%
mutate(ni = if(type == "Frequentist") {nrow(m$model)} else {nrow(m$data)})
return(es)
}5.3.1.3 Study-Specific Simple Effects Function
ipd3_study_simpeff_fun <- function(m, moder, type){
d <- if(type == "Bayesian") m$data else m$model
d <- d %>% select(-o_value, -p_value)
cols <- colnames(d)
md_cl <- class(d[,moder])
if(any(sapply(d, class) == "numeric")){
msd <- d %>%
select_if(is.numeric) %>%
pivot_longer(everything()
, names_to = "item"
, values_to = "value") %>%
group_by(item) %>%
summarize_at(vars(value), lst(mean, sd)) %>%
ungroup()
}
if(any(sapply(d, class) == "factor")){
fct_lev <- d %>%
select_if(is.factor) %>%
summarize_all(~list(unique(.)))
}
d <- d %>% select(-one_of(moder))
md_levs <- if(md_cl == "numeric"){
if(moder %in% c("age", "baseAge", "baseYear")) {
c(-10, 0, 10)
} else if (moder %in% c("predInt", "education")) {
c(-5, 0, 5)
} else {
with(msd, c(mean[item == moder] - sd[item == moder], mean[item == moder], mean[item == moder] + sd[item == moder]))
}
} else {
unique(fct_lev[,moder][[1]])
}
mod_frame <- crossing(
p_value = seq(0,10,.5)
, modvalue = md_levs
) %>% setNames(c("p_value", moder))
if(ncol(d) > 0){
if(any(sapply(d, class) == "numeric")){
mod_frame <- tibble(mod_frame, d %>% select_if(is.numeric) %>% summarize_all(mean))
}
if(any(sapply(d, class) == "factor")){
mod_frame <- tibble(mod_frame, d %>% select_if(is.factor) %>% summarize_all(~levels(.)[1]))
}
}
pred.fx <- if(type == "Bayesian"){
bind_cols(
mod_frame,
fitted(m, newdata = mod_frame) %>% data.frame
) %>%
select(one_of(colnames(m$data)), pred = Estimate, lower = Q2.5, upper = Q97.5)
} else {
bind_cols(
mod_frame,
predict(m, newdata = mod_frame, interval = "confidence") %>% data.frame
) %>%
select(one_of(colnames(m$model)), pred = fit, lower = lwr, upper = upr)
}
rm(list = c("m", "mod_frame", "d", "md_levs"))
gc()
return(pred.fx)
}5.3.2 Run Models
# done <- tibble(file = list.files(sprintf("%s/results/3_ipd_meta/Bayesian/studyModels", local_path))) %>%
# separate(file, c("Outcome", "Trait", "Moderator", "Covariate", "study"), sep = "_") %>%
# mutate(study = str_remove_all(study, ".RData"),
# done = "done")
plan(multisession(workers = 12L))
nested_ipd_reg <- tibble(files = list.files(sprintf("%s/data/two_stage", local_path))) %>%
separate(files, c("Trait", "Outcome", "study"), sep = "_") %>%
filter(!is.na(study)) %>%
mutate(study = str_remove(study, ".RData")) %>%
left_join(
crossing(
study = studies
, Trait = traits$short_name
, Outcome = outcomes$short_name
, type = c("Frequentist", "Bayesian")
, Moderator = c("age", "gender", "education")
, Covariate = c("none", "all")
) %>%
full_join(
crossing(
study = studies
, Trait = traits$short_name
, Outcome = outcomes$short_name
, type = c("Frequentist", "Bayesian")
, Moderator = "none"
, Covariate = c("none", "age", "gender", "education", "all")
)
)
) %>%
filter(!is.na(Covariate)) %>%
# full_join(done) %>% filter(is.na(done)) %>%
# filter(type == "Frequentist") %>%
filter(Trait == "N" & study == "HILDA") %>%
mutate(run =
# pmap(list(Trait, Outcome, type, Moderator, study, Covariate)
future_pmap(list(Trait, Outcome, type, Moderator, study, Covariate)
, possibly(ipd3_study_mod_fun, NA_real_)
, .progress = T
, .options = furrr_options(
globals = c("ipd3_es_fun"
, "ipd3_study_simpeff_fun"
, "read_path"
, "local_path"
, "res_path"
, "codebook"
, "covars"
, "moders"
, "outcomes"
, "studies"
, "stdyModers"
, "traits"
, "data_path")
, packages = c("lme4"
, "broom"
, "psych"
, "knitr"
, "broom.mixed"
, "brms"
#, "tidybayes"
#, "bootpredictlme4"
, "rstan"
, "estimatr"
#, "merTools"
, "plyr"
, "tidyverse"))
))
closeAllConnections()5.4 Step 3: Meta-Analyze Results
5.4.1 Functions
5.4.1.0.1 Meta-Analysis Function
ipd3_meta_fun <- function(es, type, trait, outcome, mod, cov){
print(paste(type, trait, outcome, mod, cov))
## bayesian sample models for stability and speed
if(type == "Bayesian") {
mr <- if(mod %in% stdyModers$short_name) "metareg" else "meta"
load(sprintf("%s/results/3_ipd_meta/bayes_sample_%s.RData", local_path, mr))
}
## adding meta-regression values to effect sizes
if(mod %in% stdyModers$short_name) {
load(sprintf("%s/data/two_stage/meta_data/%s_%s.RData", local_path, trait, outcome))
es <- d %>% select(study, one_of(mod)) %>%
setNames(c("study", "metamod")) %>%
full_join(es) %>%
mutate_if(is.character, factor)
es0 <- es %>% filter(grepl("Intercept", term)) %>% select(-term)
}
es <- es %>% filter(!grepl("Intercept", term)) %>% select(-term)
## base bayesian model
# brm(formula = bf(Estimate | se(SEI) ~ 1 + (1 | study))
# , save_pars = "all"
# , sample_prior = T
# , prior = prior(cauchy(0,1), class = sd)
# , iter = 4000)
## base bayesian meta-regression model
# update(mt
# , formula. = bf(~ . + metamod)
# , newdata = es
# , sample_prior = T)
# run the meta-analytic model
f <- if (mod %in% c("none", moders$short_name)) "Estimate ~ 1" else "Estimate ~ 1 + metamod"
mt <- if(type == "Frequentist"){
rma(formula(f)
, sei = SEI
, ni = ni
, slab = study
, data = es)
} else {
if(!mod %in% stdyModers$short_name) {
update(m
, newdata = es
, iter = 2000
, warmup = 1000
, cores = 1)
} else {
update(m
, formula. = bf(~ . + metamod)
, newdata = es
, iter = 2000
, warmup = 1000
, cores = 1)
}
}
if(mod %in% stdyModers$short_name) {
mt0 <- if(type == "Frequentist") {
rma(formula(f)
, sei = SEI
, ni = ni
, slab = study
, data = es0)
} else {
update(m
, formula. = bf(~ . + metamod)
, newdata = es0
, iter = 2000
, warmup = 1000
, cores = 1)
}
}
save(mt, file = sprintf("%s/results/3_ipd_meta/%s/metaModels/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
# pull out and format the fixed effects (i.e. overall effects)
fx <- ipd3_meta_fx_fun(mt, type, mod)
save(fx, file = sprintf("%s/results/3_ipd_meta/%s/metaSummary/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
# meta-analysis predictions
if(mod %in% stdyModers$short_name) {
pred.fx0 <- ipd3_meta_simpeff_fun(mt0, mod, type) %>% setNames(c(mod, "b0_pred", "b0_lower", "b0_higher"))
pred.fx1 <- ipd3_meta_simpeff_fun(mt, mod, type) %>% setNames(c(mod, "b1_pred", "b1_lower", "b1_higher"))
pred.fx <- crossing(
p_value = seq(0, 10, .25)
, pred.fx0 %>%
full_join(pred.fx1)
) %>%
mutate(pred = b0_pred + b1_pred*p_value) %>%
select(p_value, all_of(mod), pred)
save(pred.fx, file = sprintf("%s/results/3_ipd_meta/%s/metaPredicted/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
}
# pull out and format the cross-study heterogeneity estimates
het <- ipd3_meta_rx_fun(mt, type)
save(het, file = sprintf("%s/results/3_ipd_meta/%s/metaHetero/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
rm(list = c("mt", "es", "type", "trait", "outcome", "mod", "fx", "rx"))
return(NULL)
}5.4.1.0.2 Meta-Analysis Fixed Effect Function
ipd3_meta_fx_fun <- function(mt, type, mod){
trgt <- if(mod %in% c("none", stdyModers$short_name)) "p_value" else paste0("p_value:", mod)
if (type == "Frequentist"){
coef(summary(mt)) %>%
rownames_to_column("term") %>%
select(term, estimate, SE = se, conf.low = ci.lb, conf.high = ci.ub) %>%
mutate(study = "Meta-Analytic",
term = mapvalues(term, "intrcpt", trgt))
} else {
fixef(mt) %>% data.frame() %>%
rownames_to_column("term") %>%
select(term, estimate = Estimate, SE = Est.Error, conf.low = Q2.5, conf.high = Q97.5) %>%
mutate(study = "Meta-Analytic",
term = mapvalues(term, "Intercept", trgt))
}
}5.4.1.0.3 Meta-Analysis Heterogeneity Function
ipd3_meta_rx_fun <- function(mt, type){
if (type == "Frequentist"){
## for frequentist, we'll grab estimates of:
## - tau^2: estimated between-study heterogeneity
## - I^2: total hetero (tau^2) / total hetero + total var (tau^2 + sigma^2)
## - H^2: total var (tau^2 + sigma^2) / total sampling var (sigma^2)
## - QE: Chi^2 dist Cochran's Q statistic (hetero > 0)
## - QEp: associated p-value for QE for k df
mt[c("tau2", "se.tau2", "I2", "H2", "QE", "QEp")] %>%
ldply() %>%
pivot_wider(names_from = ".id", values_from = "V1")
} else {
## for Bayesian, we'll grab estimates of:
## note, for these, we must estimate some directly
## but will use Bayes Factor to estimate probability of tau^2 > 0
## this should converge with other estimates but is more appropriate for Bayes
## - tau^2: average estimated between-study hetero across Bayes samples
## - I^2: total hetero (tau^2) / total hetero + total var (tau^2 + sigma^2)
## - H^2: total var (tau^2 + sigma^2) / total sampling var (sigma^2)
## - BF: posterior prob / prior prob
tibble(tau2 = summary(mt)$random$study[,"Estimate"]^2
, se.tau2 = summary(mt)$random$study[,"Est.Error"]^2
, I2 = tau2 / (tau2 + var(resid(mt)[,"Estimate"]))
, H2 = (tau2 + var(resid(mt)[,"Estimate"])) / var(resid(mt)[,"Estimate"])
, BF = 1/hypothesis(mt, "study__Intercept^2 = 0", class = "sd")$hypothesis$Evid.Ratio
)
}
}5.4.1.0.4 Meta-Analysis Simple Effects Function
## seemingly only seems to make sense for meta regressions
## the effect sizes uses would be the difference in the association
## associated with a 1 SD change in the moderator
## or with a dummy code?\
ipd3_meta_simpeff_fun <- function(m, moder, type){
d <- if(type == "Bayesian") m$data %>% select(metamod) else m$X %>% data.frame() %>% select(-intrcpt)
md_cl <- class(d$metamod)
if(any(sapply(d, class) == "numeric")){
msd <- d %>%
select_if(is.numeric) %>%
pivot_longer(everything()
, names_to = "item"
, values_to = "value") %>%
group_by(item) %>%
summarize_at(vars(value), lst(mean, sd)) %>%
ungroup()
}
if(any(sapply(d, class) == "factor")){
fct_lev <- d %>%
select_if(is.factor) %>%
summarize_all(~list(unique(.)))
}
# d <- d %>% select(-one_of(moder))
md_levs <- if(md_cl == "numeric"){
if(moder %in% c("age", "baseAge", "baseYear")) {
c(-10, 0, 10)
} else if (moder %in% c("predInt", "education")) {
c(-5, 0, 5)
} else {
with(msd, c(mean[item == moder] - sd[item == moder], mean[item == moder], mean[item == moder] + sd[item == moder]))
}
} else {
if(type == "Bayesian"){
unique(fct_lev[,moder][[1]])
} else {
unique(d) %>% as.matrix()
}
}
if(moder %in% c("country", "continent", "scale")) rownames(md_levs) <- 1:nrow(md_levs)
probs <- tribble(
~wrong , ~correct,
"metamodNEO.FFI" , "metamodNEO-FFI",
"metamodBFI.S" , "metamodBFI-S",
"metamodIPIP.NEO", "metamodIPIP NEO",
"metamodTDA.40" , "metamodTDA-40",
"metamodThe.Netherlands" , "metamodThe Netherlands"
)
mod_frame <- if(type == "Bayesian") {
expand.grid(
SEI = 0
, metamod = md_levs
, stringsAsFactors = F
)
} else {
md_levs
}
if(moder %in% c("scale", "country")) colnames(mod_frame) <- mapvalues(colnames(mod_frame), probs$wrong, probs$correct)
pred.fx <- if(type == "Bayesian"){
bind_cols(
mod_frame,
fitted(m
, newdata = mod_frame
, re_formula = NA
) %>% data.frame
) %>%
select(metamod, pred = Estimate, lower = Q2.5, upper = Q97.5)
} else {
bind_cols(
tibble(modvalue = mod_frame),
predict(m, mod_frame) %>% data.frame()
) %>%
select(modvalue, pred, lower = ci.lb, upper = ci.ub)
}
if(moder %in% c("country", "scale", "continent")){
mssng <- if(moder == "country") "United States" else if(moder == "continent") "North America" else "NEO-FFI"
nms <- str_remove(colnames(mod_frame), "metamod")
trck <- apply(mod_frame, 1, function(x){
if(sum(x) > 0) which(x == 1) else ncol(mod_frame) + 1
}); trck <- c(nms, mssng)[trck]
rownames(mod_frame) <- trck
pred.fx <- pred.fx %>% mutate(modvalue = trck)
}
rm(list = c("m", "mod_frame", "d", "md_levs"))
gc()
return(pred.fx)
}5.4.2 Run Meta-Analysis and Meta-Regression Models
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/3_ipd_meta/%s/%s/%s", local_path, type, folder, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
## load in effect size data
## first get file names
nested_ipd3_meta <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/3_ipd_meta/%s/studyEffects", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate", "study"), sep = "_", remove = F) %>%
filter(!is.na(study)) %>%
# filter(Trait == "N") %>%
## read in the files
mutate(study = str_remove(study, ".RData"),
data = map2(file, type, ~loadRData(.x, .y, "es", "studyEffects"))) %>%
select(-file) %>%
## unnest effect sizes
unnest(data)
nested_ipd3_meta## # A tibble: 1,712 × 10
## type Outcome Trait Moderator Covariate study term Estimate SEI ni
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <int>
## 1 Frequentist crystallized A age all EAS (Intercept) 5.67 0.344 788
## 2 Frequentist crystallized A age all EAS p_value:age -0.0103 0.00881 788
## 3 Frequentist crystallized A age all GSOEP (Intercept) 8.08 0.206 647
## 4 Frequentist crystallized A age all GSOEP p_value:age -0.00168 0.00195 647
## 5 Frequentist crystallized A age all HILDA (Intercept) 5.13 0.110 7755
## 6 Frequentist crystallized A age all HILDA p_value:age 0.000451 0.000817 7755
## 7 Frequentist crystallized A age all HRS (Intercept) 5.14 0.112 8432
## 8 Frequentist crystallized A age all HRS p_value:age -0.00142 0.00185 8432
## 9 Frequentist crystallized A age all ROS (Intercept) 7.02 0.284 1372
## 10 Frequentist crystallized A age all ROS p_value:age 0.0136 0.00651 1372
## # ℹ 1,702 more rows## group and nest
nested_ipd3_meta <- nested_ipd3_meta %>%
group_by(type, Trait, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup()
## add in meta-moderators, which requires taking original models and
## modifying the Moderator column
nested_ipd3_meta <- nested_ipd3_meta %>%
filter(Moderator == "none") %>%
select(-Moderator) %>%
full_join(crossing(
Trait = traits$short_name
, Moderator = stdyModers$short_name)) %>%
full_join(nested_ipd3_meta)mod <- "age"
mr <- if(mod %in% stdyModers$short_name) "metareg" else "meta"
## adding meta-regression values to effect sizes
es <- (nested_ipd3_meta %>% filter(Trait == "N" & Covariate == "all" & Moderator == "age" & type == "Bayesian"))$data[[1]]
# base bayesian model
m <- brm(formula = bf(Estimate | se(SEI) ~ 1 + (1 | study))
, data = es
, save_pars = save_pars(all = T)
, sample_prior = T
, prior = prior(cauchy(0,1), class = sd)
, iter = 30
, warmup = 20)
save(m, file = sprintf("%s/results/3_ipd_meta/bayes_sample_meta.RData", local_path))
mod <- "scale"
load(sprintf("%s/data/two_stage/meta_data/N_episodicMem.RData", local_path))
es <- d %>% select(study, one_of(mod)) %>%
setNames(c("study", "metamod")) %>%
full_join(es)
m <- brm(formula = bf(Estimate | se(SEI) ~ 1 + metamod + (1 | study))
, data = es
, save_pars = save_pars(all = T)
, sample_prior = T
, prior = prior(cauchy(0,1), class = sd)
, iter = 30
, warmup = 20)
save(m, file = sprintf("%s/results/3_ipd_meta/bayes_sample_metareg.RData", local_path))
rm(list = c("d", "es", "mod", "m", "mr"))plan(multisession(workers = 12L))
nested_ipd3_meta <- nested_ipd3_meta %>%
# full_join(done) %>% filter(is.na(done)) %>%
filter(type == "Bayesian") %>%
# filter(Moderator %in% c("scale", "country", "continent")) %>%
mutate(metamod =
# pmap(list(data, type, Trait, Outcome, Moderator, Covariate)
future_pmap(list(data, type, Trait, Outcome, Moderator, Covariate)
, possibly(ipd3_meta_fun, NA_real_)
, .progress = T
, .options = furrr_options(
globals = c("ipd3_meta_fx_fun"
, "ipd3_meta_rx_fun"
, "ipd3_meta_simpeff_fun"
, "read_path"
, "local_path"
, "res_path"
, "codebook"
, "covars"
, "moders"
, "outcomes"
, "studies"
, "stdyModers"
, "traits"
, "data_path")
, packages = c("lme4"
, "broom"
, "psych"
, "knitr"
, "broom.mixed"
, "brms"
#, "tidybayes"
#, "bootpredictlme4"
, "rstan"
, "estimatr"
, "metafor"
, "plyr"
, "tidyverse"))
))
closeAllConnections()
nested_ipd3_metacontr_fun <- function(m, std){
cntrm <- c("p_value = 0", "p_value + p_value:genderFemale = 0")
(multcomp::glht(m, cntrm) %>% # multcomp hypothesis function
confint(., calpha = multcomp::univariate_calpha()))$confint %>%
data.frame() %>%
mutate(term = c("male", "female")) %>%
rename(estimate = Estimate, conf.low = lwr, conf.high = upr)
}
res <- nested_ipd3_meta %>%
filter(type == "Frequentist" & Covariate == "all") %>%
mutate(res = map2(model, study, contr_fun))
res %>%
select(-model) %>%
unnest(res) %>%
mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns")) %>%
filter(sig == "sig")5.4.3 Compile Results
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/3_ipd_meta/%s/%s/%s", local_path, type, folder, fileName)
# print(path)
load(path)
get(ls()[grepl(obj, ls())])
}
## load in effect size data
## first get file names
nested_ipd3_meta <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/3_ipd_meta/%s/studyEffects", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate", "study"), sep = "_", remove = F) %>%
filter(!is.na(study)) %>%
## read in the files
mutate(study = str_remove(study, ".RData"),
studyEff = map2(file, type, ~loadRData(.x, .y, "x", "studySummary")),
n = map2(file, type, ~loadRData(.x, .y, "es", "studyEffects")),
n = map_dbl(n, ~(.)$ni[1])) %>%
select(-file) %>%
unnest(studyEff) %>%
group_by(type, Trait, Outcome, Moderator, Covariate) %>%
nest(studyEff = study:n) %>%
ungroup()
## now we add in the study-level moderators (i.e. meta-regression)
nested_ipd3_meta <- nested_ipd3_meta %>%
filter(Moderator == "none") %>%
select(-Moderator) %>%
full_join(crossing(
Trait = traits$short_name
, Moderator = stdyModers$short_name)) %>%
full_join(nested_ipd3_meta) %>%
mutate(file = sprintf("%s_%s_%s_%s.RData", Outcome, Trait, Moderator, Covariate),
metaEff = map2(file, type, ~loadRData(.x, .y, "fx", "metaSummary")),
metaHet = map2(file, type, possibly(~loadRData(.x, .y, "het", "metaHetero"), NA_real_))) %>%
select(-file)This results in a nested data frame, with columns:
studyEff= standardized study-specific effects from stage 1 regressions
metaEff= standardized meta-analytic effect from stage 2 meta-analysis
metaHet= Measures of cross-study heterogeneity
## # A tibble: 410 × 8
## type Outcome Trait Covariate studyEff Moderator metaEff metaHet
## <chr> <chr> <chr> <chr> <list> <chr> <list> <list>
## 1 Frequentist crystallized A age <tibble [18 × 6]> baseAge <df [2 × 6]> <tibble [1 × 6]>
## 2 Frequentist crystallized A age <tibble [18 × 6]> baseYear <df [2 × 6]> <tibble [1 × 6]>
## 3 Frequentist crystallized A age <tibble [18 × 6]> continent <df [3 × 6]> <tibble [1 × 6]>
## 4 Frequentist crystallized A age <tibble [18 × 6]> country <df [4 × 6]> <tibble [1 × 6]>
## 5 Frequentist crystallized A age <tibble [18 × 6]> predInt <df [2 × 6]> <tibble [1 × 6]>
## 6 Frequentist crystallized A age <tibble [18 × 6]> scale <df [2 × 6]> <dbl [1]>
## 7 Frequentist crystallized A all <tibble [30 × 6]> baseAge <df [2 × 6]> <tibble [1 × 6]>
## 8 Frequentist crystallized A all <tibble [30 × 6]> baseYear <df [2 × 6]> <tibble [1 × 6]>
## 9 Frequentist crystallized A all <tibble [30 × 6]> continent <df [3 × 6]> <tibble [1 × 6]>
## 10 Frequentist crystallized A all <tibble [30 × 6]> country <df [4 × 6]> <tibble [1 × 6]>
## # ℹ 400 more rows5.4.3.0.1 Tables
First, we’ll combine study and meta-analytic results.
ipd3_meta_res <- nested_ipd3_meta %>%
mutate(comEff = map2(studyEff, metaEff, ~(.y) %>% full_join(.x))) %>%
select(type, Outcome, Trait, Moderator, Covariate, comEff) %>%
unnest(comEff) %>%
filter((Moderator == "none" & term == "p_value") |
(!Moderator %in% stdyModers$short_name & (grepl("p_value:", term))) |
(Moderator %in% stdyModers$short_name & grepl("metamod", term))) %>%
mutate(term = str_replace_all(term, "metamod", paste0("p_value:", Moderator)),
study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin"))) %>%
select(-SE)
ipd3_meta_res## # A tibble: 1,635 × 11
## type Outcome Trait Moderator Covariate term estimate conf.low conf.high study n
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 Frequentist crystallized A baseAge age p_value:baseAge 0.00400 -0.00174 0.00973 Meta-Ana… NA
## 2 Frequentist crystallized A baseYear age p_value:baseYear 0.0102 0.00399 0.0164 Meta-Ana… NA
## 3 Frequentist crystallized A continent age p_value:continentAustralia -0.0489 -0.211 0.113 Meta-Ana… NA
## 4 Frequentist crystallized A continent age p_value:continentEurope -0.177 -0.314 -0.0391 Meta-Ana… NA
## 5 Frequentist crystallized A country age p_value:countryAustralia -0.0309 -0.0775 0.0158 Meta-Ana… NA
## 6 Frequentist crystallized A country age p_value:countryGermany -0.0873 -0.133 -0.0418 Meta-Ana… NA
## 7 Frequentist crystallized A country age p_value:countrySweden -0.268 -0.377 -0.158 Meta-Ana… NA
## 8 Frequentist crystallized A predInt age p_value:predInt -0.0198 -0.0331 -0.00647 Meta-Ana… NA
## 9 Frequentist crystallized A baseAge all p_value:baseAge 0.00164 -0.00104 0.00431 Meta-Ana… NA
## 10 Frequentist crystallized A baseYear all p_value:baseYear 0.00507 0.000720 0.00943 Meta-Ana… NA
## # ℹ 1,625 more rows5.4.3.0.1.1 Study-Specific
Next, we’ll make a table of the results, separately for each moderator. To do this efficiently, we’ll make a function that creates those tables across all combinations. Before calling that function, we’ll do some reformatting and reshaping to get the data ready.
ipd3_res_tab <- ipd3_meta_res %>%
mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns")) %>% # significance marker
mutate_at(vars(estimate:conf.high),
~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f",.))) %>%
mutate(est = sprintf("%s<br>[%s, %s]", estimate, conf.low, conf.high),
est = ifelse(sig == "sig", sprintf("<strong>%s</strong>", est), est),
study = factor(study, c(studies_long, "Meta-Analytic")),
Outcome = factor(Outcome, outcomes$short_name, outcomes$long_name)) %>%
select(type:Covariate, term, study, est) %>%
pivot_wider(names_from = "Trait", values_from = "est") %>%
select(type:study, E, A, C, N, O) %>%
arrange(type, Outcome, Moderator, Covariate, study)
ipd3_res_tab## # A tibble: 414 × 11
## type Outcome Moderator Covariate term study E A C N O
## <chr> <fct> <chr> <chr> <chr> <fct> <chr> <chr> <chr> <chr> <chr>
## 1 Bayesian Crystallized Ability age all p_value:age BASE 0.002<br>[-0.04, 0.04] <NA> <NA> 0.01… -0.0…
## 2 Bayesian Crystallized Ability age all p_value:age EAS -0.005<br>[-0.02, 0.01] -0.0… <str… <str… -0.0…
## 3 Bayesian Crystallized Ability age all p_value:age GSOEP -0.000<br>[-0.004, 0.003] -0.0… 0.00… 0.00… -0.0…
## 4 Bayesian Crystallized Ability age all p_value:age HILDA <strong>0.001<br>[0.000, … 0.00… <str… <str… <str…
## 5 Bayesian Crystallized Ability age all p_value:age HRS -0.002<br>[-0.005, 0.001] -0.0… <str… 0.00… <str…
## 6 Bayesian Crystallized Ability age all p_value:age LASA <NA> <NA> <NA> -0.0… <NA>
## 7 Bayesian Crystallized Ability age all p_value:age MAP -0.002<br>[-0.009, 0.005] <NA> 0.00… -0.0… <NA>
## 8 Bayesian Crystallized Ability age all p_value:age MARS <NA> <NA> <NA> -0.0… <NA>
## 9 Bayesian Crystallized Ability age all p_value:age OCTO-Twin -0.01<br>[-0.06, 0.03] <NA> <NA> 0.03… <NA>
## 10 Bayesian Crystallized Ability age all p_value:age ROS 0.004<br>[-0.007, 0.01] <str… 0.00… <str… <str…
## # ℹ 404 more rowsipd3_stdmeta_table_fun <- function(d, type, moder, cov){
cv <- mapvalues(cov, covars$short_name, covars$long_name, warn_missing = F)
if(!grepl("djust", cv)) cv <- paste(cv, "Adjusted")
md <- mapvalues(moder, c(moders$short_name, stdyModers$short_name),
c(moders$long_name, stdyModers$long_name), warn_missing = F)
rs <- d %>% group_by(Outcome) %>% tally() %>%
mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
cs <- rep(1,6); names(cs) <- c(" ", traits$short_name)
cap <- if(moder == "none") {
sprintf("<strong>Table X.</strong><br><em>Method 3 Pooled Regression Using Random Effects: Study and Meta-Analytic Estimates of %s Personality-Cognitive Domain Relationships</em>", cv)
} else {
sprintf("<strong>Table X.</strong><br><em>Method 3 Pooled Regression Using Random Effects: Study and Meta-Analytic %s Moderation of %s Personality-Cognitive Domain Relationships</em>", md, cv)
}
tab <- d %>%
select(-Outcome) %>%
kable(., "html"
, escape = F
, booktabs = T
, col.names = c("Study", rep("$\\beta$ [CI]", 5))
, align = c("r", rep("c", 5))
, caption = cap) %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
add_header_above(cs)
for (i in 1:nrow(rs)){
tab <- tab %>%
kableExtra::group_rows(rs$Outcome[i], rs$start[i], rs$end[i])
}
save_kable(tab, file = sprintf("%s/results/3_ipd_meta/%s/tables/study specific/%s_%s.html",
local_path, type, moder, cov))
return(tab)
}
ipd3_std_tab <- ipd3_res_tab %>%
filter(!Moderator %in% stdyModers$short_name) %>%
select(-term) %>%
group_by(type, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Moderator, Covariate), ipd3_stdmeta_table_fun))5.4.3.0.1.2 Meta-Analytic
ipd3_tab_fun <- function(d, type, moder){
md <- mapvalues(moder, c(moders$long_name, stdyModers$long_name), c(moders$short_name, stdyModers$short_name), warn_missing = F)
rs <- d %>% group_by(Outcome) %>% tally() %>%
mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
cs <- if(length(unique(d$term)) == 1) rep(1,6) else c(2, rep(1,5))
names(cs) <- c(" ", traits$long_name)
cln <- if(length(unique(d$term)) == 1) c("CovariaCtes", rep("<em>b</em> [CI]", 5)) else c("Covariates", "Term", rep("<em>b</em> [CI]", 5))
# cln <- if(length(unique(d$term)) == 1) c("Covariates", rep("\\textit{b} [CI]", 5)) else c("Covariates", "Term", rep("\\textit{b} [CI]", 5))
al <- if(length(unique(d$term)) == 1) c("r", rep("c", 5)) else c("r", "r", rep("c", 5))
if(length(unique(d$term)) == 1) d <- d %>% select(-term)
cap <- if(md == "none") "<strong>Table X.</strong><br><em>Method 3 Pooled Regression Using Random Effects: Meta-Analytic Effects of Personality-Crystallized Domain Associations</em>" else sprintf("<strong>Table X.</strong><br><em>Method 3 Pooled Regression Using Random Effects: Meta-Analytic %s Moderation of Personality-Crystallized Domain Associations</em>", md)
tab <- d %>%
arrange(Outcome) %>%
select(-Outcome) %>%
kable(., "html"
# kable(., "latex"
, booktabs = T
, escape = F
, col.names = cln
, align = al
, caption = cap
) %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
add_header_above(cs)
for (i in 1:nrow(rs)) {
tab <- tab %>% kableExtra::group_rows(rs$Outcome[i], rs$start[i], rs$end[i])
}
save_kable(tab, file = sprintf("%s/results/3_ipd_meta/%s/tables/overall/%s.html"
, local_path, type, md))
return(tab)
}
ipd3_meta_res_tab <- ipd3_res_tab %>%
filter(study == "Meta-Analytic") %>%
select(-study) %>%
mutate(term = str_remove_all(term, "p_value:"),
term = mapvalues(term, str_remove_all(stdyModers$short_term, " "), stdyModers$long_term),
term = mapvalues(term, str_replace_all(stdyModers$short_term, "-", "M"), stdyModers$long_term),
term = mapvalues(term, moders$short_term, moders$long_term),
term = mapvalues(term, moders$short_name, moders$long_term),
Covariate = mapvalues(Covariate, covars$short_name, covars$long_name)) %>%
arrange(Outcome, term, Covariate)
ipd3_meta_res_tab %>%
group_by(type, Moderator) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Moderator), ipd3_tab_fun))## # A tibble: 20 × 4
## type Moderator data tab
## <chr> <chr> <list> <list>
## 1 Bayesian age <tibble [2 × 8]> <kablExtr [1]>
## 2 Frequentist age <tibble [2 × 8]> <kablExtr [1]>
## 3 Bayesian continent <tibble [10 × 8]> <kablExtr [1]>
## 4 Frequentist continent <tibble [10 × 8]> <kablExtr [1]>
## 5 Bayesian country <tibble [20 × 8]> <kablExtr [1]>
## 6 Frequentist country <tibble [20 × 8]> <kablExtr [1]>
## 7 Bayesian education <tibble [2 × 8]> <kablExtr [1]>
## 8 Frequentist education <tibble [2 × 8]> <kablExtr [1]>
## 9 Bayesian gender <tibble [2 × 8]> <kablExtr [1]>
## 10 Frequentist gender <tibble [2 × 8]> <kablExtr [1]>
## 11 Bayesian none <tibble [5 × 8]> <kablExtr [1]>
## 12 Frequentist none <tibble [5 × 8]> <kablExtr [1]>
## 13 Bayesian predInt <tibble [5 × 8]> <kablExtr [1]>
## 14 Frequentist predInt <tibble [5 × 8]> <kablExtr [1]>
## 15 Bayesian scale <tibble [30 × 8]> <kablExtr [1]>
## 16 Frequentist scale <tibble [30 × 8]> <kablExtr [1]>
## 17 Bayesian baseAge <tibble [5 × 8]> <kablExtr [1]>
## 18 Frequentist baseAge <tibble [5 × 8]> <kablExtr [1]>
## 19 Bayesian baseYear <tibble [5 × 8]> <kablExtr [1]>
## 20 Frequentist baseYear <tibble [5 × 8]> <kablExtr [1]>ipd3_tab_fun <- function(d, type, cov){
# long outcome name
covar <- mapvalues(cov, covars$long_name, covars$short_name, warn_missing = F)
# getting row numbers for later grouping
rs <- d %>% group_by(Moderator) %>% tally() %>%
mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
# number and name of columns for span columns
cs <- rep(1,6)
names(cs) <- c(" ", traits$long_name)
# cln <- if(length(unique(d$term2)) == 1) c("Covariate", rep("\\textit{b} [CI]", 5)) else c(" ", "Term", rep("\\textit{b} [CI]", 5))
cln <- c("Term", rep("<em>b</em> [CI]", 5))
al <- c("r", rep("c", 5))
# caption
cap <- sprintf("Method 3 Pooled Regression Using Random Effects: Meta-Analytic Estimates of %s Personality-Crystallized Domain Associations", cov)
# kable the table
tab <- d %>%
select(-Moderator) %>%
kable(., "html"
# kable(., "latex"
, booktabs = T
, escape = F
, col.names = cln
, align = al
, caption = cap
) %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
# kable_styling(full_width = F, font_size = 7) %>%
add_header_above(cs)
# for loop to add grouped sections
for (i in 1:nrow(rs)){
tab <- tab %>%
kableExtra::group_rows(rs$Moderator[i], rs$start[i], rs$end[i])
}
# save the resulting html table
save_kable(tab, file = sprintf("%s/results/3_ipd_meta/%s/tables/key terms/%s.html"
, local_path, type, covar))
return(tab) # return the html table
}
ipd3_fx_tab2 <- ipd3_res_tab %>%
filter(study == "Meta-Analytic") %>%
select(-study) %>%
mutate(term = str_remove_all(term, "p_value:"),
term = mapvalues(term, str_remove_all(stdyModers$short_term, " "), stdyModers$long_term),
term = mapvalues(term, str_replace_all(stdyModers$short_term, "-", "M"), stdyModers$long_term),
term = mapvalues(term, moders$short_term, moders$long_term),
term = mapvalues(term, moders$short_name, moders$long_term),
Moderator = factor(Moderator, levels = c(moders$short_name, stdyModers$short_name), labels = c(moders$long_name, stdyModers$long_name)),
Covariate = mapvalues(Covariate, covars$short_name, covars$long_name)) %>%
arrange(Outcome, Moderator, term, Covariate) %>%
filter(Covariate %in% c("Unadjusted", "Fully Adjusted")) %>%
group_by(Outcome, type, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Covariate), ipd3_tab_fun))
# save(ipd1a_fx_tab, ipd1a_fx_tab2, ipd1a_res, file = sprintf("%s/manuscript/results/ipd1b_fx_tab.RData", res_path))5.4.3.0.1.3 All Model Terms
ipd3_mod_tab <- nested_ipd3_meta %>%
select(-metaEff, -metaHet) %>%
unnest(studyEff) %>%
# keep key terms
# mark significance and prettify trait, outcome, and covariate names
mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns"),
Trait = factor(Trait, traits$short_name),
Outcome = factor(Outcome, outcomes$short_name, outcomes$long_name),
Moderator = factor(Moderator, c(moders$short_name, stdyModers$short_name), c(moders$long_name, stdyModers$long_name)),
Covariate = factor(Covariate, covars$short_name, str_wrap(covars$long_name, 15)),
term = str_replace_all(term, "metamod", paste0("p_value:", Moderator))) %>%
# format values as text, combine estimates and CI's, bold significance
mutate_at(vars(estimate, conf.low, conf.high),
~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
mutate(est = sprintf("%s<br>[%s, %s]", estimate, conf.low, conf.high),
est = ifelse(sig == "sig", sprintf("<strong>%s</strong>", est), est),
study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin"))) %>%
# mutate(est = sprintf("%s [%s, %s]", estimate, conf.low, conf.high),
# est = ifelse(sig == "sig", sprintf("\\textbf{%s}", est), est)) %>%
# final reshaping, remove extra columns, arrange values, and change to wide format
select(-estimate, -conf.low, -conf.high, -sig, -n) %>%
arrange(type, Outcome, Trait, Moderator, Covariate) %>%
pivot_wider(names_from = "Trait", values_from = "est")
ipd3_mod_tab_fun <- function(d, type, out, moder, cov){
md <- mapvalues(moder, c(moders$long_name, stdyModers$long_name), c(moders$short_name, stdyModers$short_name), warn_missing = F)
o <- mapvalues(out, outcomes$long_name, outcomes$short_name, warn_missing = F)
cv <- mapvalues(cov, covars$long_name, covars$short_name, warn_missing = F)
cs <- rep(1,6)
names(cs) <- c(" ", traits$long_name)
# cln <- if(length(unique(d$term2)) == 1) c("Covariate", rep("\\textit{b} [CI]", 5)) else c(" ", "Term", rep("\\textit{b} [CI]", 5))
cln <- c("Term", rep("<em>b</em> [CI]", 5))
al <- c("r", rep("c", 5))
# caption
cap <- if(md == "none") "3 Separate Models Followed Random Effects Meta-Analysis: All Model Estimates of Fixed Effect Personality-Crystallized Domain Associations" else sprintf("3 Separate Models Followed Random Effects Meta-Analysis: All Model Estimates of Fixed Effect %s Moderation of Personality-Crystallized Domain Associations", md)
d <- d %>% arrange(study, term)
rs <- d %>% group_by(study) %>% tally() %>%
mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
# kable the table
tab <- d %>%
select(-study) %>%
kable(., "html"
# kable(., "latex"
, booktabs = T
, escape = F
, col.names = cln
, align = al
, caption = cap
) %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
# kable_styling(full_width = F, font_size = 7) %>%
add_header_above(cs)
for(i in 1:nrow(rs)){
tab <- tab %>% kableExtra::group_rows(rs$study[i], rs$start[i], rs$end[i])
}
# save the resulting html table
save_kable(tab, file = sprintf("%s/results/3_ipd_meta/%s/tables/all terms/%s-%s-%s.html"
, local_path, type, o, md, cv))
return(tab) # return the html table
}
ipd3_mod_tab2 <- ipd3_mod_tab %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Outcome, Moderator, Covariate), ipd3_mod_tab_fun))ipd3_meta_all <- nested_ipd3_meta %>%
select(-studyEff, -metaHet) %>%
unnest(metaEff) %>%
mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns"),
Trait = factor(Trait, traits$short_name),
Outcome = factor(Outcome, outcomes$short_name, outcomes$long_name),
Moderator = factor(Moderator, c(moders$short_name, stdyModers$short_name), c(moders$long_name, stdyModers$long_name)),
Covariate = ifelse(Moderator != "None" & Covariate == "None", Moderator, Covariate),
Covariate = factor(Covariate, moders$short_name, str_wrap(moders$long_name, 15)),
term = str_replace_all(term, "metamod", paste0("p_value:", Moderator))) %>%
# format values as text, combine estimates and CI's, bold significance
mutate_at(vars(estimate, conf.low, conf.high),
~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
mutate(est = sprintf("%s<br>[%s, %s]", estimate, conf.low, conf.high),
est = ifelse(sig == "sig", sprintf("<strong>%s</strong>", est), est)) %>%
# mutate(est = sprintf("%s [%s, %s]", estimate, conf.low, conf.high),
# est = ifelse(sig == "sig", sprintf("\\textbf{%s}", est), est)) %>%
# final reshaping, remove extra columns, arrange values, and change to wide format
select(-estimate, -conf.low, -conf.high, -sig, -SE) %>%
arrange(type, Outcome, Trait, Moderator, Covariate) %>%
pivot_wider(names_from = "Trait", values_from = "est") %>%
group_by(type, Outcome, Covariate) %>%
nest() %>%
ungroup() 5.4.3.0.1.4 Heterogeneity
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/3_ipd_meta/%s/%s/%s", local_path, type, folder, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
## load in "fixed" effects
## first get file names
nested_ipd3_het <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/3_ipd_meta/%s/metaHetero", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate"), sep = "_", remove = F) %>%
## read in the files
mutate(Covariate = str_remove(Covariate, ".RData"),
het = map2(file, type, possibly(~loadRData(.x, .y, "het", "metaHetero"), NA_real_))) %>%
filter(!is.na(het)) %>%
select(-file) round_fun <- function(x){
ifelse(x < .001 & x > 0, "< .001"
, ifelse(x > -.001 & x < 0, "> -.001"
, ifelse(abs(x) < .01, sprintf("%.3f", x)
, sprintf("%.2f", x))))
}
ip3_hetero_tab_fun <- function(d, type, out, mod, cov){
moder <- mapvalues(mod, moders$short_name, moders$long_name, warn_missing = F)
d2 <- d %>%
mutate(Trait = factor(Trait, traits$short_name, traits$long_name)) %>%
mutate_at(vars(-Trait), round_fun) %>%
arrange(Trait)
cap <- if(mod == "none") "Heterogeneity Estimates of Personality-Crystallized Domain Associations" else sprintf("Heterogeneity Estimates for Overall %s Moderation of Personality-Crystallized Domain Associations", moder)
cap <- sprintf("<strong>Table SX</strong><br><em>%s</em>", cap)
tab <- d2 %>%
kable(.
, "html"
, align = c("r", rep("c", ncol(d2)-1))
, caption = cap
, escape = F
) %>%
kable_classic(full_width = F, html_font = "Times New Roman")
save_kable(tab, file = sprintf("%s/results/3_ipd_meta/%s/tables/heterogeneity/%s-%s-%s.html", local_path, type, out, mod, cov))
return(tab)
}
nested_ipd3_het_tab <- nested_ipd3_het %>%
filter(!Moderator %in% stdyModers$short_name) %>%
group_by(type) %>%
nest() %>%
ungroup() %>%
mutate(data = map(data, ~(.) %>%
unnest(het) %>%
group_by(Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup())) %>%
unnest(data) %>%
mutate(tab = pmap(list(data, type, Outcome, Moderator, Covariate), ip3_hetero_tab_fun))5.4.4 Figures
5.4.4.0.0.1 Overall Forest
ipd3_fx_plot_fun <- function(df, mod, type, cov){
m <- mapvalues(mod, c(moders$long_name, stdyModers$long_name), c(moders$short_name, stdyModers$short_name), warn_missing = F)
cv <- mapvalues(cov, covars$short_name, covars$long_name, warn_missing = F)
d <- round(max(abs(min(df$estimate)), abs(max(df$estimate))), 3)
lim <- c(0-d-(d/2.5), 0+d+(d/2.5))
brk <- if(d > .01) round(c(0-d-(d/5), 0, 0+d+(d/5)),2) else round(c(0-d-(d/5), 0, 0+d+(d/5)),3)
# lim_high <- lim[2]*4
lab <- str_replace(brk, "^0.", ".")
shapes <- c(15, 16, 17, 18)[1:length(unique(df$term))]
lt <- rep("solid", length(unique(df$term)))
titl <- if(mod == "none"){NULL} else {sprintf("%s Moderation of Personality-Outcome Associations", mod)}
titl <- if(!cov %in% c("none", "all")) paste(cv, "Adjusted", titl, collapse = " ") else paste(cv, titl, collapse = " ")
leg <- if(length(unique(df$term)) > 1){"bottom"} else {"none"}
p <- df %>%
mutate(conf.low = ifelse(conf.low < lim[1], lim[1], conf.low),
conf.high = ifelse(conf.high > lim[2], lim[2], conf.high)) %>%
ggplot(aes(x = Outcome, y = estimate)) +
scale_y_continuous(limits = lim, breaks = brk, labels = lab) +
scale_size_manual(values = c(1.2, .85)) +
scale_shape_manual(values = shapes) +
scale_color_manual(values = c("blue", "black")) +
scale_linetype_manual(values = lt) +
geom_hline(aes(yintercept = 0), size = .25, color = "gray50") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high, linetype = term)
, width = 0
, position = position_dodge(width = .9)) +
geom_point(aes(color = sig, size = sig, shape = term)
, position = position_dodge(width = .9)) +
labs(x = NULL
, y = "Estimate (POMP)"
, title = titl
, subtitle = "Method 3: Pooled Regression Using Random Effects"
) +
guides(size = "none", color = "none") +
facet_grid(~Trait, scales = "free_y", space = "free") +
coord_flip() +
theme_classic() +
theme(legend.position = leg,
plot.title = element_text(face = "bold", size = rel(1.2), hjust = .5),
plot.subtitle = element_text(size = rel(1.1), hjust = .5),
panel.background = element_rect(color = "black", fill = "white"),
strip.background = element_blank(),
strip.text = element_text(face = "bold", color = "black", size = rel(1.4)),
axis.text = element_text(color = "black"),
axis.text.y = element_text(size = rel(1)))
ht <- length(unique(df$Outcome)); ht2 <- length(unique(df$term))
local_path <- length(unique(df$Trait))
ggsave(file = sprintf("%s/results/3_ipd_meta/%s/figures/overall forest/%s_%s_fixed.png", local_path, type, m, cov), width = local_path*2, height = 1.25*ht + .75*ht2)
rm(p)
gc()
return(T)
}
nested_ipd3_meta %>%
select(-studyEff, -metaHet) %>%
unnest(metaEff) %>%
filter((Moderator == "none" & term == "p_value")|
(Moderator != "none" & grepl("^p_value:", term))) %>%
mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns"),
sig = factor(sig, levels = c("sig","ns")),
Trait = factor(Trait, levels = traits$short_name),
Outcome = factor(Outcome, levels = outcomes$short_name, labels = str_wrap(outcomes$long_name, 15)),
Moderator = factor(Moderator, levels = c(moders$short_name, stdyModers$short_name), labels = c(moders$long_name, stdyModers$long_name)),
Outcome = forcats::fct_rev(Outcome)) %>%
filter(type == "Frequentist" & Moderator == "None") %>%
group_by(type, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(pmap(list(data, Moderator, type, Covariate), ipd3_fx_plot_fun))5.4.4.0.0.2 Study-Specific Forest
ipd3_rx_plot_fun <- function(df, outcome, mod, type, cov, trait){
print(paste(outcome, mod))
trt <- mapvalues(trait, traits$short_name, traits$long_name)
m <- mapvalues(mod, moders$short_name, moders$long_name, warn_missing = F)
d <- round(max(abs(min(df$estimate)), abs(max(df$estimate))), 3)
# stds <- unique(df$study)
lim <- c(0-d-(d/2.5), 0+d+(d/2.5))
brk <- round(c(0-d-(d/5), 0, 0+d+(d/5)),2)
lab <- str_remove(round(c(0-d-(d/5), 0, 0+d+(d/5)),2), "^0")
shapes <- c(15, 16, 17, 18)[1:length(unique(df$term))]
lt <- rep("solid", length(unique(df$term)))
titl <- if(mod == "none"){trt} else {sprintf("%s x %s", trt, m)}
leg <- if(length(unique(df$term)) > 1){"bottom"} else {"none"}
df <- df %>% mutate(study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin")))
df <- df %>% full_join(tibble(study = " ", estimate = NA, n = NA))
df <- df %>% arrange(estimate)
stds <- df$study[!df$study %in% c("Meta-Analytic", " ")]
df <- df %>%
mutate(study = factor(study, rev(c(" ", stds, "Meta-Analytic")))
# , conf.low = ifelse(conf.low < lim[1], lim[1], conf.low)
# , conf.high = ifelse(conf.high > lim[2], lim[2], conf.high)
, lb = ifelse(conf.low < lim[1], "lower"
, ifelse(conf.high > lim[2], "upper", "neither"))
, conf.low2 = ifelse(conf.low < lim[1], lim[1], conf.low)
, conf.high2 = ifelse(conf.high > lim[2], lim[2], conf.high)
# , study = factor(study, levels = str_remove_all(c("Overall", studies_long), "-"), labels = c("Overall", studies_long))
# Trait = factor(Trait, levels = traits$short_name, labels = traits$long_name),
, type = ifelse(study == "Meta-Analytic", "fixed", "random"))
p1 <- df %>%
ggplot(aes(x = study, y = estimate)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high)
, position = "dodge"
, width = .2) +
geom_point(aes(shape = term, size = term)) +
geom_segment(data = df %>% filter(lb == "lower")
, aes(y = conf.high2, yend = conf.low2, xend = study)
, arrow = arrow(type = "closed", length = unit(0.1, "cm"))) +
geom_segment(data = df %>% filter(lb == "upper")
, aes(y = conf.low2, yend = conf.high2, xend = study)
, arrow = arrow(type = "closed", length = unit(0.1, "cm"))) +
geom_hline(aes(yintercept = 0), linetype = "dashed", size = .5) +
geom_vline(aes(xintercept = 1.5)) +
geom_vline(aes(xintercept = length(stds) + 1.5)) +
annotate("rect", xmin = length(stds) + 1.6, xmax = Inf, ymin = -Inf, ymax = Inf, fill = "white") +
scale_y_continuous(limits = lim, breaks = brk, labels = lab) +
scale_size_manual(values = c(3,2)) +
scale_shape_manual(values = c(15, 16)) +
labs(x = NULL
, y = "Estimate"
# , title = " "
) +
coord_flip() +
theme_classic() +
theme(legend.position = "none"
, axis.text = element_text(face = "bold")
, axis.title = element_text(face = "bold")
, plot.title = element_text(face = "bold", hjust = .5)
, axis.ticks.y = element_blank()
, axis.line.y = element_blank()
, axis.line.x.top = element_line(size = 1))
d2 <- df %>%
mutate_at(vars(estimate, conf.low, conf.high)
, ~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
mutate_at(vars(estimate, conf.low, conf.high), ~str_replace_all(., "^0.", ".")) %>%
mutate_at(vars(estimate, conf.low, conf.high), ~str_replace_all(., "^-0.", "-.")) %>%
mutate(est = ifelse(study != " ", sprintf("%s [%s, %s] ", estimate, conf.low, conf.high), "")
, n = as.character(n)
) %>%
select(study, n, est) %>%
pivot_longer(cols = c(n, est), names_to = "est", values_to = "value")
p2 <- d2 %>%
ggplot(aes(x = study, y = est)) +
geom_text(aes(label = value), hjust = .5, size = 3.5) +
annotate("text", label = "b [CI]", x = length(stds) + 1.75, y = "est", hjust = .5, vjust = 0) +
annotate("text", label = "N", x = length(stds) + 1.75, y = "n", hjust = .5, vjust = 0) +
geom_vline(aes(xintercept = 1.5)) +
geom_vline(aes(xintercept = length(stds) + 1.5)) +
coord_flip() +
theme_void() +
theme(plot.title = element_text(face = "bold", hjust = 0)
, axis.text = element_blank()
, axis.ticks = element_blank()
, axis.title = element_blank())
my_theme <- function(...) {
theme_classic() +
theme(plot.title = element_text(face = "italic"))
}
title_theme <- calc_element("plot.title", my_theme())
ttl <- ggdraw() +
draw_label(
titl,
fontfamily = title_theme$family,
fontface = title_theme$face,
size = title_theme$size-2
)
p3 <- cowplot::plot_grid(p1, p2
, rel_widths = c(.5, .5)#c(.4, .6)
, align = "h"
)
# p <- cowplot::plot_grid(ttl, subttl, p3, rel_heights = c(.05, .05, .9), nrow = 3)
p <- cowplot::plot_grid(ttl, p3, rel_heights = c(.05, .95), nrow = 2)
save(p
, file = sprintf("%s/results/3_ipd_meta/%s/figures/study specific forest/rdata/%s_%s_%s_%s.RData", local_path, type, outcome, trait, mod, cov))
gc()
return(p)
}
## fixed effects
nested_ipd3_reg_fp <- ipd3_meta_res %>%
filter(Moderator %in% moders$short_name) %>%
## filter key terms
filter((Moderator == "none" & term == "p_value")|
(Moderator != "none" & grepl("^p_value:", term) & !grepl("study", term))) %>%
## significance
mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns")
, study = mapvalues(study, studies_long, studies_sp, warn_missing = F)) %>%
## grouping for plotting
group_by(Outcome, Moderator, type, Covariate, Trait) %>%
nest() %>%
ungroup() %>%
mutate(p = pmap(list(data, Outcome, Moderator, type, Covariate, Trait), ipd3_rx_plot_fun))
ipd3_rx_plot_comb_fun <- function(outcome, cov, mod, type, d){
o <- mapvalues(outcome, outcomes$short_name, outcomes$long_name, warn_missing = F)
cv <- mapvalues(cov, covars$short_name, covars$long_name, warn_missing = F)
titl <- paste0(o, ",")
titl <- if(!cov %in% c("none", "all")) paste(titl, cv, "Adjusted", collapse = ", ") else paste(titl, cv, collapse = ", ")
p1 <- plot_grid(
d$p[[1]]
, d$p[[2]]
, d$p[[3]]
, d$p[[4]]
, d$p[[5]]
, nrow = 3
, ncol = 2
, axis = "tblr"
, align = "hv"
)
my_theme <- function(...) {
theme_classic() +
theme(plot.title = element_text(face = "bold"))
}
title_theme <- calc_element("plot.title", my_theme())
ttl <- ggdraw() +
draw_label(
titl,
fontfamily = title_theme$family,
fontface = title_theme$face,
size = title_theme$size
)
my_theme <- function(...) {
theme_classic() +
theme(plot.subtitle = element_text(hjust = 0))
}
subtitle_theme <- calc_element("subplot.title", my_theme())
subttl <- ggdraw() +
draw_label(
"Method 3: Two-Stage Individual Participant Meta-Analysis",
fontfamily = subtitle_theme$family,
fontface = subtitle_theme$face,
size = subtitle_theme$size
)
p <- cowplot::plot_grid(ttl, subttl, p1, rel_heights = c(.03, .03, .94), nrow = 3)
ggsave(p
, file = sprintf("%s/results/3_ipd_meta/%s/figures/study specific forest/%s_%s_%s.png", local_path, type, outcome, mod, cov)
, width = 10, height = 10)
ggsave(p
, file = sprintf("%s/results/3_ipd_meta/%s/figures/study specific forest/%s_%s_%s.pdf", local_path, type, outcome, mod, cov)
, width = 10, height = 10)
return(T)
}
nested_ipd3_reg_fp %>%
mutate(Trait = factor(Trait, traits$short_name)) %>%
arrange(Trait) %>%
select(-data) %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(p = pmap(list(Outcome, Covariate, Moderator, type, data), ipd3_rx_plot_comb_fun))5.4.4.1 Overall Simple Effects Plots
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/3_ipd_meta/%s/%s/%s", local_path, type, folder, fileName)
# print(path)
load(path)
get(ls()[grepl(obj, ls())])
}
## load in effect size data
## first get file names
nested_ipd3_meta <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/3_ipd_meta/%s/studyPredicted", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate", "study"), sep = "_", remove = F) %>%
filter(!is.na(study)) %>%
## read in the files
mutate(study = str_remove(study, ".RData"),
pred.rx = map2(file, type, ~loadRData(.x, .y, "pred.rx", "studyPredicted"))
, n = map2(file, type, ~loadRData(.x, .y, "es", "studyEffects"))
, n = map_dbl(n, ~(.)$ni[1])) %>%
select(-file) %>%
unnest(pred.rx) %>%
group_by(type, Trait, Outcome, Moderator, Covariate) %>%
nest(studyPred = study:n) %>%
ungroup()ipd3_pred_fx_fun <- function(d, mod, type, outcome, trait, cov){
d <- d %>% unclass %>% data.frame
d$mod_value <- d[,mod]
d <- d %>% select(-all_of(mod)) %>% as_tibble
if(class(d$mod_value) %in% c("factor", "character")){
d <- d %>% mutate(mod_fac = factor(mod_value))
} else{
d2 <- d %>%
select(study, mod_value) %>%
distinct() %>%
arrange(study, mod_value)
if(mod == "age") d2 <- d2 %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("-10 yrs", "M", "+10 yrs")))
else if(mod == "baseYear") d2 <- d2 %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("1990", "200)0", "2010")))
else if(mod == "baseAge") d2 <- d2 %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("50", "60", "70")))
else if(mod == "predInt") d2 <- d2 %>% mutate(mod_fac = factor(mod_value, levels = c(-5, 0, 5), labels = c("-5 yrs", "5 yrs", "+5 yrs")))
else if(mod == "education") d2 <- d2 %>% mutate(mod_fac = factor(mod_value, levels = c(-5, 0, 5), labels = c("-5 yrs", "12 years", "+5 yrs")))
else d2 <- d2 %>% mutate(mod_fac = factor(mod_value, levels = unique(mod_value), labels = c("-1 SD", "M", "+1 SD")))
d <- d %>% full_join(d2) %>% ungroup()
}
pred.fx <- d %>%
group_by(p_value, mod_fac) %>%
summarize_at(vars(pred, lower, upper), ~weighted.mean(., n)) %>%
ungroup()
save(pred.fx, file = sprintf("%s/results/3_ipd_meta/%s/metaPredicted/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
}
nested_ipd3_meta %>%
# filter(Moderator == "baseAge") %>%
mutate(pred = pmap(list(studyPred, Moderator, type, Outcome, Trait, Covariate), ipd3_pred_fx_fun))ipd3_se_plot_fun <- function(d, outcome, mod, cov, type){
# print(paste(int, mod, cov, random, imp))
o <- mapvalues(outcome, outcomes$short_name, outcomes$long_name, warn_missing = F)
cv <- mapvalues(cov, covars$short_name, covars$long_name, warn_missing = F)
m <- mapvalues(mod, moders$short_name, moders$long_name, warn_missing = F)
titl <- if(mod == "none"){sprintf("%s", o)} else {sprintf("%s: Personality x %s Simple Effects", o, m)}
d <- d %>%
mutate(study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin")),
study = factor(study, levels = stdcolors$studies))
std <- unique(d$study)
# cols <- (stdcolors %>% filter(studies %in% std))$colors
d <- d %>% unclass %>% data.frame
d$mod_value <- d[,mod]
d <- d %>% select(-all_of(mod)) %>% as_tibble
d <- if(class(d$mod_value) %in% c("factor", "character")){
d %>% mutate(mod_fac = factor(mod_value))
} else{
d %>%
full_join(
d %>% select(Trait,study, mod_value) %>% distinct() %>% arrange(study, mod_value) %>%
group_by(Trait, study) %>% mutate(mod_fac = factor(c("-1 SD", "M", "+1 SD"), levels = c("-1 SD", "M", "+1 SD"))) %>% ungroup()
)
}
lt <- c("dotted", "solid", "dashed")[1:length(unique(d$mod_fac))]
ht <- length(unique(d$mod_fac))
p <- d %>%
mutate(Trait = factor(Trait, levels = traits$short_name, labels = traits$long_name)) %>%
ggplot(aes(x = p_value, y = pred, group = interaction(Trait, mod_fac), linetype = mod_fac)) +
# geom_line(aes(color = study
# , group = interaction(study, mod_fac)
# , linetype = study)
# , size = 1) +
# geom_ribbon(aes(fill = study
# , group = interaction(study, mod_fac)
# , ymin = lower
# , ymax = upper)
# , alpha = .25) +
stat_smooth(aes(weight = n, linetype = mod_fac, fill = mod_fac)
, method = "lm"
, formula = y~x
, size = 1
, color = "black"
) +
facet_wrap(~Trait, ncol = 2) +
scale_y_continuous(limits = c(4,10)
, breaks = c(4, 6, 8, 10)
, labels = c(4, 6, 8, 10)) +
# scale_color_manual(values = cols) +
scale_linetype_manual(values = lt) +
labs(x = "Personality Score (POMP)"
, y = "Cognition Score (POMP)"
# , color = "Study"
, fill = m
, linetype = m
, title = titl
, subtitle = "Method 3: Two-Stage Individual Participant Meta-Analysis"
) +
theme_classic() +
theme(legend.position = "bottom"
, plot.title = element_text(face = "bold", size = rel(1.2), hjust = .5)
, plot.subtitle = element_text(size = rel(1.1), hjust = .5)
, strip.background = element_rect(fill = "black")
, strip.text = element_text(face = "bold", color = "white")
, axis.text = element_text(color = "black"))
ggsave(file = sprintf("%s/results/3_ipd_meta/%s/figures/overall simple effects/%s_%s_%s.png", local_path, type, outcome, mod, cov), width = 6, height = 6)
}
nested_ipd3_meta %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(data = map(data, ~(.) %>% unnest(studyPred))
, p = pmap(list(data, Outcome, Moderator, Covariate, type), ipd3_se_plot_fun))5.4.4.2 Study-Specific Simple Effects Plots
ipd3_std_se_plot_fun <- function(d, outcome, trait, mod, cov, type){
# print(paste(int, mod, cov, random, imp))
o <- mapvalues(outcome, outcomes$short_name, outcomes$long_name, warn_missing = F)
trt <- mapvalues(trait, traits$short_name, traits$long_name, warn_missing = F)
cv <- mapvalues(cov, covars$short_name, covars$long_name, warn_missing = F)
m <- mapvalues(mod, moders$short_name, moders$long_name, warn_missing = F)
titl <- if(mod == "none"){sprintf("%s: %s", o, trt)} else {sprintf("%s: %s x %s Simple Effects", o, trt, m)}
d <- d %>% mutate(study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin")),
study = factor(study, levels = stdcolors$studies))
std <- unique(d$study)
cols <- (stdcolors %>% filter(studies %in% std))$colors
lt <- (stdcolors %>% filter(studies %in% std))$lt
d <- d %>% unclass %>% data.frame
d$mod_value <- d[,mod]
d <- d %>% select(-all_of(mod)) %>% as_tibble
d <- if(class(d$mod_value) %in% c("factor", "character")){
d %>% mutate(mod_fac = factor(mod_value))
} else{
d %>%
full_join(
d %>% select(study, mod_value) %>% distinct() %>% arrange(study, mod_value) %>%
group_by(study) %>% mutate(mod_fac = factor(c("-1 SD", "M", "+1 SD"), levels = c("-1 SD", "M", "+1 SD"))) %>% ungroup()
)
}
ht <- length(unique(d$mod_fac))
p <- d %>%
ggplot(aes(x = p_value, y = pred)) +
geom_line(aes(color = study
, group = interaction(study, mod_fac)
, linetype = study)
, size = 1) +
# geom_ribbon(aes(fill = study
# , group = interaction(study, mod_fac)
# , ymin = lower
# , ymax = upper)
# , alpha = .25) +
stat_smooth(aes(weight = n)
, method = "lm"
, formula = y~x
, size = 1.2
, color = "darkslateblue"
) +
scale_y_continuous(limits = c(4,10)
, breaks = c(4, 6, 8, 10)
, labels = c(4, 6, 8, 10)) +
scale_color_manual(values = cols) +
scale_linetype_manual(values = lt) +
labs(x = "Personality Score (POMP)"
, y = "Cognition Score (POMP)"
, color = "Study"
, fill = "Study"
, linetype = "Study"
, title = titl
, subtitle = "Method 3: Two-Stage Individual Participant Meta-Analysis"
) +
facet_wrap(~mod_fac) +
theme_classic() +
theme(legend.position = "bottom"
, plot.title = element_text(face = "bold", hjust = .5, size = rel(.95))
, plot.subtitle = element_text(size = rel(1.1), hjust = .5)
, strip.background = element_rect(fill = "black")
, strip.text = element_text(face = "bold", color = "white")
, axis.text = element_text(color = "black"))
ggsave(file = sprintf("%s/results/3_ipd_meta/%s/figures/study specific simple effects/%s_%s_%s_%s.png", local_path, type, outcome, trt, mod, cov), width = 3*ht, height = 5)
}
nested_ipd3_meta %>%
mutate(p = pmap(list(studyPred, Outcome, Trait, Moderator, Covariate, type), ipd3_std_se_plot_fun))load("~/Documents/projects/data synthesis/crystallized/results/3_ipd_meta/Frequentist/metaModels/crystallized_C_country_all.RData")
coef(mt)## intrcpt metamodAustralia metamodGermany metamodSweden
## 0.06433784 -0.06106685 -0.06322996 -0.18061653cntrm <- rbind(
c(1,0,0,0)
, c(1,1,0,0)
, c(1,0,1,0)
, c(1,0,0,1)
); rownames(cntrm) <- c("United States", "Australia", "Germany", "Sweden")
(multcomp::glht(mt, cntrm) %>% # multcomp hypothesis function
confint(., calpha = multcomp::univariate_calpha()))$confint %>%
data.frame() %>%
data.frame() %>%
rownames_to_column("cntr") %>%
mutate(term = rownames(cntrm)) %>%
select(-cntr) %>%
mutate(est = ifelse(abs(Estimate) < .01, sprintf("b = %.3f, 95%% CI [%.3f, %.3f]", Estimate, lwr, upr), sprintf("b = %.2f, 95%% CI [%.2f, %.2f]", Estimate, lwr, upr)))## Estimate lwr upr term est
## 1 0.064337841 0.04167694 0.08699874 United States b = 0.06, 95% CI [0.04, 0.09]
## 2 0.003270996 -0.02297616 0.02951815 Australia b = 0.003, 95% CI [-0.023, 0.030]
## 3 0.001107879 -0.05502738 0.05724314 Germany b = 0.001, 95% CI [-0.055, 0.057]
## 4 -0.116278686 -0.20058736 -0.03197002 Sweden b = -0.12, 95% CI [-0.20, -0.03]load("~/Documents/projects/data synthesis/crystallized/results/3_ipd_meta/Frequentist/metaModels/crystallized_C_scale_all.RData")
coef(mt)## intrcpt metamodBFI-S metamodEysenck metamodIPIP NEO metamodMIDI metamodTDA-40
## 0.08889284 -0.08778496 -0.20517153 -0.12764963 -0.02316003 -0.08562185cntrm <- rbind(
c(1,0,0,0,0,0) # NEO-FFI
, c(1,1,0,0,0,0) # BFI-S
, c(1,0,1,0,0,0) # Eysenck
, c(1,0,0,1,0,0) # IPIP NEO
, c(1,0,0,0,1,0) # MIDI
, c(1,0,0,0,0,1) # TDA-40
); rownames(cntrm) <- c("NEO-FFI", "BFI-S", "Eysenck", "IPIP NEO", "MIDI", "TDA-40")
(multcomp::glht(mt, cntrm) %>% # multcomp hypothesis function
confint(., calpha = multcomp::univariate_calpha()))$confint %>%
data.frame() %>%
data.frame() %>%
rownames_to_column("cntr") %>%
mutate(term = rownames(cntrm)) %>%
select(-cntr) %>%
mutate(est = ifelse(abs(Estimate) < .01, sprintf("b = %.3f, 95%% CI [%.3f, %.3f]", Estimate, lwr, upr), sprintf("b = %.2f, 95%% CI [%.2f, %.2f]", Estimate, lwr, upr)))## Estimate lwr upr term est
## 1 0.088892841 0.03486531 0.14292037 NEO-FFI b = 0.09, 95% CI [0.03, 0.14]
## 2 0.001107879 -0.05488928 0.05710503 BFI-S b = 0.001, 95% CI [-0.055, 0.057]
## 3 -0.116278686 -0.20049547 -0.03206191 Eysenck b = -0.12, 95% CI [-0.20, -0.03]
## 4 -0.038756793 -0.13781345 0.06029986 IPIP NEO b = -0.04, 95% CI [-0.14, 0.06]
## 5 0.065732812 0.04024951 0.09121611 MIDI b = 0.07, 95% CI [0.04, 0.09]
## 6 0.003270996 -0.02267949 0.02922148 TDA-40 b = 0.003, 95% CI [-0.023, 0.029]load("~/Documents/projects/data synthesis/crystallized/results/3_ipd_meta/Frequentist/metaModels/crystallized_N_scale_all.RData")
coef(mt)## intrcpt metamodBFI-S metamodDPQ metamodEysenck metamodIPIP NEO metamodMIDI metamodTDA-40
## -0.11303919 0.07722236 0.06822839 0.07197000 0.06484409 0.09060320 -0.00330891cntrm <- rbind(
c(1,0,0,0,0,0,0) # NEO-FFI
, c(1,1,0,0,0,0,0) # BFI-S
, c(1,0,1,0,0,0,0) # DPQ
, c(1,0,0,1,0,0,0) # Eysenck
, c(1,0,0,0,1,0,0) # IPIP NEO
, c(1,0,0,0,0,1,0) # MIDI
, c(1,0,0,0,0,0,1) # TDA-40
); rownames(cntrm) <- c("NEO-FFI", "BFI-S", "DPQ", "Eysenck", "IPIP NEO", "MIDI", "TDA-40")
(multcomp::glht(mt, cntrm) %>% # multcomp hypothesis function
confint(., calpha = multcomp::univariate_calpha()))$confint %>%
data.frame() %>%
data.frame() %>%
rownames_to_column("cntr") %>%
mutate(term = rownames(cntrm),
sig = ifelse(sign(lwr) == sign(upr), "sig", "ns")) %>%
select(-cntr) %>%
mutate(est = ifelse(abs(Estimate) < .01, sprintf("b = %.3f, 95%% CI [%.3f, %.3f]", Estimate, lwr, upr), sprintf("b = %.2f, 95%% CI [%.2f, %.2f]", Estimate, lwr, upr))) %>%
arrange(sig, Estimate)## Estimate lwr upr term sig est
## 1 -0.04819510 -0.14247152 0.046081326 IPIP NEO ns b = -0.05, 95% CI [-0.14, 0.05]
## 2 -0.04106919 -0.10615703 0.024018652 Eysenck ns b = -0.04, 95% CI [-0.11, 0.02]
## 3 -0.03581683 -0.08170762 0.010073957 BFI-S ns b = -0.04, 95% CI [-0.08, 0.01]
## 4 -0.11634810 -0.14179851 -0.090897687 TDA-40 sig b = -0.12, 95% CI [-0.14, -0.09]
## 5 -0.11303919 -0.15448612 -0.071592256 NEO-FFI sig b = -0.11, 95% CI [-0.15, -0.07]
## 6 -0.04481079 -0.08127859 -0.008343002 DPQ sig b = -0.04, 95% CI [-0.08, -0.01]
## 7 -0.02243598 -0.04344011 -0.001431855 MIDI sig b = -0.02, 95% CI [-0.04, -0.00]5.4.5 Table Meta-Analytic Heterogeneity
Relative to previous tables, this is slightly complicated because Frequentist and Bayesian meta-analysis don’t use the same methods for estimating heterogeneity.
ipd3_het_tab <- nested_ipd3_meta %>%
select(-studyEff, -metaEff) %>%
mutate(Trait = factor(Trait, traits$short_name, traits$long_name),
Outcome = factor(Outcome, outcomes$short_name, outcomes$long_name)) %>%
arrange(type, Outcome, Trait, Moderator) %>%
group_by(type, Moderator) %>%
nest() %>%
ungroup() %>%
mutate(data = map(data, ~(.) %>% unnest(metaHet)))
ipd3_het_tab
## frequentist
ipd3_het_tab$data[[nrow(ipd3_het_tab)]]
## bayesian
ipd3_het_tab$data[[1]]ipd3_het_tab_fun <- function(d, type, moder){
rs <- d %>% group_by(Outcome) %>% tally() %>%
mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
if(type == "Frequentist"){
d %>%
mutate_at(vars(tau2, QEp), ~ifelse(. < .001, "< 0.001", sprintf("%.3f", .))) %>%
mutate_at(vars(I2:QE), ~sprintf("%.2f", .)) %>%
select(Trait, tau2, I2, H2, QE, QEp) %>%
kable(., "html"
, escape = F
, digits = 2
, col.names = c("Trait", "$\\tau^2$", "$I^2$", "$H^2$", "<em>Q</em>", "<em>p</em>")
, align = c("r", rep("c",5))
, caption = sprintf("Table X. Heterogeneity estimates for %s Models with %s Moderator", type, moder)) %>%
kable_classic(full_width = F, html_font = "Times New Roman")
} else{
d %>%
mutate_at(vars(tau2, BF), ~ifelse(. < .001, "< 0.001", sprintf("%.3f", .))) %>%
mutate_at(vars(I2, H2), ~sprintf("%.2f", .)) %>%
select(Trait, tau2, I2, H2, BF) %>%
kable(., "html"
, escape = F
, digits = 2
, col.names = c("Trait", "$\\tau^2$", "$I^2$", "$H^2$", "BF")
, align = c("r", rep("c",4))
, caption = sprintf("Table X. Heterogeneity estimates for %s Models with %s Moderator", type, moder)) %>%
kable_classic(full_width = F, html_font = "Times New Roman")
}
}5.5 Sample Results Section
We examined estimates of overall prospective associations between Big Five personality characteristics and crystallized abilities as well as participant and sample-level moderators of those associations using a two-stage individual participant meta-analysis of 11 studies. This allowed us to estimate both sample-specific and overall prospective associations between personality characteristics and cognitive ability as well as heterogeneity in those estimates. Fully adjusted Big Five personality trait-crystallized ability associations as well as participant- and sample-level moderators of these associations can be found in Table S5. Of the 20 key participant-level personality-crystallized ability associations and moderators of those associations, two (10%) were significant. Specifically, Openness was positively associated with later crystallized abilities, and gender moderated this association.
|
Extraversion
|
Agreeableness
|
Conscientiousness
|
Neuroticism
|
Openness to Experience
|
|
|---|---|---|---|---|---|
| Term | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| None | |||||
| Personality |
0.01 [-0.02, 0.05] |
0.02 [-0.001, 0.03] |
0.02 [-0.03, 0.07] |
-0.07 [-0.10, -0.04] |
0.17 [0.10, 0.25] |
| Age | |||||
| Age |
-0.000 [-0.003, 0.002] |
-0.001 [-0.003, 0.002] |
-0.000 [-0.005, 0.004] |
0.000 [-0.004, 0.005] |
0.001 [-0.006, 0.007] |
| Gender | |||||
| Gender (Male v Female) |
0.04 [-0.004, 0.09] |
0.02 [-0.02, 0.05] |
0.005 [-0.03, 0.04] |
-0.000 [-0.03, 0.03] |
-0.04 [-0.07, -0.01] |
| Education | |||||
| Education (Years) |
-0.004 [-0.02, 0.008] |
-0.008 [-0.02, 0.008] |
0.003 [-0.01, 0.02] |
-0.001 [-0.01, 0.008] |
-0.007 [-0.02, 0.005] |
| Continent | |||||
| Continent (North America v Australia) |
0.03 [-0.06, 0.12] |
-0.03 [-0.15, 0.08] |
-0.06 [-0.17, 0.06] |
-0.05 [-0.13, 0.04] |
0.06 [-0.21, 0.33] |
| Continent (North America v Europe) |
0.06 [-0.02, 0.14] |
-0.08 [-0.19, 0.03] |
-0.11 [-0.21, -0.006] |
0.02 [-0.04, 0.08] |
0.08 [-0.13, 0.28] |
| Country | |||||
| Country (United States v Australia) |
0.03 [-0.05, 0.12] |
-0.03 [-0.11, 0.06] |
-0.06 [-0.10, -0.03] |
-0.04 [-0.14, 0.06] |
0.06 [-0.22, 0.34] |
| Country (United States v Germany) |
0.09 [0.002, 0.18] |
-0.03 [-0.12, 0.07] |
-0.06 [-0.12, -0.003] |
0.02 [-0.08, 0.12] |
0.04 [-0.20, 0.28] |
| Country (United States v Sweden) |
0.03 [-0.06, 0.13] |
-0.15 [-0.28, -0.02] |
-0.18 [-0.27, -0.09] |
0.03 [-0.07, 0.13] |
0.15 [-0.16, 0.46] |
| Country (United States v The Netherlands) |
0.03 [-0.08, 0.14] |
||||
| Personality Scale | |||||
| Scale (NEO-FFI v Eysenck) |
0.02 [-0.07, 0.12] |
-0.21 [-0.31, -0.11] |
0.07 [-0.005, 0.15] |
0.07 [-0.37, 0.51] |
|
| Scale (NEO-FFI v BFI-S) |
0.10 [0.01, 0.18] |
-0.09 [-0.17, -0.010] |
0.08 [0.02, 0.14] |
-0.14 [-0.57, 0.28] |
|
| Scale (NEO-FFI v DPQ) |
0.07 [0.01, 0.12] |
||||
| Scale (NEO-FFI v IPIP NEO) |
0.08 [-0.04, 0.20] |
-0.13 [-0.24, -0.01] |
0.06 [-0.04, 0.17] |
0.04 [-0.39, 0.48] |
|
| Scale (NEO-FFI v MIDI) |
-0.04 [-0.11, 0.03] |
-0.02 [-0.08, 0.04] |
0.09 [0.04, 0.14] |
-0.13 [-0.55, 0.29] |
|
| Scale (NEO-FFI v TDA-40) |
0.03 [-0.05, 0.10] |
-0.09 [-0.15, -0.03] |
-0.003 [-0.05, 0.05] |
-0.01 [-0.43, 0.41] |
|
| Baseline Age | |||||
| Study Baseline Age |
-0.001 [-0.003, 0.002] |
0.002 [-0.001, 0.004] |
0.002 [-0.001, 0.006] |
-0.000 [-0.003, 0.002] |
0.001 [-0.005, 0.007] |
| Baseline Year | |||||
| Study Baseline Year |
0.001 [-0.004, 0.006] |
0.005 [0.001, 0.009] |
0.004 [-0.003, 0.01] |
-0.000 [-0.005, 0.004] |
-0.006 [-0.01, 0.003] |
| Prediction Interval | |||||
| Prediction Interval |
-0.003 [-0.01, 0.009] |
-0.007 [-0.02, 0.001] |
-0.008 [-0.02, 0.006] |
0.001 [-0.010, 0.01] |
0.007 [-0.01, 0.03] |
To better understand the moderation, we examined fully adjusted simple effects plots for tests of moderators of personality-cognitive ability associations in Figure S7. Although the figure demonstrates the simple effects for all of the Big Five, gender only significantly moderated openness to experience (b = -0.04, 95% CI [-0.07, -0.01]). In addition, the overall association between openness-crystallized ability association was positive for males and females. However, the interaction indicates that females had a significantly lower association than males.
Figure 5.1: Figure S7. Prospective sample-specific and overall associations between extraversion (in POMP units, 0-10) and crystallized abilities (in POMP units, 0-10) across genders (male, female). Different colors and line types indicate different samples. Thicker, black lines indicate the average association across samples, while thinner lines indicate sample-specific associations.
Finally, we examined the consistency of the prospective associations across samples for both personality-cognitive ability associations and the simple effects across samples for each moderator. Figure S8 presents the forest plots of the fully adjusted sample-specific and meta-analytic estimates of Big Five personality characteristic-crystallized ability associations across samples. These associations were very consistent for openness, with all samples except for ROS demonstrating a positive association between openness and crystallized abilities. Other estimates were less consistent, with three samples (HRS, MAP, and ROS) showing a positive association between conscientiousness and crystallized abilities, one sample (SATSA) showing a negative association, and two samples (HILDA, GSOEP) showing a null association.
Figure 5.2: Figure S8. Forest plot of fully adjusted prospective associations between Big Five personality characteristics and crystallized abilities across samples for using one-stage pooled integrative data analysis via effects coded contrasts. Overall point estimates (squares) represent the grand-mean estimates of the association across samples, while sample point estimates represent regression terms or linear combinations of regression terms. Error bars capture the 95% CI around the point estimate. Arrows indicate the confidence band was truncated to better visualize the estimates.
(Although gender did not moderate the association between extraversion and crystallized abilities in Method 3, for consistency with other methods, we include interpretation here.) Next, we examined the consistency of moderator associations across samples. Figure S7 displays the overall fully adjusted predicted crystallized abilities levels at different levels of extraversion across samples for those males and females for comparison to other methods. As is clear in the Figure, the overall trend suggests a null association for women and slight negative association for men, such that men who were higher in extraversion tended to score slightly worse on crystallized domain tasks. HRS (b = 0.04, 95% CI [-0.07, -0.02]) showed a negative association between extraversion and crystallized abilities for women, HILDA (b = 0.05, 95% CI [0.02, 0.08]) and GSOEP (b = 0.10, 95% CI [0.03, 0.17]) showed positive associations, and all other samples showed null associations. For men, HRS (b = -0.05, 95% CI [-0.08, -0.02]) and SATSA (b = -0.17, 95% CI [-0.30, -0.03]) demonstrated a negative association, only GSOEP demonstrated a positive association (b = 0.08, 95% CI [0.01, 0.14]), and all other samples demonstrated a null effect.
Among sample-level moderators, 15 of the 62 (24.19%) tested sample-level moderators were significant (fully adjusted). The most notable of these were country-level differences in the association between Conscientiousness and crystallized / knowledge domain cognitive ability and differences in the associations between Neuroticism and Conscientiousness’s associations with cognitive ability across different personality scales. Specifically, although there was a positive association between Conscientiousness and crystallized / knowledge domain cognitive ability in the United States (b = 0.06, 95% CI [0.04, 0.09]), that association was null (Germany, b = 0.001, 95% CI [-0.055, 0.057]; Australia, b = 0.003, 95% CI [-0.023, 0.030]) or negative (Sweden, b = -0.12, 95% CI [-0.20, -0.03]). Similarly, the association for Conscientiousness was positive when using the NEO-FFI (b = 0.09, 95% CI [0.03, 0.14]) and MIDI (b = 0.07, 95% CI [0.04, 0.09, negative when using the measures in SATSA and OCTO-TWIN (b = -0.12, 95% CI [-0.20, -0.03]). and null when using IPIP NEO (b = -0.04, 95% CI [-0.14, 0.06]), BFI-S (b = 0.001, 95% CI [-0.06, 0.06]), and TDA-40 (b = 0.003, 95% CI [-0.023, 0.029]). For Neuroticism, the association was negative for the NEO-FFI (b = -0.11, 95% CI [-0.15, -0.07]), MIDI (b = -0.02, 95% CI [-0.04, -0.00]), DPQ (b = -0.04, 95% CI [-0.08, -0.01]), and TDA-40 (b = -0.12, 95% CI [-0.14, -0.09]), and null for the IPIP NEO (b = -0.05, 95% CI [-0.14, 0.05]), BFI-S (b = -0.04, 95% CI [-0.08, 0.01]), and the measure used in SATSA and OCTO-TWIN (b = -0.04, 95% CI [-0.11, 0.02]). Simple effects plots for all sample-level moderators can be found in the online materials and web app.