Chapter 4 Method 2: Pooled One Stage Models with Study-Specific Effects
The structure of the data for models with and without study-specific effects is identical, so we can directly proceed to step 2 and run the models. However, for the purposes of making each chapter somewhat standalone, the code from Method 1, Step 1 setup will be duplicated below.
4.1 Step 1: Combine Data
loadRData <- function(fileName, type){
#loads an RData file, and returns it
path <- sprintf("%s/data/clean/%s_cleaned.RData", local_path, fileName)
load(path)
get(ls()[grepl(type, ls())])
}
ipd_reg_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))
ipd_reg_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> <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 rows4.1.1 Study-Level Moderators
First, we need to bring in a few study-level moderators from Table 4, which includes information about the continent and country of origin, and the type of personality scale. Then we’ll join those data with cleaned data from each study in order to calculate the average age at baseline, the average baseline year, and the average interval between personality and cogntiive measurements. Because these vary within studies across personality traits and outcomes, we calculate baseline age, year, and prediction interval for all combinations of these separately.
url <- "https://github.com/emoriebeck/data-synthesis-tutorial/raw/main/codebooks/crystallized_tables.xlsx"
destfile <- "tables.xlsx"
curl::curl_download(url, destfile)
ipd_reg_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")) %>%
# mutate(p_year = as.numeric(p_year)) %>%
right_join(ipd_reg_data) %>%
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 = mean(p_year) - 2000) %>% # center at 2000
ungroup()4.1.2 Harmonize Data
Next, we need to harmonize the data across studies. As we preregistered, continuous variables (personality, cognition, and self-rated health) will be calclulated as percentages of maximum possible separately for each study, Trait and outcome. Unlike standardization procedures, that have a mean of zero and unit variance and can be misleading when data are skewed, POMP does not rescale sample variance based on the observed data, which overly relies on deviations from the mean. Instead, POMP relies on the ratio between the difference between a score and the minimum and the maximum and minimum, or
POMP = \(\frac{observed-minimum}{maximum-minimum}\)*10.
In addition, gender will be dummy coded with male as the reference group, education will be centered at 12 years of education, and age will be grand mean-centered in each study (also for each trait and outcome combination).
ipd_reg_data <- ipd_reg_data %>%
group_by(study, Trait, Outcome) %>%
mutate_at(vars(p_value, o_value),
~((. - 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, # center at 12 years of education
age = age - mean(age, na.rm = T)) %>% # center
ungroup() %>%
select(-yearBrth)4.1.3 Save Data Files
Now, we’ll save these data files into separate data files for each trait, outcome combination. This makes it easier to track each data set that will be used in subsequent analyses.
save_fun <- function(d, trait, outcome){
save(d, file = sprintf("%s/data/one_stage/%s_%s.RData", local_path, trait, outcome))
}
ipd_reg_data %>%
group_by(Trait, Outcome) %>%
nest() %>%
ungroup() %>%
mutate(data = pmap(list(data, Trait, Outcome), save_fun))## # A tibble: 5 × 3
## Trait Outcome data
## <chr> <chr> <list>
## 1 A crystallized <NULL>
## 2 C crystallized <NULL>
## 3 E crystallized <NULL>
## 4 N crystallized <NULL>
## 5 O crystallized <NULL>4.2 Step 2: Run Models and Extract Results
Once the data are prepped, we are ready to begin running models! For Method 2: Pooled One Stage Models with Study-Specific Effects, we will run two variants. The first mimics a dummy coding procedure often seen in the economics literature where study-specific effects are captured by including studies as fixed effects. The second is a multilevel model in which particiants are nested within studies, and studies are treated as random effects, rather than fixed effects, as in Method 2A. Both of of use when study-specific effects are of interest, but Method 2B is of particular use when wanting to mimic heterogeneity estimates and meta regression techniques of traditional meta analyses.
4.2.1 Method 2A: Pooled One Stage Models with Dummy Codes
4.2.1.1 Analytic Plan
In the present study, we estimate associations between the Big Five personality traits and crystallized cognitive abilities using a one-stage, fully pooled regression approach with effects coded contrasts for study and individual participant data. The procedure is as follows:
To estimate sample-specific effects using regression contrasts, there are a variety of coding procedures that could be used. We used effects codes, which result in a term that captures the overall estimated effect and k-1 terms capturing sample-specific deviations from the overall estimate (all sample estimates can then be recovered via linear combinations). Using effects codes allows us to statistically test whether each sample (minus one) differs from the overall effect. Although dummy codes are most frequently used in regression-based approaches, effects code are the means by which regression is linked to analysis of variance in order to estimate marginal means. For example, if there were only three samples under investigation the matrix of effect codes would look like:
\[ \begin{array}{rr} d_1 & d_2 \\ .5 & 0 \\ 0 & .5 \\ -.5 & -.5 \end{array} \],
where each column \(j\) is an indicator in the model and each row \(k\) represents a different sample. The first coefficient would capture the difference between the first sample and the overall estimate, while the second coefficient would capture the difference between the second sample and the overall estimate. By additionally adding an interaction between the effects coded terms and the key term (i.e. personality traits in the present study), we can additionally test for sample-specific estimates of key associations in the regression. 2. Statistical Modeling. Models were run separately for each personality characteristic, outcome, covariate, and moderator (10; 3 participant-level) combination. Using the R programming language, we created a series of functions to (1) set up and run the model and extract model coefficients and (2) extract overall simple-effects predictions for participant-level moderator models (i.e. predicted values across levels of the moderator values). In addition, to estimate sample-specific effects using effects codes, we created additional functions to (3) set up the effects codes and ensure interpretable variable labels, (4) set up and run a series of linear combinations of the effect-coded terms to get personality-crystallized abilities associations or moderators of their association for each sample, and (5) extract sample-specific simple-effects prediction for the participant-level moderators. Each of these functions were designed to flexibly handle different classes of covariates and moderators (e.g., nominal, numeric) and to wrangle the results into easily combined and understandable data frames. The functions and more detail on them can be found in the online materials. The basic form of the model is as follows:
\[Y_{ij} = b_0 + b_1\ast predictor_{ij}+b2\ast study{\ 1}_{ij}+\ldots+b_k\ast study{k}_{ij}+ b_{k+1}\ast predictor_{ij}\ast study\ 1_{ij}+\ldots+b_{2\ast k}\ast predictor_{ij}\ast study{k}_{ij}+\epsilon_{ij}\] \(\varepsilon_{ij}\sim\mathcal{N}(0, \sigma^2)\),
where \(k\) indicates the number of samples - 1. Of interest are two key sets of terms. \(b_1\) indicates the average personality-cognition relationship, and \(b_{k+1}\) to \(b_{2k}\) represent effect coded sample-specific differences in outcome associations (i.e. the estimate for a sample is \(b_1+b_{2k}\)). All other terms capture sample-specific differences in overall cognitive ability levels, if any. There is no measure of cross-sample heterogeneity, and residuals are assumed to be homogeneous.
Models were tested using the base R lm() function in R, and the tidy() function from the broom package to extract model coefficients and confidence intervals (CI). Inferences were made based on the 95% confidence intervals. Overall simple-effects predictions were calculated by providing the full range of personality (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. For sample-specific effects, effects coded linear combinations for each sample were provided to the glht() function from the multcomp package (version 1.4.20; Hothorn et al., 2008).
4.2.1.2 Model Function
The first thing we need is a function that will bring in the data, create a formula in the model based on input on the type (Frequentist or Bayesian), moderators (none, age, gender, SRhealth, and education), and combinations of covariates (single or fully adjusted based on age, gender, SRhealth, and education). This case is slightly complicated given that we are using dummy codes for each study. We want to get an overall estimate as well as study-specific estimates. So we’ll have to build contrasts that give us these. Then we run the model, extract its fixed and study-specific effect estimates, and save both for later. By saving the results, it will make it easier and faster for us to extract the necessary model results later while still retaining all information from the original model.
ipd2a_mod_fun <- function(trait, outcome, type, mod, cov){
## load the data
load(sprintf("%s/data/one_stage/%s_%s.RData", local_path, trait, outcome))
## Applt effects codes
d <- contr_fun(d)
## compiled Bayesian model to speed up processing and avoid crashing
if(type == "Bayesian") load(sprintf("%s/results/2a_ipd_dc/bayes_sample_mod.RData", local_path))
# get the model formula
## model formula
if (cov == "all") cv <- c("age", "gender", "education")
if (!cov %in% c("all", "none")) cv <- cov
rhs <- c("p_value", "study", "p_value:study")
rhs <- if(cov != "none") c(rhs, cv) else rhs
if(!mod %in% c("none", stdyModers$short_name)){rhs <- c(rhs, paste("p_value", mod, "study", sep = "*"))}
if(mod %in% stdyModers$short_name){rhs <- c(rhs, mod, 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, iter = 2000
, warmup = 1000)}
closeAllConnections()
save(m, file = sprintf("%s/results/2a_ipd_dc/%s/models/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
## extract model terms and confidence intervals & save
fx <- tidy(m, conf.int = T) %>%
select(term, estimate, conf.low, conf.high)
rx <- if(mod != stdyModers$short_name) std_eff_fun(m, type, mod) else NA
save(fx, rx, file = sprintf("%s/results/2a_ipd_dc/%s/summary/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
if(mod != "none"){
load(sprintf("%s/results/2a_ipd_dc/%s/models/%s_%s_%s_%s.RData", local_path, type, outcome, trait, mod, cov))
pred.fx <- ipd2a_simpeff_fun(m, mod, type)
pred.rx <- if(mod %in% stdyModers$short_name) NULL else ipd2a_rx_pred_fun(m, mod, type)
save(pred.fx, pred.rx, file = sprintf("%s/results/2a_ipd_dc/%s/predicted/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
}
## clean up the local function environment
rm(list = c("d", "f", "rhs", "m", "fx", "rx"))
gc()
}4.2.1.3 Contrasts Function
One function we will call in the modeling function for Method 2A is one that creates the contrasts that we’ll use for the models we’ll run. Because the contrasts necessary to get study-specific estimates and an overall estimate are not orthogonal, we will necessarily have to use two steps to get all these estimates. As a result, the method we use for contrasts doesn’t matter that much. We’ll use effects coding, as we pre-registered, because it will give us the overall estimate in the main model, as well as estimates of how S - 1 studies differ from that overall estimate.
Basically, teh function below will:
1. Get names of all those studies.
2. Use the contr.sum() function to create effects codes (i.e. a S-1 x S matrix of codes)
3. Set the row and column names to the studies we are examining
4. Save that contrast matrix to the study column of teh data frame
4.2.1.4 Study-Specific Effects Function
As noted previously, once we run the model, we will have to use a second step to get the study-specific estimates for all studies. To make this flexible for both Frequentist and Bayesian estimates, overall personality-outcome associations, and moderator associations is somewhat complicated. So we’ll start by getting names of the terms from the model. Then we’ll find the key term (either personality-outcome association or moderator effect) and the studies linked to these terms. For those studies, then getting the contrast formula is straightforward (overall effect + study effect). For the study without a model term, the formula is more complicated (overall effect - (S-1 study effects)).
Once we have these, they can easily be fed to the glht() function in the multcomp package (frequentist) or the hypothesis() function in the brms package (bayesian) to get all study-specific estimates and confidence / certainty intervals.
Before we do this, though, let’s try to better understand where these terms come from. Consider the case of four samples used to estimate an effect. This would require 3 effects coded contrast terms (\(d_1\) to \(d_3\))in the model coded as follows:
\[ \begin{array}{rr} d_1 & d_2 & d_3 \\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ -1 & -1 & -1 \end{array} \],
where \(d_1\) to \(d_{k-1}\) indicates the difference between samples 1 to \(k-1\) and the grand mean. Thus, to recover the study-specific estimates, we can return to the equation for the model for the four sample case:
\(Y_{ij}= b_0 + b_1*predictor_{ij} + b2*study1_{ij} + b3*study2_{ij} + b4*study3_{ij} + b_{5}*predictor_{ij}*study1_{ij} + b_{6}*predictor_{ij}*study2_{ij} + b_{7}*predictor_{ij}*study3_{ij} + \epsilon_{ij}\),
where \(b_5\) to b_7 indicates the difference in the persoanlity-cognitive ability association between the average association across studies and samples 1-3, respectively. Then, to recover the estimate for samples 1-3, we simply need to sum the term for the overall association with the term for the difference in the association between that and samples 1-3, or:
Sample 1: \(b_1 + b_5\)
Sample 2: \(b_1 + b_6\)
Sample 3: \(b_1 + b_7\)
Lastly, to recover the estiamte for sample 4, we must subtract the estimates of differences for samples 1-3 from the grand mean, or:
Sample 4: \(b_1 - b_5 - b_6 - b_7\)
The function below extracts each of these and saves the estimates and their confidence intervals according to study names for easy combination and understanding.
std_eff_fun <- function(m, type, mod){
if(type == "Bayesian"){
std <- as.character(unique(m$data$study)) # vector of studies
nc <- nrow(fixef(m)); nr <- length(std) # number of rows and columns
trms <- rownames(fixef(m))
} else{
std <- as.character(unique(m$model$study)) # vector of studies
nc <- length(coef(m)); nr <- length(std) # number of rows and columns
trms <- names(coef(m))
}
modtrm <- if(mod %in% c("none", stdyModers$short_name)) "p_value" else trms[grepl(paste0("p_value:", mod), trms)]
stdtrms <- trms[grepl("p_value:study", trms)]
stdtrms <- if (mod %in% c("none", stdyModers$short_name)) stdtrms else stdtrms[grepl(mod, stdtrms)]
# create character contrasts
cntrm <- paste(stdtrms, "+", modtrm, " = 0")
cntrm <- c(cntrm, paste0(" - ", stdtrms, collapse = "") %>%
paste(modtrm, ., collapse = "") %>%
paste(., "= 0", collapse = ""))
names(cntrm) <- std
# Run the contrasts and rename to match overall model coefficient tables
h <- if(type == "Bayesian") {
hypothesis(m, cntrm)$hypothesis %>% # brms hypothesis function
select(study = Hypothesis, estimate = Estimate,
conf.low = CI.Lower, conf.high = CI.Upper) %>%
mutate(term = modtrm) %>%
as_tibble()
} else {
(multcomp::glht(m, cntrm) %>% # multcomp hypothesis function
confint(., calpha = multcomp::univariate_calpha()))$confint %>%
data.frame() %>%
mutate(study = std,
term = modtrm) %>%
rename(estimate = Estimate, conf.low = lwr, conf.high = upr) %>%
as_tibble()
}
return(h)
}4.2.1.5 Simple Effects Function
4.2.1.5.1 Overall Effects
fx_lm_pred_fun <- function(m, newdata){
tt <- terms(m)
Terms <- delete.response(tt)
mf <- model.frame(Terms, newdata, xlev = m$xlevels)
X <- model.matrix(Terms, mf, contrasts.arg = m$contrasts)
X[,grepl("study", colnames(X))] <- 0
p <- m$rank
p1 <- seq_len(p)
piv <- m$qr$pivot[p1]
beta <- m$coefficients
predictor <- drop(X[, piv, drop = FALSE] %*% beta[piv])
w <- m$weights
res.var <- {
r <- m$residuals
rss <- sum(r^2)
df <- m$df.residual
rss/df
}
XRinv <- X[, piv] %*% qr.solve(qr.R(m$qr)[p1, p1])
ip <- drop(XRinv^2 %*% rep(res.var, p))
tfrac <- qt((.05)/2, df)
hwid <- tfrac * sqrt(ip)
predictor <- cbind(predictor, predictor + hwid %o% c(1, -1))
colnames(predictor) <- c("fit", "lwr", "upr")
return(predictor)
}
fx_bayes_pred_fun <- function(m, newdata){
ps <- posterior_samples(m, "^b", as.matrix = T)
X <- prepare_predictions(
m, newdata %>% mutate(o_value = 1)
)$dpars$mu$fe$X[,1:ncol(ps), drop = FALSE]
X[,grepl("study", colnames(X))] <- 0
predictor <- X %*% t(ps)
predictor <- posterior_summary(t(predictor))
return(predictor)
}
ipd2a_simpeff_fun <- function(m, moder, type){
d <- if(type == "Bayesian") m$data else m$model
std_lev <- rev(levels(d$study))[1]
d <- d %>% select(-o_value, -p_value, -study)
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)) %>%
mutate(study = std_lev)
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 %>% select(-study),
fx_bayes_pred_fun(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 %>% select(-study),
fx_lm_pred_fun(m, newdata = mod_frame) %>% 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)
}4.2.1.5.2 Study-Specific Effects
# short function to get crossings of continuous personality and moderator levels
# crossing_fun <- function(df, mod, mod_lev){
# pred.rx <- crossing(
# p_value = seq(0,10,.5),
# mod_value = mod_lev
# ) %>%
# setNames(c("p_value", mod))
# return(pred.rx)
# }
ipd2a_rx_pred_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 %>%
group_by(study) %>%
select_if(is.numeric) %>%
pivot_longer(-study
, names_to = "item"
, values_to = "value") %>%
group_by(study, item) %>%
summarize_at(vars(value), lst(mean, sd), na.rm = T) %>%
ungroup()
}
if(any(sapply(d, class) == "factor")){
fct_lev <- d %>%
group_by(study) %>%
select_if(is.factor) %>%
summarize_all(~list(unique(.))) %>%
ungroup() %>%
data.frame()
}
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 <- if(md_cl == "numeric") {
crossing(
p_value = seq(0,10,.5),
modvalue = md_levs,
study = unique(d$study)
) %>% setNames(c("p_value", moder, "study"))#%>% full_join(
# msd %>%
# filter(item == moder) %>%
# mutate(lower = mean - sd, upper = mean + sd) %>%
# select(-sd) %>%
# pivot_longer(cols = c(mean, lower, upper)
# , names_to = "meas"
# , values_to = "modvalue") %>%
# pivot_wider(names_from = "item", values_from = "modvalue") %>%
# select(study, one_of(moder))
# )
} else {
crossing(
p_value = seq(0,10,.5)
, mod_value = md_levs#unique(fct_lev[,moder][[1]])
, study = unique(d$study)
) %>%
setNames(c("p_value", moder, "study"))
}
if(ncol(d) > 0){
if(any(sapply(d, class) == "numeric")){
mod_frame <- d %>%
group_by(study) %>%
select_if(is.numeric) %>%
summarize_all(mean, na.rm = T) %>%
ungroup() %>%
full_join(mod_frame)
}
if(any(sapply(d, class) == "factor")){
mod_frame <- d %>%
group_by(study) %>%
select_if(is.factor) %>%
summarize_all(~levels(.)) %>%
ungroup() %>%
full_join(mod_frame)
}
}
pred.rx <- 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.rx)
}4.2.1.6 Run Models and Summaries
The code below generates all possible preregistrered combinations of traits, outcomes, types, moderators, and covariates. Then these combinations are fed serially to the model function written previously, which will run and save the results.
load(sprintf("%s/data/one_stage/E_crystallized.RData", local_path))
# clean data & keep only needed columns and a subset of the used variables
d <- d %>%
group_by(study) %>%
nest() %>%
ungroup() %>%
mutate(data = map(data, ~(.) %>% filter(row_number() %in% sample(1:nrow(.), 100, replace = F)))) %>%
unnest(data)
d <- contr_fun(d)
cv <- c("age", "gender", "education")
rhs <- c("p_value", "study", "p_value:study")
rhs <- c(rhs, cv)
rhs <- paste(rhs, collapse = " + ")
f <- paste("o_value ~ ", rhs, collapse = "")
# set priors & model specifications
Prior <- c(set_prior("student_t(3, 0, 2)", class = "b"),
set_prior("student_t(3, 0, 5)", class = "Intercept"))
Iter <- 30; Warmup <- 21; treedepth <- 20
f <- bf(f)
m <- brm(formula = f
, data = d
, prior = Prior
, iter = Iter
, warmup = Warmup
, cores = 4)
save(m, file = sprintf("%s/results/2a_ipd_dc/bayes_sample_mod.RData", local_path))
rm(list = c("d", "Prior", "Iter", "Warmup", "treedepth", "f", "m", "rhs", "cv"))# plan(multisession(workers = 8L))
nested_ipd2a_reg <-
# moderator result combinations
crossing(
Trait = traits$short_name
, Outcome = outcomes$short_name
, type = c("Frequentist", "Bayesian")
, Moderator = c("age", "gender", "education")
, Covariate = c("none", "all")
) %>%
full_join(
# personality-outcome association results across covariates
crossing(
Trait = traits$short_name
, Outcome = outcomes$short_name
, type = c("Frequentist", "Bayesian")
, Moderator = c("none", unique(stdyModers$short_name))
, Covariate = c("none", "age", "gender", "education", "all")
)
) %>%
# full_join(done) %>% filter(is.na(done)) %>%
# filter(Trait == "N") %>%
filter(type == "Bayesian") %>%
mutate(run =
# future_pmap(list(Trait, Outcome, type, Moderator, Covariate)
pmap(list(Trait, Outcome, type, Moderator, Covariate)
, possibly(ipd2a_mod_fun, "uh oh")
# , .progress = T
# , .options = furrr_options(
# globals = c("ipd2a_mod_fun"
# , "std_eff_fun"
# , "fx_lm_pred_fun"
# , "fx_bayes_pred_fun"
# , "ipd2a_simpeff_fun"
# , "crossing_fun"
# , "contr_fun"
# , "ipd2a_rx_pred_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()
nested_ipd2a_reg %>%
filter(!Moderator %in% stdyModers$short_name) %>%
write.table(.
, file = sprintf("%s/scripts/cluster/args/bayesian/ipd2a_bayesian.txt", local_path)
, row.names = F)4.2.1.7 Compile Results
Once all the models are run, we are ready to compile all their results. By saving the fixed and study-level effects results previously, we are able to simply load those results and ignore the models. However, because we also saved the models, we can also recall and extract information from them if and when needed.
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/2a_ipd_dc/%s/%s/%s", local_path, type, folder, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
n_fun <- function(fileName, type){
m <- loadRData(fileName, type, "^m", "models")
d <- if(type == "Bayesian") m$data else m$model
n <- d %>% group_by(study) %>% tally() %>% ungroup()
return(n)
}
## load in "fixed" effects
## first get file names
nested_ipd2a_reg <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/2a_ipd_dc/%s/models", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate"), sep = "_", remove = F) %>%
filter(!Moderator %in% stdyModers$short_name) %>%
## read in the files
mutate(Covariate = str_remove(Covariate, ".RData"),
fx = map2(file, type, ~loadRData(.x, .y, "fx", "summary")),
rx = map2(file, type, ~loadRData(.x, .y, "rx", "summary")),
n = map2(file, type, n_fun)) %>%
select(-file) 4.2.1.7.1 Tables
Next, we want to format the study results in APA table format. In this case, we are interested in the fixed and study-specific effects of personality predicting cognitive ability when there were no moderators, and the personality x moderator interaction when there was a moderator. We’ll anticipate a need to present both just fixed effects as well as fixed and study-specific effects by creating tables for each.
First, let’s format the data.
ipd2a_reg_tab <-
### fixed effects
nested_ipd2a_reg %>%
select(-rx, -n) %>%
unnest(fx) %>% # unnesting
# keep key terms
filter((Moderator == "none" & term == "p_value") |
(Moderator != "none" & grepl("p_value:", term) & !grepl("p_value:study", term))) %>%
mutate(study = "Overall") %>%
### study specific effects
full_join(
nested_ipd2a_reg %>%
select(-fx, -n) %>%
unnest(rx) %>% # unnesting
filter((Moderator == "none" & term == "p_value") |
(Moderator != "none" & grepl("p_value:", term) & !grepl("p_value:study", term))) %>%
mutate( )
) %>%
filter(!is.na(term)) %>%
# reformatting: mark significance, prettify Trait, covariate, and moderator 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_remove_all(term, "p_value:"),
term = mapvalues(term, c("scaleBFIMS", "scaleIPIPNEO", "scaleTDAM40", "countryTheNetherlands")
, c("scaleBFI-S", "scaleIPIP NEO", "scaleTDA-40", "countryThe Netherlands")),
term2 = factor(term, c(moders$short_term, stdyModers$short_term),
c(moders$long_term, stdyModers$long_term))) %>%
# prettify the number format
mutate_at(vars(estimate, conf.low, conf.high),
~ifelse(abs(.) < .001, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
# combine the effects, bold significance, factor and label study-specfic effects
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),
study = factor(study, levels = c(studies_long, "Overall"),
labels = c(studies_long, "Overall"))) %>%
# reshaping: remove extra columns, arrange by key variables, and make wide
select(type, Outcome, Trait, Moderator, Covariate, study, term2, est) %>%
arrange(type, Outcome, Trait, Moderator, Covariate, study, term2) %>%
pivot_wider(names_from = "Trait", values_from = "est")
# filter(row_number() %in% c(9, 10))
# spread(Trait, est)
ipd2a_reg_tab## # A tibble: 242 × 11
## type Outcome Moderator Covariate study term2 E A C N O
## <chr> <fct> <fct> <fct> <fct> <fct> <list> <list> <list> <list> <list>
## 1 Bayesian Crystallized Ability None Unadjusted EAS Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## 2 Bayesian Crystallized Ability None Unadjusted GSOEP Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## 3 Bayesian Crystallized Ability None Unadjusted HILDA Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## 4 Bayesian Crystallized Ability None Unadjusted HRS Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## 5 Bayesian Crystallized Ability None Unadjusted MAP Personality <chr [1]> <NULL> <chr [1]> <chr [1]> <NULL>
## 6 Bayesian Crystallized Ability None Unadjusted ROS Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## 7 Bayesian Crystallized Ability None Unadjusted SATSA Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## 8 Bayesian Crystallized Ability None Unadjusted Overall Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## 9 Bayesian Crystallized Ability None Unadjusted <NA> Personality <chr [2]> <NULL> <NULL> <chr [2]> <chr [1]>
## 10 Bayesian Crystallized Ability None Age EAS Personality <chr [1]> <chr [1]> <chr [1]> <chr [1]> <chr [1]>
## # ℹ 232 more rowsipd2a_res <- nested_ipd2a_reg %>%
select(-rx, -n) %>%
unnest(fx) %>% # unnesting
# keep key terms
filter((Moderator == "none" & term == "p_value") |
(Moderator != "none" & grepl("p_value:", term) & !grepl("p_value:study", term))) %>%
mutate(study = "Overall") %>%
mutate_at(vars(estimate, conf.low, conf.high),
~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
mutate(est = sprintf("\\textit{b} = %s, 95\\%% CI [%s, %s]", estimate, conf.low, conf.high)) Now that we’ve formatted the values, we can group by moderators and save results as separate tables. Even though additional information could be included given that we have one outcome, we’ll stick with this split because it will make it easier for those using this tutorial who multiple traits, outcomes, covariates, and moderators.
4.2.1.7.1.1 Fixed Effects
## table function
ipd2a_fx_tab_fun <- function(d, type, moder){
md <- mapvalues(moder, moders$long_name, moders$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$term2)) == 1) rep(1,6) else c(2, rep(1,5))
names(cs) <- c(" ", paste0("<strong>", traits$long_name, "</strong>"))
cln <- if(length(unique(d$term2)) == 1) c("Covariates", rep("<em>b</em> [CI]", 5)) else c("Covariates", "Term", rep("<em>b</em> [CI]", 5))
# cln <- if(length(unique(d$term2)) == 1) c("Covariates", rep("\\textit{b} [CI]", 5)) else c("Covariates", "Term", rep("\\textit{b} [CI]", 5))
al <- if(length(unique(d$term2)) == 1) c("r", rep("c", 5)) else c("r", "r", rep("c", 5))
if(length(unique(d$term2)) == 1) d <- d %>% select(-term2)
cap <- if(md == "none") "2A Pooled One Stage Models with Dummy Codes: Overall Effects of Personality-Crystallized Domain Associations" else sprintf("2A Pooled One Stage Models with Dummy Codes: Overall %s Moderation of Personality-Crystallized Domain Associations", 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") %>%
# kable_styling(full_width = F, font_size = 7) %>%
add_header_above(cs, escape = F)
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/2a_ipd_dc/%s/tables/overall/%s.html"
, local_path, type, md))
return(tab)
}
ipd2a_fx_tab <- ipd2a_reg_tab %>%
filter(study == "Overall" & !Moderator %in% stdyModers$long_name) %>%
select(-study) %>%
group_by(type, Moderator) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Moderator), ipd2a_fx_tab_fun))
ipd2a_fx_tab## # A tibble: 8 × 4
## type Moderator data tab
## <chr> <fct> <list> <list>
## 1 Bayesian None <tibble [5 × 8]> <kablExtr [1]>
## 2 Bayesian Age <tibble [2 × 8]> <kablExtr [1]>
## 3 Bayesian Gender <tibble [2 × 8]> <kablExtr [1]>
## 4 Bayesian Education <tibble [2 × 8]> <kablExtr [1]>
## 5 Frequentist None <tibble [5 × 8]> <kablExtr [1]>
## 6 Frequentist Age <tibble [2 × 8]> <kablExtr [1]>
## 7 Frequentist Gender <tibble [2 × 8]> <kablExtr [1]>
## 8 Frequentist Education <tibble [2 × 8]> <kablExtr [1]># save(ipd2a_fx_tab, ipd2a_res, file = sprintf("%s/manuscript/results/ipd2a_fx_tab.RData", res_path))
## Frequentist
(ipd2a_fx_tab %>% filter(Moderator == "None" & type == "Frequentist"))$tab[[1]]|
Extraversion
|
Agreeableness
|
Conscientiousness
|
Neuroticism
|
Openness to Experience
|
|
|---|---|---|---|---|---|
| Covariates | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| Crystallized Ability | |||||
| Unadjusted |
0.05 [0.02, 0.08] |
-0.00 [-0.03, 0.03] |
0.04 [0.02, 0.07] |
-0.13 [-0.15, -0.11] |
0.29 [0.25, 0.32] |
| Age |
0.05 [0.02, 0.08] |
0.01 [-0.02, 0.04] |
0.04 [0.01, 0.07] |
-0.13 [-0.16, -0.11] |
0.30 [0.26, 0.33] |
| Gender |
0.05 [0.02, 0.08] |
0.01 [-0.02, 0.04] |
0.04 [0.01, 0.07] |
-0.13 [-0.15, -0.11] |
0.29 [0.26, 0.33] |
| Education |
0.02 [-0.01, 0.05] |
-0.001 [-0.03, 0.03] |
0.01 [-0.02, 0.04] |
-0.07 [-0.09, -0.04] |
0.17 [0.14, 0.21] |
| Fully Adjusted |
0.02 [-0.01, 0.05] |
-0.01 [-0.04, 0.02] |
0.01 [-0.02, 0.04] |
-0.07 [-0.09, -0.05] |
0.17 [0.14, 0.21] |
|
Extraversion
|
Agreeableness
|
Conscientiousness
|
Neuroticism
|
Openness to Experience
|
|
|---|---|---|---|---|---|
| Covariates | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| Crystallized Ability | |||||
| Unadjusted |
0.05 [0.02, 0.08] |
-0.00 [-0.03, 0.02] |
0.04 [0.02, 0.06] |
-0.13 [-0.15, -0.11] |
0.29 [0.25, 0.32] |
| Age |
0.05 [0.02, 0.08] |
0.01 [-0.02, 0.04] |
0.04 [0.01, 0.07] |
-0.13 [-0.16, -0.11] |
0.30 [0.26, 0.33] |
| Gender |
0.05 [0.02, 0.08] |
0.01 [-0.02, 0.04] |
0.04 [0.01, 0.07] |
-0.13 [-0.15, -0.11] |
0.29 [0.26, 0.33] |
| Education |
0.02 [-0.01, 0.05] |
-0.00 [-0.03, 0.03] |
0.01 [-0.02, 0.04] |
-0.07 [-0.09, -0.04] |
0.17 [0.14, 0.21] |
| Fully Adjusted |
0.02 [-0.01, 0.05] |
-0.01 [-0.04, 0.02] |
0.01 [-0.02, 0.04] |
-0.07 [-0.09, -0.05] |
0.17 [0.14, 0.21] |
ipd2a_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(" ", paste0("<strong>", traits$short_name, "</strong>"))
# 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("2A Pooled One Stage Models with Dummy Codes: Fixed Effect Estimates of %s Personality-Crystallized Domain Associations", cov)
# kable the table
tab <- d %>%
select(-Moderator, -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, escape = F)
# 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/2a_ipd_dc/%s/tables/key terms/%s.html"
, local_path, type, covar))
return(tab) # return the html table
}
ipd2a_fx_tab2 <- ipd2a_reg_tab %>%
filter(study == "Overall") %>%
arrange(Moderator, term2) %>%
filter(Covariate %in% c("Unadjusted", "Fully Adjusted")) %>%
group_by(Outcome, type, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Covariate), ipd2a_tab_fun))
## Frequentist, no moderator
(ipd2a_fx_tab2 %>% filter(type == "Frequentist" & Covariate == "Fully Adjusted"))$tab[[1]]|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Term | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| None | |||||
| Personality |
0.02 [-0.01, 0.05] |
-0.01 [-0.04, 0.02] |
0.01 [-0.02, 0.04] |
-0.07 [-0.09, -0.05] |
0.17 [0.14, 0.21] |
| Age | |||||
| Age |
-0.00 [-0.01, 0.00] |
-0.000 [-0.00, 0.00] |
-0.00 [-0.01, 0.001] |
0.00 [-0.000, 0.01] |
0.000 [-0.00, 0.01] |
| Gender | |||||
| Gender (Male v Female) |
0.06 [0.001, 0.12] |
0.02 [-0.04, 0.08] |
-0.00 [-0.06, 0.05] |
-0.04 [-0.09, 0.01] |
-0.01 [-0.08, 0.06] |
| Education | |||||
| Education (Years) |
-0.00 [-0.01, 0.01] |
-0.00 [-0.02, 0.01] |
0.01 [-0.001, 0.02] |
-0.00 [-0.01, 0.01] |
-0.00 [-0.02, 0.01] |
|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Term | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| None | |||||
| Personality |
0.02 [-0.01, 0.05] |
-0.01 [-0.04, 0.02] |
0.01 [-0.02, 0.04] |
-0.07 [-0.09, -0.05] |
0.17 [0.14, 0.21] |
| Age | |||||
| Age |
-0.00 [-0.01, 0.00] |
-0.000 [-0.00, 0.00] |
-0.00 [-0.01, 0.001] |
0.00 [-0.001, 0.01] |
0.000 [-0.00, 0.01] |
| Gender | |||||
| Gender (Male v Female) |
0.06 [0.00, 0.12] |
0.02 [-0.04, 0.08] |
-0.00 [-0.06, 0.05] |
-0.04 [-0.09, 0.01] |
-0.01 [-0.08, 0.07] |
| Education | |||||
| Education (Years) |
-0.00 [-0.01, 0.01] |
-0.00 [-0.02, 0.01] |
0.01 [-0.001, 0.02] |
-0.00 [-0.01, 0.01] |
-0.00 [-0.02, 0.01] |
4.2.1.7.1.2 Study-Specific Effects
## table function
ipd2a_rx_tab_fun <- function(d, type, moder){
md <- mapvalues(moder, moders$long_name, moders$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$term2)) == 1) c(2, rep(1,5)) else c(3, rep(1,5))
names(cs) <- c(" ", traits$short_name)
cln <- if(length(unique(d$term2)) == 1) c("Covariates", "Study", rep("<em>b</em> [CI]", 5)) else c("Covariates", "Study", "Term", rep("<em>b</em> [CI]", 5))
# cln <- if(length(unique(d$term2)) == 1) c("Covariates", rep("\\textit{b} [CI]", 5)) else c("Covariates", "Term", rep("\\textit{b} [CI]", 5))
al <- if(length(unique(d$term2)) == 1) c("r", "r", rep("c", 5)) else c("r", "r", "r", rep("c", 5))
if(length(unique(d$term2)) == 1) d <- d %>% select(-term2)
# caption
cap <- if(md == "none") "2A Pooled One Stage Models with Dummy Codes: Overall Effects of Personality-Crystallized Domain Associations" else sprintf("2A Pooled One Stage Models with Dummy Codes: Overall %s Moderation of Personality-Crystallized Domain Associations", moder)
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") %>%
collapse_rows(1, "top") %>%
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/2a_ipd_dc/%s/tables/study-specific/%s.html"
, local_path, type, md))
return(tab)
}
ipd2a_rx_tab <- ipd2a_reg_tab %>%
filter(study != "Overall") %>%
arrange(type, Outcome, Moderator, Covariate, study, term2) %>%
group_by(type, Moderator) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Moderator), ipd2a_rx_tab_fun))
ipd2a_fx_tab## # A tibble: 8 × 4
## type Moderator data tab
## <chr> <fct> <list> <list>
## 1 Bayesian None <tibble [5 × 8]> <kablExtr [1]>
## 2 Bayesian Age <tibble [2 × 8]> <kablExtr [1]>
## 3 Bayesian Gender <tibble [2 × 8]> <kablExtr [1]>
## 4 Bayesian Education <tibble [2 × 8]> <kablExtr [1]>
## 5 Frequentist None <tibble [5 × 8]> <kablExtr [1]>
## 6 Frequentist Age <tibble [2 × 8]> <kablExtr [1]>
## 7 Frequentist Gender <tibble [2 × 8]> <kablExtr [1]>
## 8 Frequentist Education <tibble [2 × 8]> <kablExtr [1]>|
E
|
A
|
C
|
N
|
O
|
||
|---|---|---|---|---|---|---|
| Covariates | Study | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| Crystallized Ability | ||||||
| Unadjusted | EAS |
0.13 [0.08, 0.19] |
0.02 [-0.03, 0.07] |
0.02 [-0.03, 0.07] |
-0.13 [-0.18, -0.08] |
0.37 [0.32, 0.42] |
| GSOEP |
0.08 [0.03, 0.13] |
-0.02 [-0.07, 0.03] |
-0.03 [-0.08, 0.03] |
-0.06 [-0.11, -0.02] |
0.10 [0.05, 0.14] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.05 [0.02, 0.08] |
0.07 [0.04, 0.09] |
-0.15 [-0.18, -0.13] |
0.27 [0.25, 0.30] |
|
| HRS |
-0.00 [-0.03, 0.02] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.07 [-0.10, -0.05] |
0.21 [0.19, 0.24] |
|
| LASA | NULL | NULL | NULL |
-0.13 [-0.17, -0.09] |
NULL | |
| MAP |
0.03 [-0.03, 0.09] |
NULL |
0.10 [0.02, 0.19] |
-0.16 [-0.22, -0.09] |
NULL | |
| MARS | NULL | NULL | NULL |
-0.28 [-0.39, -0.18] |
NULL | |
| ROS |
-0.01 [-0.07, 0.06] |
0.12 [0.04, 0.20] |
0.13 [0.05, 0.21] |
-0.09 [-0.15, -0.02] |
0.18 [0.10, 0.27] |
|
| SATSA |
0.04 [-0.05, 0.14] |
-0.23 [-0.35, -0.12] |
-0.17 [-0.27, -0.07] |
-0.09 [-0.19, -0.00] |
0.44 [0.30, 0.58] |
|
| Age | EAS |
0.12 [0.03, 0.21] |
0.09 [0.01, 0.18] |
0.01 [-0.08, 0.09] |
-0.13 [-0.21, -0.05] |
0.42 [0.34, 0.51] |
| GSOEP |
0.08 [0.03, 0.13] |
-0.02 [-0.07, 0.03] |
-0.03 [-0.08, 0.02] |
-0.06 [-0.11, -0.02] |
0.10 [0.05, 0.14] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.04 [0.01, 0.07] |
0.07 [0.04, 0.09] |
-0.16 [-0.19, -0.14] |
0.28 [0.25, 0.30] |
|
| HRS |
-0.00 [-0.03, 0.02] |
0.06 [0.03, 0.09] |
0.16 [0.14, 0.19] |
-0.08 [-0.10, -0.05] |
0.22 [0.19, 0.24] |
|
| LASA | NULL | NULL | NULL |
-0.12 [-0.17, -0.08] |
NULL | |
| MAP |
0.03 [-0.03, 0.09] |
NULL |
0.10 [0.02, 0.19] |
-0.16 [-0.22, -0.09] |
NULL | |
| MARS | NULL | NULL | NULL |
-0.28 [-0.38, -0.18] |
NULL | |
| ROS |
-0.01 [-0.08, 0.06] |
0.12 [0.04, 0.20] |
0.13 [0.06, 0.21] |
-0.09 [-0.16, -0.02] |
0.19 [0.11, 0.27] |
|
| SATSA |
0.04 [-0.05, 0.14] |
-0.23 [-0.35, -0.12] |
-0.17 [-0.27, -0.07] |
-0.09 [-0.18, -0.00] |
0.45 [0.30, 0.59] |
|
| Gender | EAS |
0.12 [0.03, 0.21] |
0.10 [0.01, 0.18] |
0.01 [-0.08, 0.09] |
-0.13 [-0.22, -0.05] |
0.42 [0.34, 0.51] |
| GSOEP |
0.08 [0.03, 0.13] |
-0.02 [-0.07, 0.03] |
-0.03 [-0.08, 0.02] |
-0.06 [-0.11, -0.02] |
0.10 [0.05, 0.14] |
|
| HILDA |
0.01 [-0.02, 0.04] |
0.05 [0.02, 0.08] |
0.07 [0.04, 0.09] |
-0.15 [-0.18, -0.13] |
0.27 [0.25, 0.30] |
|
| HRS |
-0.00 [-0.03, 0.02] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.07 [-0.10, -0.05] |
0.21 [0.19, 0.24] |
|
| LASA | NULL | NULL | NULL |
-0.13 [-0.17, -0.09] |
NULL | |
| MAP |
0.03 [-0.03, 0.09] |
NULL |
0.10 [0.02, 0.19] |
-0.16 [-0.22, -0.10] |
NULL | |
| MARS | NULL | NULL | NULL |
-0.28 [-0.39, -0.18] |
NULL | |
| ROS |
-0.01 [-0.08, 0.06] |
0.12 [0.04, 0.20] |
0.13 [0.06, 0.21] |
-0.09 [-0.16, -0.02] |
0.19 [0.11, 0.27] |
|
| SATSA |
0.05 [-0.05, 0.14] |
-0.23 [-0.35, -0.12] |
-0.17 [-0.27, -0.07] |
-0.09 [-0.18, -0.00] |
0.44 [0.30, 0.58] |
|
| Education | EAS |
0.07 [-0.01, 0.16] |
0.05 [-0.03, 0.13] |
-0.04 [-0.13, 0.04] |
-0.06 [-0.14, 0.02] |
0.19 [0.11, 0.27] |
| GSOEP |
0.09 [0.01, 0.17] |
0.03 [-0.06, 0.12] |
0.01 [-0.08, 0.10] |
-0.03 [-0.10, 0.05] |
0.06 [-0.02, 0.13] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.04 [0.01, 0.07] |
0.03 [0.01, 0.06] |
-0.14 [-0.16, -0.11] |
0.19 [0.17, 0.22] |
|
| HRS |
-0.03 [-0.06, -0.01] |
0.04 [0.01, 0.06] |
0.09 [0.06, 0.11] |
-0.02 [-0.04, -0.00] |
0.10 [0.08, 0.13] |
|
| LASA | NULL | NULL | NULL |
-0.05 [-0.09, -0.01] |
NULL | |
| MAP |
0.01 [-0.05, 0.07] |
NULL |
0.06 [-0.02, 0.14] |
-0.03 [-0.09, 0.03] |
NULL | |
| MARS | NULL | NULL | NULL |
-0.13 [-0.23, -0.03] |
NULL | |
| ROS |
-0.08 [-0.14, -0.01] |
0.03 [-0.04, 0.11] |
0.08 [0.00, 0.15] |
-0.01 [-0.07, 0.06] |
-0.01 [-0.09, 0.07] |
|
| SATSA |
0.02 [-0.07, 0.12] |
-0.19 [-0.30, -0.08] |
-0.14 [-0.24, -0.05] |
-0.07 [-0.15, 0.02] |
0.31 [0.17, 0.45] |
|
| Fully Adjusted | EAS |
0.07 [-0.01, 0.16] |
0.05 [-0.03, 0.13] |
-0.04 [-0.13, 0.04] |
-0.06 [-0.14, 0.02] |
0.19 [0.11, 0.27] |
| GSOEP |
0.08 [0.01, 0.16] |
0.02 [-0.06, 0.11] |
0.00 [-0.09, 0.10] |
-0.03 [-0.10, 0.05] |
0.06 [-0.01, 0.13] |
|
| HILDA |
0.01 [-0.01, 0.03] |
0.03 [-0.00, 0.05] |
0.03 [0.00, 0.05] |
-0.14 [-0.16, -0.11] |
0.20 [0.17, 0.22] |
|
| HRS |
-0.03 [-0.06, -0.01] |
0.03 [0.00, 0.06] |
0.08 [0.06, 0.11] |
-0.02 [-0.05, -0.00] |
0.10 [0.08, 0.13] |
|
| LASA | NULL | NULL | NULL |
-0.05 [-0.09, -0.01] |
NULL | |
| MAP |
0.01 [-0.05, 0.07] |
NULL |
0.06 [-0.02, 0.14] |
-0.03 [-0.09, 0.03] |
NULL | |
| MARS | NULL | NULL | NULL |
-0.13 [-0.23, -0.03] |
NULL | |
| ROS |
-0.08 [-0.14, -0.01] |
0.04 [-0.04, 0.11] |
0.08 [0.01, 0.15] |
-0.00 [-0.07, 0.06] |
-0.01 [-0.08, 0.07] |
|
| SATSA |
0.03 [-0.07, 0.12] |
-0.20 [-0.31, -0.09] |
-0.14 [-0.24, -0.05] |
-0.07 [-0.16, 0.02] |
0.31 [0.17, 0.45] |
|
4.2.1.7.1.3 Heterogeneity Estimates
This header is here to simply emphasize that this method does not provide heterogeneity estimates because it does not provide study-specific estimates, unlike Methods 2b, 3, and 4.
4.2.1.7.1.4 All Model Terms
ipd2a_mod_tab <- nested_ipd2a_reg %>%
select(-rx, -n) %>%
unnest(fx) %>%
# 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))) %>%
# 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) %>%
arrange(type, Outcome, Trait, Moderator, Covariate, term) %>%
pivot_wider(names_from = "Trait", values_from = "est")
ipd2a_mod_tab_fun <- function(d, type, out, moder, cov){
print(paste(type, out, moder, cov))
md <- mapvalues(moder, moders$long_name, moders$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(" ", paste0("<strong>", traits$long_name, "</strong>"))
# 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") "2A Pooled One Stage Models with Dummy Codes: All Model Estimates of Fixed Effect Personality-Crystallized Domain Associations" else sprintf("2A Pooled One Stage Models with Dummy Codes: All Model Estimates of Fixed Effect %s Moderation of Personality-Crystallized Domain Associations", md)
# kable the table
tab <- d %>%
arrange(term) %>%
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, escape = F)
# save the resulting html table
save_kable(tab, file = sprintf("%s/results/2a_ipd_dc/%s/tables/all terms/%s-%s-%s.html"
, local_path, type, o, md, cv))
return(tab) # return the html table
}
ipd2a_mod_tab <- ipd2a_mod_tab %>%
filter(!Moderator %in% stdyModers$long_name) %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Outcome, Moderator, Covariate), ipd2a_mod_tab_fun))## [1] "Bayesian Crystallized Ability None Unadjusted"
## [1] "Bayesian Crystallized Ability None Age"
## [1] "Bayesian Crystallized Ability None Gender"
## [1] "Bayesian Crystallized Ability None Education"
## [1] "Bayesian Crystallized Ability None Fully Adjusted"
## [1] "Bayesian Crystallized Ability Age Unadjusted"
## [1] "Bayesian Crystallized Ability Age Fully Adjusted"
## [1] "Bayesian Crystallized Ability Gender Unadjusted"
## [1] "Bayesian Crystallized Ability Gender Fully Adjusted"
## [1] "Bayesian Crystallized Ability Education Unadjusted"
## [1] "Bayesian Crystallized Ability Education Fully Adjusted"
## [1] "Frequentist Crystallized Ability None Unadjusted"
## [1] "Frequentist Crystallized Ability None Age"
## [1] "Frequentist Crystallized Ability None Gender"
## [1] "Frequentist Crystallized Ability None Education"
## [1] "Frequentist Crystallized Ability None Fully Adjusted"
## [1] "Frequentist Crystallized Ability Age Unadjusted"
## [1] "Frequentist Crystallized Ability Age Fully Adjusted"
## [1] "Frequentist Crystallized Ability Gender Unadjusted"
## [1] "Frequentist Crystallized Ability Gender Fully Adjusted"
## [1] "Frequentist Crystallized Ability Education Unadjusted"
## [1] "Frequentist Crystallized Ability Education Fully Adjusted"4.2.1.7.2 Figures
4.2.1.7.2.1 Overall Forest
ipd2a_fx_plot_fun <- function(df, mod, type, cov){
print(paste(mod, type, cov))
m <- mapvalues(mod, moders$short_name, moders$long_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", m)}
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.8, 1.3)) +
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 2A: Pooled Regression Using Dummy Codes") +
guides(shape = "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/2a_ipd_dc/%s/figures/overall forest/%s_%s_fixed.png", local_path, type, mod, cov), width = local_path*2, height = ht+ht2)
rm(p)
gc()
return(T)
}
nested_ipd2a_reg %>%
filter(!Moderator %in% stdyModers$short_name) %>%
select(-rx) %>%
unnest(fx) %>%
filter((Moderator == "none" & term == "p_value")|
(Moderator != "none" & grepl("^p_value:", term) & !grepl("study", 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)),
Outcome = forcats::fct_rev(Outcome)) %>%
group_by(type, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(pmap(list(data, Moderator, type, Covariate), possibly(ipd2a_fx_plot_fun, NA_real_)))4.2.1.7.2.2 Study-Specific Forest
ipd2a_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 %>% full_join(tibble(study = " ", estimate = NA, n = NA))
df <- df %>% arrange(estimate)
stds <- df$study[!df$study %in% c(" Overall", " ")]
df <- df %>%
mutate(study = factor(study, rev(c(" ", stds, " Overall")))
# , 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", "no")
, ub = ifelse(conf.high > lim[2], "upper", "no")
, 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 == " Overall", "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(ub == "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)
, axis = "tblr"
, 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)
gc()
save(p
, file = sprintf("%s/results/2a_ipd_dc/%s/figures/study specific forest/rdata/%s_%s_%s_%s.RData", local_path, type, outcome, trait, mod, cov))
return(p)
}
## fixed effects
nested_ipd2a_reg_fp <- nested_ipd2a_reg %>%
filter(Moderator %in% moders$short_name) %>%
select(-rx, -n) %>%
unnest(fx) %>%
mutate(study = " Overall") %>%
## random effects
full_join(
nested_ipd2a_reg %>%
filter(Moderator %in% moders$short_name) %>%
mutate(rx = map2(rx, n, ~(.x) %>% full_join(.y))) %>%
select(-fx, -n) %>%
unnest(rx) %>%
mutate(study = mapvalues(study, c("OCTOTWIN", "BASEI"), c("OCTO-Twin", "BASE"))
, study = mapvalues(study, studies_long, studies_sp, warn_missing = F))
# mutate(term = ifelse(Moderator != "none", paste(term, mapvalues(Moderator, moders$short_name, moders$short_term, warn_missing = F), sep = ":"), term))
) %>%
## 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")) %>%
## grouping for plotting
group_by(Outcome, Moderator, type, Covariate, Trait) %>%
nest() %>%
ungroup() %>%
mutate(p = pmap(list(data, Outcome, Moderator, type, Covariate, Trait), ipd2a_rx_plot_fun))
ipd2a_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 2A: Pooled Regression Using Dummy Codes",
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/2a_ipd_dc/%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/2a_ipd_dc/%s/figures/study specific forest/%s_%s_%s.pdf", local_path, type, outcome, mod, cov)
, width = 10, height = 10)
return(T)
}
nested_ipd2a_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), ipd2a_rx_plot_comb_fun))4.2.1.7.2.3 Overall Simple Effects
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/2a_ipd_dc/%s/%s/%s", local_path, type, folder, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
## load in "fixed" effects
## first get file names
nested_ipd2a_simp <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/2a_ipd_dc/%s/predicted", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate"), sep = "_", remove = F) %>%
## read in the files
mutate(Covariate = str_remove(Covariate, ".RData"),
pred.fx = map2(file, type, ~loadRData(.x, .y, "pred.fx", "predicted")),
pred.rx = map2(file, type, ~loadRData(.x, .y, "pred.rx", "predicted"))) %>%
select(-file)
nested_ipd2a_simpipd2a_se_plot_fun <- function(df, outcome, mod, type, cov){
print(paste(outcome, mod))
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)
d <- round(max(abs(min(df$pred)), abs(max(df$pred))), 3)
# mini <- if(d > 2) .05 else 0-(d+(d/5))
# maxi <- if(d > 2) 2.05 else 0+d+(d/5)
# lim <- c(mini, maxi)
# brk <- if(d > 2) c(0, 1, 2) else{round(c(0-d-(d/10), 0, 0+d+(d/10)),2)}
# lab <- if(d > 2){c("0", "1", "2")} else{str_remove(c(round(0-d-(d/10),2), 0, round(0+d+(d/10),2)), "^0")}
titl <- if(mod == "none"){o} else {sprintf("%s: Personality x %s Simple Effects", o, m)}
# colnames(df)[colnames(df) == mod] <- "mod_value"
df <- df %>% unclass %>% data.frame
df$mod_value <- df[,mod]
df <- df %>% select(-all_of(mod)) %>% as_tibble
if(class(df$mod_value) == "factor"){df <- df %>% mutate(mod_fac = factor(mod_value))}
else{df <- df %>%
group_by(Trait) %>%
mutate(mod_fac = factor(mod_value, levels = unique(mod_value), labels = c("-1 SD", "M", "+1 SD"))) %>%
ungroup()
}
lt <- c("dotted", "solid", "dashed")[1:length(unique(df$mod_fac))]
df %>%
mutate(Trait = factor(Trait, levels = traits$short_name, labels = traits$long_name),
lower = ifelse(lower < 4, 4, lower),
upper = ifelse(upper > 10, 10, upper)) %>%
ggplot(aes(x = p_value
, y = pred
, group = mod_fac)) +
scale_y_continuous(limits = c(4,10)
, breaks = c(4, 6, 8, 10)
, labels = c(4, 6, 8, 10)) +
scale_linetype_manual(values = lt) +
geom_ribbon(aes(ymin = lower
, ymax = upper
, fill = mod_fac)
, alpha = .25) +
geom_line(aes(linetype = mod_fac)) +
labs(x = "Personality (POMP)"
, y = paste(o, "(POMP)")
, title = titl
, linetype = m
, fill = m
, subtitle = "Method 2A: Pooled Regression Using Dummy Codes") +
facet_wrap(~Trait, nrow = 3) +
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/2a_ipd_dc/%s/figures/overall simple effects/%s_%s_%s.png", local_path, type, outcome, mod, cov), width = 6, height = 6)
}
ipd2a_se_plot <- nested_ipd2a_simp %>%
select(-pred.rx) %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(data = map(data, ~(.) %>% unnest(pred.fx)),
plot = pmap(list(data, Outcome, Moderator, type, Covariate), ipd2a_se_plot_fun))4.2.1.7.2.4 Study-Specific Simple Effects
ipd2a_std_se_plot_fun <- function(df, outcome, trait, mod, cov, type){
print(paste(outcome, mod))
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)
d <- round(max(abs(min(df$pred)), abs(max(df$pred))), 3)
# mini <- if(d > 2) .05 else 0-(d+(d/5))
# maxi <- if(d > 2) 2.05 else 0+d+(d/5)
# lim <- c(mini, maxi)
# brk <- if(d > 2) c(0, 1, 2) else{round(c(0-d-(d/10), 0, 0+d+(d/10)),2)}
# lab <- if(d > 2){c("0", "1", "2")} else{str_remove(c(round(0-d-(d/10),2), 0, round(0+d+(d/10),2)), "^0")}
titl <- if(mod == "none"){sprintf("%s: %s", o, trt)} else {sprintf("%s: %s x %s Simple Effects", o, trt, m)}
# colnames(df)[colnames(df) == mod] <- "mod_value"
df <- df %>% mutate(study = mapvalues(study, c("BASEI", "OCTOTWIN"), c("BASE", "OCTO-Twin")))
df <- df %>% unclass %>% data.frame
df$mod_value <- df[,mod]
df <- df %>% select(-all_of(mod)) %>% as_tibble
if(class(df$mod_value) == "factor"){df <- df %>% mutate(mod_fac = factor(mod_value))}
else{df <- df %>%
group_by(study) %>%
mutate(mod_fac = factor(mod_value, levels = unique(mod_value), labels = c("-1 SD", "M", "+1 SD"))) %>%
ungroup()
}
# lt <- c("dotted", "solid", "dashed")[1:length(unique(df$mod_fac))]
std <- unique(df$study)
cols <- (stdcolors %>% filter(studies %in% std))$colors
lt <- (stdcolors %>% filter(studies %in% std))$lt
ht <- length(unique(df$mod_fac))
df %>%
as_tibble() %>%
mutate(gr = ifelse(study == "Overall", "Overall", "study")) %>%
group_by(p_value, mod_fac, study, gr) %>%
summarize_at(vars(pred, lower, upper), mean, na.rm = T) %>%
mutate(study = factor(study
, levels = c("Overall", studies_long),
, labels = c("Overall", studies_long)),
lower = ifelse(lower < 4, 4, lower),
upper = ifelse(upper > 10, 10, upper)) %>%
ggplot(aes(x = p_value
, y = pred
, group = study)) +
scale_y_continuous(limits = c(4,10)
, breaks = c(4, 6, 8, 10)
, labels = c(4, 6, 8, 10)) +
scale_linetype_manual(values = lt) +
scale_color_manual(values = cols) +
scale_size_manual(values = c(2,.8)) +
# geom_ribbon(aes(ymin = lower
# , ymax = upper
# , fill = study)
# , alpha = .25) +
geom_line(aes(linetype = study, color = study, size = gr)) +
labs(x = paste(trt, "(POMP)")
, y = paste(o, "(POMP)")
, title = titl
, linetype = "Study"
, fill = "Study"
, color = "Study"
, subtitle = "Method 2A: Pooled Regression Using Dummy Codes") +
guides(size = "none") +
facet_wrap(~mod_fac, nrow = 1) +
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/2a_ipd_dc/%s/figures/study specific simple effects/%s_%s_%s_%s.png", local_path, type, outcome, trt, mod, cov), width = 3*ht, height = 5)
}
nested_ipd2a_simp %>%
filter(!Moderator %in% stdyModers$short_name) %>%
mutate(pred.fx = map(pred.fx, ~(.) %>% mutate(study = "Overall")),
comb.fx = map2(pred.fx, pred.rx, full_join)) %>%
select(-pred.fx, -pred.rx) %>%
mutate(pmap(list(comb.fx, Outcome, Trait, Moderator, Covariate, type), possibly(ipd2a_std_se_plot_fun, NA_real_)))4.2.1.8 Sample Results Section
We examined estimates of overall prospective associations between Big Five personality characteristics and crystallized abilities as well as participant-level moderators of those associations using fully pooled, one-stage regression models with effects codes to capture sample-specific estimates using individual participant data in 11 studies. Heterogeneity was indirectly estimated by extracting sample-level estimates of associations and examining dispersion measures.
Table S3 presents the fully adjusted estimates of the key terms from all prospective Big Five personality characteristic-crystallized abilities associations and participant-level moderators of those associations. Across the 20 key terms, four (20%) were significant. For unmoderated personality-cognition associations, neuroticism (-) and openness (+) were significantly associated with later crystallized abilities. In other words, higher neuroticism was prospectively associated with lower crystallized abilities across samples, while higher openness was prospectively associated with higher crystallized abilities across samples.
|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Term | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| None | |||||
| Personality |
0.02 [-0.01, 0.05] |
-0.01 [-0.04, 0.02] |
0.01 [-0.02, 0.04] |
-0.07 [-0.09, -0.05] |
0.17 [0.14, 0.21] |
| Age | |||||
| Age |
-0.00 [-0.01, 0.00] |
-0.000 [-0.00, 0.00] |
-0.00 [-0.01, 0.001] |
0.00 [-0.000, 0.01] |
0.000 [-0.00, 0.01] |
| Gender | |||||
| Gender (Male v Female) |
0.06 [0.001, 0.12] |
0.02 [-0.04, 0.08] |
-0.00 [-0.06, 0.05] |
-0.04 [-0.09, 0.01] |
-0.01 [-0.08, 0.06] |
| Education | |||||
| Education (Years) |
-0.00 [-0.01, 0.01] |
-0.00 [-0.02, 0.01] |
0.01 [-0.001, 0.02] |
-0.00 [-0.01, 0.01] |
-0.00 [-0.02, 0.01] |
In addition to the overall effects, we also examined those same effects at the sample level. Figure S3 presents the forest plot of fully adjusted sample-specific and overall prospective associations between the Big Five and crystallized abilities across all methods. Forest plots for all participant-level moderators are available in the online materials and the web app. As shown in the forest plot, openness was most consistently associated with crystallized ability, with five of seven samples showing a positive and significant prospective association between openness and crystallized abilities. Other associations, such as those with extraversion and conscientiousness, were less consistent, with at least one significant association in different directions across samples. For example, although there was no overall association between conscientiousness and crystallized abilities, SATSA showed a significant negative association, and HILDA, ROS, and HRS showed positive associations.
Figure 4.1: Figure S3. 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.
Finally, we also examined the consistency of simple effects across samples for each moderator. Figure S4 displays the fully adjusted predicted crystallized abilities levels at different levels of extraversion across samples for those who identified as male or female. The pattern suggests a positive association for females and null association for males. This holds for females in EAS (b = 0.22, 95% CI [0.09, 0.34]) and GSOEP (b = 0.14, 95% CI [0.02, 0.25]), while all other samples demonstrate null associations. For men, the associations seem somewhat less consistent, with three samples demonstrating significant negative associations (HILDA, b = -0.04, 95% CI [-0.08, -0.005]; HRS, b = -0.04, 95% CI [-0.08, -0.006]; ROS, b = -0.08, 95% CI [-0.16, -0.008]) one sample demonstrating a positive association (EAS, b = 0.16, 95% CI [0.03, 0.29]), and four samples demonstrating null associations (GSOEP, b = 0.08, 95% CI [-0.04, 0.19]; MAP, b = 0.009, 95% CI [-0.06, 0.08]; OCTO-Twin, b = 0.08, 95% CI [-0.07, 0.22]; BASE, b = -0.14, 95% CI [-0.37, 0.10]).
Figure 4.2: Figure S4. 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.
4.2.2 Method 2B: Pooled One Stage Models with Random Effects
For Method 2B, we’ll be estimating overall and study-specific effects using multilevel models in which participants are nested within studies. In a basic sense, these differ from using dummy codes in that those models are estimated using OLS, while multilevel models are estimated using maximum likelihood. In addition, unlike OLS, linear models, ML MLM treats nesting units as a distribution centered around the overall effect and shrinks random effects toward 0.
4.2.2.1 Analytic Plan
In the present study, we estimate associations between the Big Five personality traits and crystallized cognitive abilities using a one-stage, fully pooled, multilevel regression approach with participants nested within sample using individual participant data. The procedure is as follows:
4.2.2.1.1 Statistical Modeling
Models are estimated separately for each personality characteristic, outcome, covariate, and moderator combination. Because these multilevel models provide both sample-specific estimates and estimates of cross-sample heterogeneity, we created a series of functions in the R programming language to (1) set up and run the model and extract model coefficients and (2) extract overall simple-effects predictions for participant- and sample-level moderator models (i.e. predicted values across levels of the moderator values), (3) extract the sample-specific effects, (4) extract sample-specific simple-effects prediction for the participant-level moderators, and (5) extract estimates of heterogeneity. Each of these functions were designed to flexibly handle different classes of covariates and moderators (e.g., nominal, numeric) and to wrangle the results into easily combined and understandable data frames. The functions and more detail on them can be found in the online materials.
The simple unadjusted personality characteristic-crystallized abilities model is as follows:
Level 1 \(Y_{ij}=\beta_{0J}+\beta_{1J}*predictor_{ij}+\epsilon_{ij}\)
Level 2 \(\beta_{0j}=\gamma_{00}+u_{0j}\) \(\beta_{1J}=\gamma_{10}+u_{1j}\)
\(\epsilon_{ij}\sim\mathcal{N}\left(0,\sigma^2\right)\) \(\begin{matrix}u_{0j}\\\ u_{1j}\\\end{matrix}\ \sim\mathcal{N}\left(\begin{matrix}\tau_{00}^2\ \ &\tau_{10}\\\ \tau_{10}&\tau_{11}^2\\\end{matrix}\ \right)\),
where the prospective associations for each sample will be captured by \(\beta_{1j}\), \(j\) indicates each sample 1 to \(j\), and \(u_{1j}\) captures the sample specific deviation from the overall estimate of the personality-cognition relationship \(\gamma_{10}\). The multilevel models decompose the residual variance into different sources, in this case, participant-level residual variance (\(\sigma^2\)) and sample-level residual variance (\(\tau_{00}^2\) and \(\tau_{11}^2\)). The models also use partial pooling, which pools information across samples and shrinks sample-level estimates toward the fixed effect. Participant-level moderators will be added at Level 1 with random effects at Level 2. For example, two new Level 1 terms were added for baseline age (centered at 60): \(\beta_{2J}\), which indicates the sample-specific estimate of the main effect of age, and \(\beta_{3J}\), which indicates the difference in the sample-specific estimate of the prospective personality-outcome association as a function of age. For each, we will include random slopes (\(u_{2j}\) and \(u_{3j}\)) to allow these to vary across samples. The overall main effect of the moderator and the moderator estimate will be captured by \(\gamma_{20}\) and \(\gamma_{30}\), respectively.
Heterogeneity estimates are captured in \(\tau\) matrix, where \(\tau_{00}^2\) captures the variance in the random intercepts (i.e. variance in average levels of cognitive ability across samples), \(\tau_{11}^2\) captures the variance in random slope (i.e. variance in the personality-cognitive ability association across samples), and \(\tau_{01}\) captures the correlation between average levels of cognitive ability and personality-cognitive ability associations across samples (e.g., is the association stronger in samples with higher average crystallized ability). Moderators, but not covariates, were included as random slopes and allowed to covary with other random terms.
Models were tested using the lme4 package in R, and the tidy() function from the broom.mixed package to extract model coefficients and confidence intervals (CI) of fixed effects and heterogeneity estimates. . In order to extract the random effects (i.e. Level 2 variances) and interval estimates, the effects argument was set to "ran_pars," the conf.int argument was set to TRUE, the conf.method argument was set to "boot", and the n.sim argument was set to 500. Inferences were made based on the 95% bootstrapped confidence intervals. As before, overall simple-effects predictions were calculated by providing the full range of personality (0-10) and average levels of the covariates from the data used to estimate the model as the "newdata" argument, and setting the "re.form" argument as FALSE in the lme4 predict.merMod() function in conjunction with the bootMer() function to provide intervals around the estimates. For sample-specific effects, averages and levels of covariates were determined using data from each sample separately and provided to the lme4 predict.merMod() function as "newdata."
4.2.2.2 Model Function
The first thing we need is a function that will bring in the data, create a formula in the model based on input on the type (Frequentist or Bayesian), moderators (none, age, gender, and education), and combinations of covariates (single or fully adjusted based on age, gender, and education). This case is slightly less complicated than the case of dummy codes, because we simply include study as the nesting unit and personality effect of the personality and moderator effect as random slopes. This results in a intercept and target effect for each study without the need to build contrasts. Then we run the model, extract its fixed and study-specific effect estimates, and save both for later. By saving the results, it will make it easier and faster for us to extract the necessary model results later while still retaining all information from the original model.
ipd2b_mod_fun <- function(trait, outcome, type, mod, cov){
## load the data
load(sprintf("%s/data/one_stage/%s_%s.RData", local_path, trait, outcome))
## compiled Bayesian model to speed up processing and avoid crashing
if(type == "Bayesian") load(sprintf("%s/results/2b_ipd_mlm/bayes_sample_mod.RData", local_path))
## 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 = "*"))}
re <- if(mod == "none" | mod %in% stdyModers$short_name) "(p_value | study)" else paste(paste("(p_value", mod, sep = " * "), "| study)")
rhs <- paste(c(rhs, re), collapse = " + ")
f <- paste("o_value ~ ", rhs, collapse = "")
## run the models & save
m <- if(type == "Frequentist"){lmer(f, data = d)} else {update(m, formula = f, newdata = d, iter = 2000
, warmup = 1000)}
save(m, file = sprintf("%s/results/2b_ipd_mlm/%s/models/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
## extract model terms and confidence intervals & save
fx <- tidy(m, conf.int = T) %>%
select(term, estimate, conf.low, conf.high)
rx <- std_eff_fun(m, type)
save(fx, rx, file = sprintf("%s/results/2b_ipd_mlm/%s/summary/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
## extract heterogeneity estimates
het <- ipd2b_hetero_fun(m, type)
save(het, file = sprintf("%s/results/2b_ipd_mlm/%s/heterogeneity/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
## simple effects for moderators
if(mod != "none"){
pred.fx <- ipd2b_fx_pred_fun(m, mod, type)
pred.rx <- if(mod %in% stdyModers$short_name) NULL else ipd2b_rx_pred_fun(m, mod, type)
save(pred.fx, pred.rx, file = sprintf("%s/results/2b_ipd_mlm/%s/predicted/%s_%s_%s_%s.RData"
, local_path, type, outcome, trait, mod, cov))
}
## clean up the local function environment
rm(list = c("d", "f", "rhs", "m", "fx", "rx", "het"))
gc()
}4.2.2.3 Study-Specific Effects Function
As noted previously, once we run the model, we will have to use a second step to get the study-specific estimates for all studies. Unlike with dummy codes, doing so is much more straightforward. We just have to pull study-specific effects using the coef() for both Bayesian and Frequentist approaches. However, to get confidence intervals for Frequentist approaches, we will additionally have to use the standard_error() function from the parameters package to get standard errors.
std_eff_fun <- function(m, type){
if(type == "Frequentist"){
# coef function gives fixed + random effect estimates but not SE's
# so we'll take those estimates and get their SE's from the parameters package
coef(m)$study %>%
data.frame() %>%
rownames_to_column("study") %>%
mutate(term = "estimate") %>%
full_join(
parameters::standard_error(m, effects = "random")$study %>%
data.frame() %>%
rownames_to_column("study") %>%
mutate(term = "SE")) %>%
rename(Intercept = X.Intercept.) %>%
# select(study, term, , p_value) %>%
pivot_longer(c(-study, -term), names_to = "names", values_to = "estimate") %>%
pivot_wider(names_from = "term", values_from = "estimate") %>%
rename(term = names) %>%
mutate(conf.low = estimate - 2*SE, conf.high = estimate + 2*SE,
term = ifelse(grepl("p_value.", term), str_replace_all(term, "p_value.", "p_value:"), term))
} else {
coef(m, probs = c(0.025, 0.975))[[1]] %>% array_tree(3) %>%
tibble(term = names(.), data = .) %>%
mutate(data = map(data, ~(.) %>% data.frame %>%
rownames_to_column("study"))) %>%
unnest(data) %>%
select(term, study, estimate = Estimate, conf.low = Q2.5, conf.high = Q97.5)
}
}4.2.2.4 Heterogeneity Estimates Function
The Final pieces of information we need to extract from these models are estimates of the heterogeneity of effects across studies.
ipd2b_hetero_fun <- function(m, type){
args <- list(x = m, effects = "ran_pars", conf.int = T)
if(type == "Frequentist")args$conf.method = "boot"; args$nsim <- 100; args$parallel = "multicore"; args$ncpus = 4
do.call(tidy, args) %>%
select(group, term, estimate, conf.low, conf.high) %>%
separate(term, c("est", "term"), sep = "__") %>%
mutate_at(vars(estimate:conf.high), ~ifelse(est == "sd", .^2, .)) %>%
mutate(est = ifelse(est == "sd", "var", est))
}4.2.2.5 Simple Effects Function
4.2.2.5.1 Fixed Effects
predIntlme4 <- function(m, mod_frame, ref){
b <- bootMer(m
, FUN = function(x) lme4:::predict.merMod(x, newdata = mod_frame , re.form = ref)
, nsim = 100
, parallel = "multicore"
, ncpus = 8
)
ci <- apply(b$t, 2, quantile, probs = c(.05/2, 1 - .05/2)) %>% t()
data.frame(pred = b$t0, ci) %>%
setNames(c("pred", "lower", "upper")) %>% as_tibble()
}
ipd2b_fx_pred_fun <- function(m, moder, type){
d <- if(type == "Bayesian") m$data else m@frame
d <- d %>% select(-o_value, -p_value, -study)
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_vars <- colnames(d)[sapply(d, class) == "factor"]
fct_lev <- d %>%
select_if(is.factor) %>%
summarize_all(~list(levels(.)))
}
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 {
(fct_lev %>% select(one_of(moder)) %>% unnest(one_of(moder)) %>% data.frame())[,moder]
}
mod_frame <- expand.grid(
p_value = seq(0,10,.5)
, modvalue = md_levs
, stringsAsFactors = F
) %>% 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
, re_formula = NA) %>% data.frame
) %>%
select(one_of(colnames(m$data)), pred = Estimate, lower = Q2.5, upper = Q97.5)
} else {
bind_cols(
mod_frame,
predIntlme4(m, mod_frame, NA)
) #%>%
# select(one_of(colnames(m@frame)), pred, lower, upper)
}
rm(list = c("m", "mod_frame", "d", "md_levs"))
gc()
return(pred.fx)
}4.2.2.5.2 Study-Specific Effects
ipd2b_rx_pred_fun <- function(m, moder, type){
d <- if(type == "Bayesian") m$data else m@frame
d <- d %>% select(-o_value, -p_value)
cols <- colnames(d)
md_cl <- class(d[,moder])
if(any(sapply(d, class) == "numeric")){
msd <- d %>%
group_by(study) %>%
select_if(is.numeric) %>%
pivot_longer(-study
, names_to = "item"
, values_to = "value") %>%
group_by(study, item) %>%
summarize_at(vars(value), lst(mean, sd), na.rm = T) %>%
ungroup()
}
if(any(sapply(d, class) == "factor")){
fct_vars <- colnames(d)[sapply(d %>% select(-study), class) == "factor"]
fct_lev <- d %>%
group_by(study) %>%
select_if(is.factor) %>%
summarize_all(~list(unique(.))) %>%
ungroup() #%>%
# unnest(one_of(fct_vars))
}
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 {
(fct_lev %>% select(one_of(moder)) %>% unnest(one_of(moder)) %>% data.frame())[,moder]
}
mod_frame <- if(md_cl == "numeric") {
crossing(
p_value = seq(0,10,.5),
modvalue = md_levs,
study = unique(d$study)
) %>% setNames(c("p_value", moder, "study"))#%>% full_join(
# msd %>%
# filter(item == moder) %>%
# mutate(lower = mean - sd, upper = mean + sd) %>%
# select(-sd) %>%
# pivot_longer(cols = c(mean, lower, upper)
# , names_to = "meas"
# , values_to = "modvalue") %>%
# pivot_wider(names_from = "item", values_from = "modvalue") %>%
# select(study, one_of(moder))
# )
} else {
crossing(
p_value = seq(0,10,.5)
, mod_value = md_levs#unique(fct_lev[,moder][[1]])
, study = unique(d$study)
) %>%
setNames(c("p_value", moder, "study"))
}
if(ncol(d) > 0){
if(any(sapply(d, class) == "numeric")){
mod_frame <- d %>%
group_by(study) %>%
select_if(is.numeric) %>%
summarize_all(mean, na.rm = T) %>%
ungroup() %>%
full_join(mod_frame)
}
if(any(sapply(d, class) == "factor")){
mod_frame <- d %>%
group_by(study) %>%
select_if(is.factor) %>%
summarize_all(levels) %>%
ungroup() %>%
full_join(mod_frame)
}
}
pred.rx <- 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,
tibble(pred = predict(m, newdata = mod_frame))#predIntlme4(m, mod_frame, NULL)
)
}
rm(list = c("m", "mod_frame", "d"))
gc()
return(pred.rx)
}4.2.2.6 Run Models and Summaries
# Sample Bayesian Model
# load data
load(sprintf("%s/data/one_stage/N_crystallized.RData", local_path))
# clean data & keep only needed columns and a subset of the used variables
d <- d %>%
group_by(study) %>%
nest() %>%
ungroup() %>%
mutate(data = map(data, ~(.) %>% filter(row_number() %in% sample(1:nrow(.), 100, replace = F)))) %>%
unnest(data)
# set priors & model specifications
Prior <- c(set_prior("cauchy(0,1)", class = "sd"),
set_prior("student_t(3, 0, 2)", class = "b"),
set_prior("student_t(3, 0, 5)", class = "Intercept"))
Iter <- 30; Warmup <- 21; treedepth <- 20
f <- formula(o_value ~ p_value*age + gender + education + (age*p_value | study))
m <- brm(formula = f
, data = d
, prior = Prior
, iter = Iter
, warmup = Warmup
, cores = 4)
save(m, file = sprintf("%s/results/2b_ipd_mlm/bayes_sample_mod.RData", local_path))
rm(list = c("d", "Prior", "Iter", "Warmup", "treedepth", "f", "m"))# done <- tibble(type = c("Bayesian", "Frequentist"),
# file = map(type, ~list.files(sprintf("%s/results/2b_ipd_mlm/%s/predicted", local_path, .)))) %>%
# unnest(file) %>%
# separate(file, c("Outcome", "Trait", "Moderator", "Covariate"), sep = "_") %>%
# mutate(Covariate = str_remove_all(Covariate, ".RData")
# , done = "done") %>%
# filter(!is.na(Covariate))
nested_ipd_reg <- crossing(
Trait = traits$short_name
, Outcome = outcomes$short_name
, type = c("Frequentist", "Bayesian")
, Moderator = c("age", "gender", "education")
, Covariate = c("none", "all")
) %>%
full_join(
crossing(
Trait = traits$short_name
, Outcome = outcomes$short_name
, type = c("Frequentist", "Bayesian")
, Moderator = c("none", stdyModers$short_name)
, Covariate = c("none", "age", "gender", "education", "all")
)
) %>%
# full_join(done) %>% filter(is.na(done) & type != "Frequentist") %>%
filter(Moderator %in% c("age", "education", "baseAge", "baseYear", "predInt")) %>%
filter(Trait == "N" & type == "Frequentist")
mutate(run = pmap(list(Trait, Outcome, type, Moderator, Covariate), ipd2b_mod_fun))
nested_ipd_reg %>%
write.table(.
, file = sprintf("%s/scripts/cluster/args/frequentist/ipd2b_frequentist_pred.txt", local_path)
, row.names = F)4.2.2.7 Compile Results
Once all the models are run, we are ready to compile all their results. By saving the fixed and study-level effects results previously, we are able to simply load those results and ignore the models. However, because we also saved the models, we can also recall and extract information from them if and when needed.
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/2b_ipd_mlm/%s/%s/%s", local_path, type, folder, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
n_fun <- function(fileName, type){
m <- loadRData(fileName, type, "^m", "models")
d <- if(type == "Bayesian") m$data else m@frame
n <- d %>% group_by(study) %>% tally() %>% ungroup()
return(n)
}
fix_terms <- function(rx){
if(!"term" %in% colnames(rx)) rx <- rx %>% mutate(term = names)
return(rx)
}
## load in "fixed" effects
## first get file names
nested_ipd2b_reg <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/2b_ipd_mlm/%s/summary", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate"), sep = "_", remove = F) %>%
## read in the files
mutate(Covariate = str_remove(Covariate, ".RData"),
fx = map2(file, type, ~loadRData(.x, .y, "fx", "summary")),
rx = map2(file, type, possibly(~loadRData(.x, .y, "rx", "summary"), NA_real_)),
n = map2(file, type, n_fun)) %>%
select(-file) %>%
filter(!is.na(rx)) %>%
mutate(rx = map(rx, fix_terms))4.2.2.7.1 Tables
Next, we want to format the study results in APA table format. In this case, we are interested in the fixed and study-specific effects of personality predicting cognitive ability when there were no moderators, and the personality x moderator interaction when there was a moderator. We’ll anticipate a need to present both just fixed effects as well as fixed and study-specific effects by creating tables for each.
First, let’s format the data.
ipd2b_reg_tab <-
### fixed effects
nested_ipd2b_reg %>%
select(-rx) %>%
unnest(fx) %>% # unnesting
# keep key terms
filter((Moderator == "none" & term == "p_value") |
(Moderator != "none" & grepl("p_value:", term)) &
!grepl("cor", term) & !grepl("sd", term)) %>%
mutate(study = "Overall") %>%
### study specific effects
full_join(
nested_ipd2b_reg %>%
select(-fx) %>%
unnest(rx) %>% # unnesting
mutate(term = ifelse(is.na(term), names, term)) %>%
select(-names) %>%
filter((Moderator == "none" & term == "p_value") |
(Moderator != "none" & grepl("p_value:", term)))
) %>%
# reformatting: mark significance, prettify Trait, covariate, and moderator 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_remove_all(term, "p_value:"),
term = mapvalues(term, c("scaleBFIMS", "scaleIPIPNEO", "scaleTDAM40", "countryTheNetherlands",
"scaleBFI.S", "scaleIPIP.NEO", "scaleTDA.40", "countryThe.Netherlands")
, c("scaleBFI-S", "scaleIPIP NEO", "scaleTDA-40", "countryThe Netherlands",
"scaleBFI-S", "scaleIPIP NEO", "scaleTDA-40", "countryThe Netherlands")),
term = factor(term, c(moders$short_term, stdyModers$short_term),
c(moders$long_term, stdyModers$long_term))) %>%
# prettify the number format
mutate_at(vars(estimate, conf.low, conf.high),
~ifelse(abs(.) < .0014, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
# combine the effects, bold significance, factor and label study-specfic effects
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),
study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin")),
study = factor(study, levels = c(studies_long, "Overall"),
labels = c(studies_long, "Overall"))) %>%
# reshaping: remove extra columns, arrange by key variables, and make wide
select(type, Outcome, Trait, Moderator, Covariate, study, term, est) %>%
arrange(type, Outcome, Trait, Moderator, Covariate, study, term, est) %>%
pivot_wider(names_from = "Trait", values_from = "est")
ipd2b_reg_tab## # A tibble: 1,824 × 11
## type Outcome Moderator Covariate study term E A C N O
## <chr> <fct> <fct> <fct> <fct> <fct> <chr> <chr> <chr> <chr> <chr>
## 1 Bayesian Crystallized Ability None Unadjusted BASE Personality 0.05<br>[-0.04, 0.16] <NA> <NA> <str… <str…
## 2 Bayesian Crystallized Ability None Unadjusted EAS Personality <strong>0.10<br>[0.04, 0… 0.02… 0.02… <str… <str…
## 3 Bayesian Crystallized Ability None Unadjusted GSOEP Personality <strong>0.07<br>[0.03, 0… -0.0… -0.0… <str… <str…
## 4 Bayesian Crystallized Ability None Unadjusted HILDA Personality 0.01<br>[-0.01, 0.04] <str… <str… <str… <str…
## 5 Bayesian Crystallized Ability None Unadjusted HRS Personality 0.001<br>[-0.02, 0.03] <str… <str… <str… <str…
## 6 Bayesian Crystallized Ability None Unadjusted MAP Personality 0.04<br>[-0.01, 0.09] <NA> <str… <str… <NA>
## 7 Bayesian Crystallized Ability None Unadjusted OCTO-Twin Personality 0.05<br>[-0.02, 0.11] <NA> <NA> <str… <NA>
## 8 Bayesian Crystallized Ability None Unadjusted ROS Personality 0.02<br>[-0.05, 0.07] <str… <str… <str… <str…
## 9 Bayesian Crystallized Ability None Unadjusted SATSA Personality 0.04<br>[-0.02, 0.12] <str… <str… <str… <str…
## 10 Bayesian Crystallized Ability None Unadjusted Overall Personality <strong>0.04<br>[0.00, 0… 0.00… 0.04… <str… <str…
## # ℹ 1,814 more rowsipd2b_res <- nested_ipd2b_reg %>%
select(-rx) %>%
unnest(fx) %>% # unnesting
# keep key terms
filter((Moderator == "none" & term == "p_value") |
(Moderator != "none" & grepl("p_value:", term)) &
!grepl("cor", term) & !grepl("sd", term)) %>%
mutate(study = "Overall") %>%
mutate_at(vars(estimate, conf.low, conf.high),
~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
mutate(est = sprintf("\\textit{b} = %s, 95\\%% CI [%s, %s]", estimate, conf.low, conf.high)) Now that we’ve formatted the values, we can group by moderators and save results as separate tables. Even though additional information could be included given that we have one outcome, we’ll stick with this split because it will make it easier for those using this tutorial who multiple traits, outcomes, covariates, and moderators.
4.2.2.7.1.1 Fixed Effects
## table function
ipd2b_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(" ", paste0("<strong>", traits$short_name, "</strong>"))
cln <- if(length(unique(d$term)) == 1) c("Covariates", 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") "2B Pooled One Stage Models with Random Effects: Overall Effects of Personality-Crystallized Domain Associations" else sprintf("2B Pooled One Stage Models with Random Effects: Overall %s Moderation of Personality-Crystallized Domain Associations", 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, escape = F)
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/2b_ipd_mlm/%s/tables/overall/%s.html"
, local_path, type, md))
return(tab)
}
ipd2b_fx_tab <- ipd2b_reg_tab %>%
filter(study == "Overall") %>%
select(-study) %>%
group_by(type, Moderator) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Moderator), ipd2b_tab_fun))
# save(ipd2b_reg_tab, ipd2b_res, file = sprintf("%s/manuscript/results/ipd2b_fx_tab.RData", res_path))
## Frequentist
(ipd2b_fx_tab %>% filter(Moderator == "None" & type == "Frequentist"))$tab[[1]]|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Covariates | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| Crystallized Ability | |||||
| Unadjusted |
0.04 [0.01, 0.08] |
0.01 [-0.08, 0.09] |
0.05 [-0.03, 0.12] |
-0.12 [-0.15, -0.09] |
0.27 [0.18, 0.36] |
| Age |
0.03 [0.01, 0.06] |
0.02 [-0.08, 0.11] |
0.04 [-0.04, 0.12] |
-0.12 [-0.16, -0.09] |
0.28 [0.19, 0.38] |
| Gender |
0.03 [0.01, 0.06] |
0.02 [-0.07, 0.11] |
0.04 [-0.04, 0.12] |
-0.12 [-0.15, -0.09] |
0.28 [0.18, 0.38] |
| Education |
0.01 [-0.02, 0.05] |
0.01 [-0.02, 0.04] |
0.02 [-0.04, 0.07] |
-0.06 [-0.09, -0.03] |
0.16 [0.07, 0.25] |
| Fully Adjusted |
0.01 [-0.02, 0.05] |
0.00 [-0.03, 0.03] |
0.01 [-0.04, 0.07] |
-0.06 [-0.09, -0.03] |
0.16 [0.07, 0.25] |
|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Covariates | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| Crystallized Ability | |||||
| Unadjusted |
0.04 [0.00, 0.09] |
0.00 [-0.13, 0.12] |
0.04 [-0.07, 0.15] |
-0.12 [-0.16, -0.09] |
0.28 [0.16, 0.40] |
| Age |
0.03 [0.00, 0.07] |
0.01 [-0.11, 0.14] |
0.04 [-0.07, 0.16] |
-0.12 [-0.16, -0.08] |
0.29 [0.17, 0.42] |
| Gender |
0.03 [-0.001, 0.07] |
0.01 [-0.12, 0.14] |
0.04 [-0.07, 0.15] |
-0.12 [-0.17, -0.08] |
0.28 [0.15, 0.43] |
| Education |
0.01 [-0.03, 0.06] |
0.01 [-0.06, 0.08] |
0.02 [-0.06, 0.09] |
-0.06 [-0.10, -0.02] |
0.16 [0.02, 0.29] |
| Fully Adjusted |
0.01 [-0.03, 0.06] |
0.00 [-0.07, 0.07] |
0.02 [-0.07, 0.09] |
-0.06 [-0.10, -0.03] |
0.16 [0.04, 0.29] |
ipd2b_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$short_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("2B Pooled One Stage Models with Random Effects: Fixed Effect 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/2b_ipd_mlm/%s/tables/key terms/%s.html"
, local_path, type, covar))
return(tab) # return the html table
}
ipd2b_fx_tab2 <- ipd2b_reg_tab %>%
filter(study == "Overall") %>%
select(-study) %>%
arrange(Moderator, term) %>%
filter(Covariate %in% c("Unadjusted", "Fully Adjusted")) %>%
group_by(Outcome, type, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Covariate), ipd2b_tab_fun))
## Frequentist, no moderator
(ipd2b_fx_tab2 %>% filter(type == "Frequentist" & Covariate == "Fully Adjusted"))$tab[[1]]|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Term | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| None | |||||
| Personality |
0.01 [-0.02, 0.05] |
0.00 [-0.03, 0.03] |
0.01 [-0.04, 0.07] |
-0.06 [-0.09, -0.03] |
0.16 [0.07, 0.25] |
| Age | |||||
| Age |
0.001 [-0.03, 0.03] |
0.000 [-0.09, 0.09] |
-0.001 [-0.01, 0.01] |
0.00 [-0.00, 0.01] |
0.001 [-0.09, 0.09] |
| Gender | |||||
| Gender (Male v Female) |
0.05 [0.01, 0.09] |
0.02 [-0.02, 0.06] |
0.01 [-0.03, 0.06] |
-0.02 [-0.05, 0.01] |
-0.000 [-0.06, 0.06] |
| Education | |||||
| Education (Years) |
-0.01 [-0.07, 0.06] |
-0.01 [-0.02, 0.01] |
0.01 [-0.02, 0.04] |
0.00 [-0.01, 0.01] |
-0.01 [-0.12, 0.09] |
| Continent | |||||
| Continent (North America v Europe) |
0.07 [0.00, 0.14] |
-0.12 [-0.26, 0.03] |
-0.12 [-0.24, 0.00] |
-0.03 [-0.06, 0.01] |
0.14 [-0.08, 0.36] |
| Continent (North America v Australia) |
0.03 [-0.05, 0.10] |
-0.01 [-0.18, 0.15] |
-0.02 [-0.16, 0.11] |
-0.11 [-0.14, -0.08] |
0.11 [-0.19, 0.40] |
| Country | |||||
| Country (United States v Germany) |
0.08 [-0.02, 0.18] |
-0.01 [-0.10, 0.08] |
-0.05 [-0.19, 0.09] |
-0.02 [-0.09, 0.05] |
0.12 [-0.15, 0.38] |
| Country (United States v Sweden) |
0.07 [-0.03, 0.16] |
-0.23 [-0.35, -0.12] |
-0.20 [-0.34, -0.06] |
-0.05 [-0.11, 0.01] |
0.22 [-0.13, 0.58] |
| Country (United States v The Netherlands) |
-0.02 [-0.06, 0.03] |
||||
| Country (United States v Australia) |
0.02 [-0.06, 0.11] |
-0.01 [-0.05, 0.03] |
-0.02 [-0.13, 0.08] |
-0.11 [-0.14, -0.08] |
0.11 [-0.21, 0.43] |
| Personality Scale | |||||
| Scale (NEO-FFI v DPQ) |
0.01 [-0.09, 0.11] |
||||
| Scale (NEO-FFI v Eysenck) |
0.09 [-0.01, 0.18] |
-0.24 [-9.45, 8.97] |
-0.21 [-0.33, -0.10] |
-0.02 [-0.11, 0.08] |
0.14 [-0.50, 0.77] |
| Scale (NEO-FFI v MIDI) |
-0.00 [-0.08, 0.08] |
-0.01 [-9.22, 9.20] |
0.01 [-0.06, 0.08] |
0.03 [-0.06, 0.12] |
-0.07 [-0.70, 0.55] |
| Scale (NEO-FFI v BFI-S) |
0.12 [0.01, 0.23] |
-0.01 [-9.22, 9.20] |
-0.07 [-0.18, 0.05] |
0.03 [-0.09, 0.14] |
-0.12 [-0.75, 0.51] |
| Scale (NEO-FFI v IPIP NEO) |
0.11 [-0.01, 0.22] |
0.01 [-9.20, 9.22] |
-0.12 [-0.22, -0.01] |
-0.00 [-0.12, 0.12] |
0.01 [-0.61, 0.64] |
| Scale (NEO-FFI v TDA-40) |
0.04 [-0.04, 0.13] |
-0.01 [-9.22, 9.20] |
-0.04 [-0.11, 0.02] |
-0.08 [-0.17, 0.01] |
0.02 [-0.60, 0.65] |
| Baseline Age | |||||
| Study Baseline Age |
-0.001 [-0.00, 0.00] |
0.00 [-0.00, 0.01] |
0.00 [-0.00, 0.01] |
0.001 [-0.001, 0.00] |
-0.000 [-0.01, 0.01] |
| Baseline Year | |||||
| Study Baseline Year |
-0.000 [-0.01, 0.00] |
0.01 [0.00, 0.01] |
0.01 [-0.00, 0.01] |
0.001 [-0.00, 0.01] |
-0.01 [-0.02, 0.00] |
| Prediction Interval | |||||
| Prediction Interval |
-0.00 [-0.01, 0.01] |
-0.01 [-0.02, -0.00] |
-0.01 [-0.02, 0.01] |
-0.00 [-0.01, 0.01] |
0.01 [-0.02, 0.04] |
## bayesian
(ipd2b_fx_tab2 %>% filter(type == "Bayesian" & Covariate == "Fully Adjusted"))$tab[[1]]# save(ipd1a_fx_tab, ipd1a_fx_tab2, ipd1a_res, file = sprintf("%s/manuscript/results/ipd1b_fx_tab.RData", res_path))|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Term | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| None | |||||
| Personality |
0.01 [-0.03, 0.06] |
0.00 [-0.07, 0.07] |
0.02 [-0.07, 0.09] |
-0.06 [-0.10, -0.03] |
0.16 [0.04, 0.29] |
| Age | |||||
| Age |
-0.00 [-0.01, 0.00] |
-0.001 [-0.01, 0.00] |
-0.00 [-0.01, 0.00] |
0.001 [-0.00, 0.01] |
-0.001 [-0.01, 0.01] |
| Gender | |||||
| Gender (Male v Female) |
0.04 [-0.01, 0.10] |
0.00 [-0.08, 0.07] |
-0.02 [-0.09, 0.04] |
-0.01 [-0.06, 0.03] |
-0.04 [-0.09, 0.04] |
| Education | |||||
| Education (Years) |
-0.00 [-0.01, 0.01] |
-0.00 [-0.02, 0.02] |
0.01 [-0.02, 0.04] |
0.000 [-0.01, 0.01] |
-0.00 [-0.02, 0.01] |
| Continent | |||||
| Continent (North America v Europe) |
0.07 [-0.04, 0.16] |
-0.08 [-0.35, 0.15] |
-0.11 [-0.30, 0.12] |
-0.03 [-0.07, 0.02] |
0.14 [-0.09, 0.40] |
| Continent (North America v Australia) |
0.02 [-0.12, 0.14] |
-0.03 [-0.31, 0.27] |
-0.03 [-0.27, 0.22] |
-0.11 [-0.17, -0.04] |
0.09 [-0.24, 0.41] |
| Country | |||||
| Country (United States v Germany) |
0.08 [-0.06, 0.21] |
0.01 [-0.15, 0.21] |
-0.04 [-0.29, 0.27] |
-0.02 [-0.13, 0.08] |
0.11 [-0.14, 0.37] |
| Country (United States v Sweden) |
0.07 [-0.07, 0.19] |
-0.22 [-0.40, -0.03] |
-0.18 [-0.42, 0.11] |
-0.04 [-0.13, 0.05] |
0.24 [-0.09, 0.58] |
| Country (United States v The Netherlands) |
-0.01 [-0.11, 0.09] |
||||
| Country (United States v Australia) |
0.02 [-0.14, 0.16] |
-0.01 [-0.15, 0.13] |
-0.03 [-0.29, 0.20] |
-0.11 [-0.19, -0.01] |
0.08 [-0.23, 0.39] |
| Personality Scale | |||||
| Scale (NEO-FFI v DPQ) |
0.01 [-0.17, 0.22] |
||||
| Scale (NEO-FFI v Eysenck) |
0.08 [-0.08, 0.28] |
-0.21 [-1.49, 1.03] |
-0.21 [-0.58, 0.14] |
-0.02 [-0.16, 0.13] |
0.19 [-0.68, 1.10] |
| Scale (NEO-FFI v MIDI) |
-0.01 [-0.24, 0.21] |
-0.12 [-1.40, 1.05] |
0.00 [-0.40, 0.30] |
0.04 [-0.14, 0.25] |
-0.07 [-0.90, 0.84] |
| Scale (NEO-FFI v BFI-S) |
0.12 [-0.09, 0.32] |
-0.04 [-1.55, 1.11] |
-0.06 [-0.41, 0.28] |
0.03 [-0.14, 0.24] |
-0.01 [-0.89, 0.97] |
| Scale (NEO-FFI v IPIP NEO) |
0.10 [-0.11, 0.32] |
0.03 [-1.11, 1.27] |
-0.11 [-0.48, 0.23] |
-0.00 [-0.20, 0.20] |
0.05 [-0.84, 1.04] |
| Scale (NEO-FFI v TDA-40) |
0.04 [-0.19, 0.23] |
-0.12 [-1.52, 0.95] |
-0.05 [-0.43, 0.38] |
-0.08 [-0.23, 0.12] |
0.000 [-0.91, 0.83] |
| Baseline Age | |||||
| Study Baseline Age |
-0.001 [-0.00, 0.00] |
0.00 [-0.01, 0.01] |
0.00 [-0.00, 0.01] |
0.001 [-0.00, 0.00] |
-0.000 [-0.01, 0.01] |
| Baseline Year | |||||
| Study Baseline Year |
-0.000 [-0.01, 0.01] |
0.01 [0.00, 0.01] |
0.01 [-0.01, 0.02] |
0.001 [-0.00, 0.01] |
-0.01 [-0.02, 0.01] |
| Prediction Interval | |||||
| Prediction Interval |
-0.00 [-0.02, 0.01] |
-0.01 [-0.03, 0.000] |
-0.01 [-0.03, 0.01] |
-0.001 [-0.01, 0.01] |
0.01 [-0.02, 0.05] |
4.2.2.7.1.2 Study-Specific Effects
## table function
ipd2b_tab_fun <- function(d, type, moder){
print(paste(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) c(2, rep(1,5)) else c(3, rep(1,5))
names(cs) <- c(" ", paste0("<strong>", traits$short_name, "</strong>"))
cln <- if(length(unique(d$term)) == 1) c("Covariate", "Study", rep("<em>b</em> [CI]", 5)) else c("Covariate", "Study", "Term", rep("<em>b</em> [CI]", 5))
al <- if(length(unique(d$term)) == 1) c("r", "r", rep("c", 5)) else c("r", "r", "r", rep("c", 5))
if(length(unique(d$term)) == 1) d <- d %>% select(-term)
cap <- if(md == "none") "Overall Effects of Personality-Crystallized Domain Associations" else sprintf("Overall %s Moderation of Personality-Crystallized Domain Associations", md)
tab <- d %>%
arrange(Outcome) %>%
select(-Outcome) %>%
kable(., "html"
, escape = F
, col.names = cln
, align = al
, caption = cap
) %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
add_header_above(cs, escape = F) %>%
collapse_rows(1, valign = "top")
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/2b_ipd_mlm/%s/tables/study specific/%s.html"
, local_path, type, md))
return(tab)
}
ipd2b_rx_tab <- ipd2b_reg_tab %>%
filter(!Moderator %in% stdyModers$long_name) %>%
arrange(type, Outcome, Moderator, Covariate, study, term) %>%
group_by(type, Moderator) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Moderator), ipd2b_tab_fun))## [1] "Bayesian None"
## [1] "Bayesian Age"
## [1] "Bayesian Gender"
## [1] "Bayesian Education"
## [1] "Frequentist None"
## [1] "Frequentist Age"
## [1] "Frequentist Gender"
## [1] "Frequentist Education"## # A tibble: 8 × 4
## type Moderator data tab
## <chr> <fct> <list> <list>
## 1 Bayesian None <tibble [60 × 9]> <kablExtr [1]>
## 2 Bayesian Age <tibble [24 × 9]> <kablExtr [1]>
## 3 Bayesian Gender <tibble [24 × 9]> <kablExtr [1]>
## 4 Bayesian Education <tibble [24 × 9]> <kablExtr [1]>
## 5 Frequentist None <tibble [60 × 9]> <kablExtr [1]>
## 6 Frequentist Age <tibble [24 × 9]> <kablExtr [1]>
## 7 Frequentist Gender <tibble [24 × 9]> <kablExtr [1]>
## 8 Frequentist Education <tibble [24 × 9]> <kablExtr [1]>|
E
|
A
|
C
|
N
|
O
|
||
|---|---|---|---|---|---|---|
| Covariate | Study | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| Crystallized Ability | ||||||
| Unadjusted | BASE |
0.05 [-0.02, 0.12] |
-0.13 [-0.20, -0.06] |
0.39 [0.27, 0.52] |
||
| EAS |
0.10 [0.06, 0.15] |
0.02 [-0.03, 0.07] |
0.02 [-0.03, 0.07] |
-0.13 [-0.17, -0.08] |
0.37 [0.32, 0.42] |
|
| GSOEP |
0.07 [0.03, 0.11] |
-0.02 [-0.07, 0.03] |
-0.02 [-0.07, 0.03] |
-0.08 [-0.12, -0.04] |
0.10 [0.06, 0.15] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.05 [0.01, 0.08] |
0.07 [0.04, 0.09] |
-0.15 [-0.17, -0.12] |
0.27 [0.25, 0.30] |
|
| HRS |
0.001 [-0.02, 0.03] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.08 [-0.10, -0.06] |
0.22 [0.19, 0.24] |
|
| LASA |
-0.13 [-0.17, -0.09] |
|||||
| MAP |
0.04 [-0.01, 0.09] |
0.09 [0.02, 0.17] |
-0.15 [-0.20, -0.10] |
|||
| MARS |
-0.18 [-0.25, -0.12] |
|||||
| OCTO-Twin |
0.05 [-0.01, 0.11] |
-0.10 [-0.16, -0.04] |
||||
| ROS |
0.02 [-0.03, 0.07] |
0.10 [0.02, 0.17] |
0.12 [0.04, 0.19] |
-0.10 [-0.15, -0.05] |
0.19 [0.11, 0.27] |
|
| SATSA |
0.04 [-0.02, 0.11] |
-0.17 [-0.27, -0.07] |
-0.12 [-0.21, -0.03] |
-0.11 [-0.17, -0.05] |
0.37 [0.25, 0.49] |
|
| Overall |
0.04 [0.01, 0.08] |
0.01 [-0.08, 0.09] |
0.05 [-0.03, 0.12] |
-0.12 [-0.15, -0.09] |
0.27 [0.18, 0.36] |
|
| Age | BASE |
0.03 [-0.02, 0.07] |
-0.14 [-0.21, -0.06] |
0.40 [0.28, 0.53] |
||
| EAS |
0.05 [0.01, 0.09] |
0.08 [0.01, 0.16] |
0.01 [-0.07, 0.10] |
-0.13 [-0.19, -0.07] |
0.41 [0.33, 0.49] |
|
| GSOEP |
0.06 [0.03, 0.09] |
-0.02 [-0.07, 0.03] |
-0.02 [-0.08, 0.03] |
-0.08 [-0.12, -0.04] |
0.10 [0.06, 0.15] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.04 [0.01, 0.07] |
0.07 [0.04, 0.09] |
-0.16 [-0.19, -0.14] |
0.28 [0.25, 0.30] |
|
| HRS |
0.00 [-0.02, 0.02] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.08 [-0.10, -0.06] |
0.22 [0.19, 0.24] |
|
| LASA |
-0.12 [-0.16, -0.09] |
|||||
| MAP |
0.04 [0.01, 0.08] |
0.09 [0.02, 0.17] |
-0.15 [-0.20, -0.10] |
|||
| MARS |
-0.19 [-0.26, -0.12] |
|||||
| OCTO-Twin |
0.03 [-0.01, 0.07] |
-0.10 [-0.16, -0.04] |
||||
| ROS |
0.03 [-0.01, 0.06] |
0.10 [0.03, 0.18] |
0.12 [0.05, 0.19] |
-0.10 [-0.16, -0.05] |
0.19 [0.12, 0.27] |
|
| SATSA |
0.04 [-0.00, 0.08] |
-0.18 [-0.28, -0.08] |
-0.13 [-0.22, -0.04] |
-0.11 [-0.17, -0.05] |
0.38 [0.26, 0.50] |
|
| Overall |
0.03 [0.01, 0.06] |
0.02 [-0.08, 0.11] |
0.04 [-0.04, 0.12] |
-0.12 [-0.16, -0.09] |
0.28 [0.19, 0.38] |
|
| Gender | BASE |
0.03 [-0.02, 0.07] |
-0.13 [-0.21, -0.06] |
0.40 [0.27, 0.53] |
||
| EAS |
0.05 [0.01, 0.09] |
0.09 [0.01, 0.16] |
0.01 [-0.07, 0.10] |
-0.13 [-0.19, -0.07] |
0.41 [0.33, 0.49] |
|
| GSOEP |
0.06 [0.03, 0.09] |
-0.02 [-0.07, 0.03] |
-0.02 [-0.08, 0.03] |
-0.08 [-0.12, -0.04] |
0.10 [0.06, 0.14] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.05 [0.02, 0.08] |
0.07 [0.04, 0.09] |
-0.15 [-0.17, -0.12] |
0.27 [0.25, 0.30] |
|
| HRS |
0.00 [-0.02, 0.02] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.08 [-0.10, -0.06] |
0.22 [0.19, 0.24] |
|
| LASA |
-0.13 [-0.17, -0.09] |
|||||
| MAP |
0.04 [0.01, 0.08] |
0.09 [0.02, 0.17] |
-0.15 [-0.20, -0.10] |
|||
| MARS |
-0.19 [-0.25, -0.12] |
|||||
| OCTO-Twin |
0.03 [-0.01, 0.07] |
-0.10 [-0.16, -0.04] |
||||
| ROS |
0.03 [-0.01, 0.06] |
0.10 [0.03, 0.18] |
0.12 [0.05, 0.19] |
-0.10 [-0.16, -0.05] |
0.19 [0.11, 0.27] |
|
| SATSA |
0.04 [-0.00, 0.08] |
-0.18 [-0.28, -0.08] |
-0.13 [-0.22, -0.04] |
-0.11 [-0.17, -0.05] |
0.38 [0.26, 0.50] |
|
| Overall |
0.03 [0.01, 0.06] |
0.02 [-0.07, 0.11] |
0.04 [-0.04, 0.12] |
-0.12 [-0.15, -0.09] |
0.28 [0.18, 0.38] |
|
| Education | BASE |
0.00 [-0.06, 0.07] |
-0.09 [-0.15, -0.02] |
0.30 [0.18, 0.42] |
||
| EAS |
0.03 [-0.02, 0.09] |
0.03 [0.02, 0.03] |
-0.02 [-0.09, 0.05] |
-0.07 [-0.12, -0.01] |
0.19 [0.11, 0.27] |
|
| GSOEP |
0.06 [0.01, 0.12] |
-0.04 [-0.04, -0.03] |
0.00 [-0.07, 0.08] |
-0.03 [-0.08, 0.02] |
0.06 [-0.01, 0.13] |
|
| HILDA |
0.01 [-0.01, 0.03] |
0.04 [0.04, 0.04] |
0.03 [0.01, 0.06] |
-0.13 [-0.15, -0.11] |
0.19 [0.17, 0.22] |
|
| HRS |
-0.03 [-0.05, -0.01] |
0.04 [0.04, 0.04] |
0.08 [0.06, 0.11] |
-0.03 [-0.05, -0.01] |
0.10 [0.08, 0.13] |
|
| LASA |
-0.05 [-0.09, -0.01] |
|||||
| MAP |
0.02 [-0.03, 0.07] |
0.04 [-0.03, 0.10] |
-0.04 [-0.08, 0.01] |
|||
| MARS |
-0.08 [-0.14, -0.02] |
|||||
| OCTO-Twin |
0.04 [-0.02, 0.09] |
-0.07 [-0.13, -0.02] |
||||
| ROS |
-0.04 [-0.09, 0.01] |
0.01 [0.01, 0.02] |
0.06 [-0.00, 0.12] |
-0.03 [-0.08, 0.02] |
0.01 [-0.07, 0.08] |
|
| SATSA |
0.03 [-0.03, 0.09] |
-0.02 [-0.03, -0.02] |
-0.08 [-0.16, -0.01] |
-0.05 [-0.11, 0.00] |
0.26 [0.14, 0.38] |
|
| Overall |
0.01 [-0.02, 0.05] |
0.01 [-0.02, 0.04] |
0.02 [-0.04, 0.07] |
-0.06 [-0.09, -0.03] |
0.16 [0.07, 0.25] |
|
| Fully Adjusted | BASE |
0.00 [-0.06, 0.07] |
-0.09 [-0.16, -0.02] |
0.31 [0.18, 0.43] |
||
| EAS |
0.03 [-0.02, 0.09] |
0.02 [0.01, 0.02] |
-0.02 [-0.09, 0.05] |
-0.07 [-0.12, -0.01] |
0.19 [0.11, 0.27] |
|
| GSOEP |
0.06 [0.01, 0.12] |
-0.04 [-0.04, -0.04] |
0.000 [-0.08, 0.08] |
-0.03 [-0.08, 0.02] |
0.06 [-0.00, 0.13] |
|
| HILDA |
0.01 [-0.02, 0.03] |
0.03 [0.03, 0.03] |
0.03 [0.00, 0.05] |
-0.13 [-0.16, -0.11] |
0.20 [0.17, 0.22] |
|
| HRS |
-0.03 [-0.05, -0.01] |
0.03 [0.03, 0.03] |
0.08 [0.06, 0.11] |
-0.03 [-0.05, -0.01] |
0.10 [0.08, 0.13] |
|
| LASA |
-0.05 [-0.09, -0.01] |
|||||
| MAP |
0.02 [-0.02, 0.07] |
0.04 [-0.03, 0.11] |
-0.03 [-0.08, 0.01] |
|||
| MARS |
-0.08 [-0.14, -0.02] |
|||||
| OCTO-Twin |
0.04 [-0.02, 0.09] |
-0.07 [-0.13, -0.02] |
||||
| ROS |
-0.04 [-0.09, 0.01] |
0.01 [0.00, 0.01] |
0.06 [-0.000, 0.13] |
-0.03 [-0.08, 0.02] |
0.01 [-0.06, 0.09] |
|
| SATSA |
0.03 [-0.03, 0.09] |
-0.03 [-0.03, -0.02] |
-0.09 [-0.16, -0.01] |
-0.05 [-0.11, 0.00] |
0.26 [0.14, 0.38] |
|
| Overall |
0.01 [-0.02, 0.05] |
0.00 [-0.03, 0.03] |
0.01 [-0.04, 0.07] |
-0.06 [-0.09, -0.03] |
0.16 [0.07, 0.25] |
|
|
E
|
A
|
C
|
N
|
O
|
||
|---|---|---|---|---|---|---|
| Covariate | Study | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| Crystallized Ability | ||||||
| Unadjusted | BASE |
0.05 [-0.04, 0.16] |
-0.14 [-0.23, -0.06] |
0.39 [0.26, 0.53] |
||
| EAS |
0.10 [0.04, 0.16] |
0.02 [-0.03, 0.07] |
0.02 [-0.03, 0.07] |
-0.13 [-0.17, -0.08] |
0.37 [0.32, 0.42] |
|
| GSOEP |
0.07 [0.03, 0.12] |
-0.02 [-0.07, 0.04] |
-0.02 [-0.08, 0.03] |
-0.08 [-0.12, -0.03] |
0.10 [0.06, 0.15] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.04 [0.01, 0.08] |
0.07 [0.04, 0.09] |
-0.15 [-0.17, -0.12] |
0.27 [0.25, 0.30] |
|
| HRS |
0.001 [-0.02, 0.03] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.08 [-0.10, -0.06] |
0.22 [0.19, 0.24] |
|
| LASA |
-0.13 [-0.17, -0.09] |
|||||
| MAP |
0.04 [-0.01, 0.09] |
0.10 [0.02, 0.17] |
-0.15 [-0.20, -0.09] |
|||
| MARS |
-0.19 [-0.29, -0.11] |
|||||
| OCTO-Twin |
0.05 [-0.02, 0.11] |
-0.10 [-0.17, -0.03] |
||||
| ROS |
0.02 [-0.05, 0.07] |
0.10 [0.02, 0.18] |
0.12 [0.05, 0.20] |
-0.10 [-0.15, -0.04] |
0.19 [0.12, 0.27] |
|
| SATSA |
0.04 [-0.02, 0.12] |
-0.17 [-0.30, -0.04] |
-0.13 [-0.22, -0.03] |
-0.11 [-0.17, -0.04] |
0.38 [0.26, 0.52] |
|
| Overall |
0.04 [0.00, 0.09] |
0.00 [-0.13, 0.12] |
0.04 [-0.07, 0.15] |
-0.12 [-0.16, -0.09] |
0.28 [0.16, 0.40] |
|
| Age | BASE |
0.03 [-0.03, 0.12] |
-0.14 [-0.24, -0.06] |
0.39 [0.26, 0.53] |
||
| EAS |
0.05 [0.00, 0.14] |
0.08 [0.00, 0.17] |
0.01 [-0.07, 0.09] |
-0.13 [-0.19, -0.06] |
0.41 [0.32, 0.49] |
|
| GSOEP |
0.06 [0.02, 0.11] |
-0.02 [-0.07, 0.04] |
-0.02 [-0.07, 0.03] |
-0.08 [-0.12, -0.03] |
0.10 [0.06, 0.15] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.04 [0.01, 0.07] |
0.07 [0.04, 0.09] |
-0.16 [-0.19, -0.14] |
0.28 [0.25, 0.30] |
|
| HRS |
0.00 [-0.02, 0.03] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.08 [-0.10, -0.06] |
0.22 [0.19, 0.24] |
|
| LASA |
-0.12 [-0.16, -0.08] |
|||||
| MAP |
0.04 [-0.01, 0.09] |
0.10 [0.02, 0.17] |
-0.15 [-0.20, -0.09] |
|||
| MARS |
-0.19 [-0.29, -0.11] |
|||||
| OCTO-Twin |
0.03 [-0.02, 0.10] |
-0.10 [-0.17, -0.03] |
||||
| ROS |
0.02 [-0.04, 0.07] |
0.10 [0.03, 0.18] |
0.12 [0.05, 0.20] |
-0.10 [-0.15, -0.04] |
0.20 [0.12, 0.27] |
|
| SATSA |
0.04 [-0.01, 0.10] |
-0.18 [-0.30, -0.05] |
-0.13 [-0.23, -0.03] |
-0.11 [-0.17, -0.04] |
0.40 [0.26, 0.54] |
|
| Overall |
0.03 [0.00, 0.07] |
0.01 [-0.11, 0.14] |
0.04 [-0.07, 0.16] |
-0.12 [-0.16, -0.08] |
0.29 [0.17, 0.42] |
|
| Gender | BASE |
0.03 [-0.03, 0.12] |
-0.14 [-0.24, -0.06] |
0.39 [0.26, 0.53] |
||
| EAS |
0.05 [0.01, 0.14] |
0.08 [0.01, 0.16] |
0.01 [-0.07, 0.09] |
-0.13 [-0.19, -0.06] |
0.41 [0.32, 0.49] |
|
| GSOEP |
0.06 [0.02, 0.11] |
-0.02 [-0.07, 0.03] |
-0.02 [-0.07, 0.03] |
-0.08 [-0.12, -0.03] |
0.10 [0.06, 0.15] |
|
| HILDA |
0.01 [-0.01, 0.04] |
0.05 [0.02, 0.08] |
0.07 [0.04, 0.09] |
-0.15 [-0.17, -0.12] |
0.27 [0.25, 0.30] |
|
| HRS |
0.00 [-0.02, 0.03] |
0.06 [0.03, 0.09] |
0.16 [0.13, 0.19] |
-0.08 [-0.10, -0.06] |
0.22 [0.19, 0.24] |
|
| LASA |
-0.13 [-0.17, -0.09] |
|||||
| MAP |
0.04 [-0.01, 0.09] |
0.10 [0.02, 0.17] |
-0.15 [-0.20, -0.09] |
|||
| MARS |
-0.19 [-0.29, -0.11] |
|||||
| OCTO-Twin |
0.03 [-0.02, 0.10] |
-0.10 [-0.17, -0.03] |
||||
| ROS |
0.02 [-0.03, 0.07] |
0.10 [0.03, 0.18] |
0.12 [0.05, 0.20] |
-0.10 [-0.16, -0.04] |
0.19 [0.12, 0.27] |
|
| SATSA |
0.04 [-0.01, 0.10] |
-0.18 [-0.29, -0.06] |
-0.13 [-0.23, -0.03] |
-0.11 [-0.17, -0.04] |
0.39 [0.26, 0.53] |
|
| Overall |
0.03 [-0.001, 0.07] |
0.01 [-0.12, 0.14] |
0.04 [-0.07, 0.15] |
-0.12 [-0.17, -0.08] |
0.28 [0.15, 0.43] |
|
| Education | BASE |
0.01 [-0.08, 0.10] |
-0.09 [-0.17, -0.02] |
0.30 [0.17, 0.44] |
||
| EAS |
0.04 [-0.03, 0.11] |
0.04 [-0.02, 0.10] |
-0.02 [-0.10, 0.05] |
-0.06 [-0.12, -0.00] |
0.18 [0.11, 0.26] |
|
| GSOEP |
0.06 [-0.00, 0.14] |
0.01 [-0.06, 0.09] |
0.01 [-0.07, 0.09] |
-0.03 [-0.09, 0.03] |
0.07 [-0.00, 0.14] |
|
| HILDA |
0.01 [-0.01, 0.03] |
0.04 [0.01, 0.06] |
0.03 [0.01, 0.06] |
-0.13 [-0.15, -0.11] |
0.19 [0.17, 0.22] |
|
| HRS |
-0.03 [-0.05, -0.01] |
0.04 [0.01, 0.06] |
0.08 [0.06, 0.11] |
-0.03 [-0.05, -0.00] |
0.10 [0.08, 0.12] |
|
| LASA |
-0.05 [-0.09, -0.01] |
|||||
| MAP |
0.02 [-0.03, 0.07] |
0.05 [-0.02, 0.12] |
-0.04 [-0.09, 0.02] |
|||
| MARS |
-0.09 [-0.17, -0.02] |
|||||
| OCTO-Twin |
0.04 [-0.02, 0.11] |
-0.07 [-0.13, -0.01] |
||||
| ROS |
-0.04 [-0.10, 0.01] |
0.02 [-0.03, 0.08] |
0.06 [-0.00, 0.13] |
-0.03 [-0.08, 0.03] |
0.01 [-0.07, 0.08] |
|
| SATSA |
0.03 [-0.04, 0.09] |
-0.08 [-0.22, 0.03] |
-0.08 [-0.19, 0.01] |
-0.05 [-0.12, 0.01] |
0.27 [0.14, 0.41] |
|
| Overall |
0.01 [-0.03, 0.06] |
0.01 [-0.06, 0.08] |
0.02 [-0.06, 0.09] |
-0.06 [-0.10, -0.02] |
0.16 [0.02, 0.29] |
|
| Fully Adjusted | BASE |
0.00 [-0.08, 0.09] |
-0.09 [-0.18, -0.01] |
0.30 [0.17, 0.44] |
||
| EAS |
0.04 [-0.03, 0.11] |
0.03 [-0.03, 0.10] |
-0.02 [-0.10, 0.05] |
-0.07 [-0.13, -0.01] |
0.19 [0.11, 0.26] |
|
| GSOEP |
0.06 [-0.00, 0.13] |
0.01 [-0.07, 0.09] |
0.01 [-0.07, 0.09] |
-0.03 [-0.09, 0.03] |
0.07 [-0.00, 0.14] |
|
| HILDA |
0.01 [-0.02, 0.03] |
0.03 [-0.001, 0.05] |
0.03 [0.00, 0.05] |
-0.13 [-0.16, -0.11] |
0.20 [0.17, 0.22] |
|
| HRS |
-0.03 [-0.06, -0.01] |
0.03 [0.00, 0.05] |
0.08 [0.05, 0.11] |
-0.03 [-0.05, -0.01] |
0.10 [0.08, 0.13] |
|
| LASA |
-0.05 [-0.09, -0.01] |
|||||
| MAP |
0.02 [-0.03, 0.07] |
0.05 [-0.02, 0.12] |
-0.03 [-0.09, 0.02] |
|||
| MARS |
-0.09 [-0.16, -0.02] |
|||||
| OCTO-Twin |
0.04 [-0.02, 0.11] |
-0.07 [-0.14, -0.02] |
||||
| ROS |
-0.04 [-0.10, 0.01] |
0.02 [-0.03, 0.09] |
0.06 [-0.001, 0.13] |
-0.03 [-0.08, 0.03] |
0.01 [-0.06, 0.09] |
|
| SATSA |
0.03 [-0.04, 0.10] |
-0.10 [-0.23, 0.02] |
-0.09 [-0.19, 0.01] |
-0.06 [-0.12, 0.01] |
0.27 [0.14, 0.41] |
|
| Overall |
0.01 [-0.03, 0.06] |
0.00 [-0.07, 0.07] |
0.02 [-0.07, 0.09] |
-0.06 [-0.10, -0.03] |
0.16 [0.04, 0.29] |
|
4.2.2.7.1.3 Heterogeneity Estimates
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/2b_ipd_mlm/%s/%s/%s", local_path, type, folder, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
## load in "fixed" effects
## first get file names
nested_ipd2b_het <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/2b_ipd_mlm/%s/heterogeneity", 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, ~loadRData(.x, .y, "het", "heterogeneity"))) %>%
select(-file) ip2b_hetero_tab_fun <- function(d, type, out, mod, cov){
moder <- mapvalues(mod, moders$short_name, moders$long_name, warn_missing = F)
d2 <- d %>%
separate(term, c("V1", "V2"), sep = "[.]") %>%
mutate(V2 = ifelse(is.na(V2), V1, V2)
, est = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high)) %>%
mutate_at(vars(V1, V2), ~mapvalues(., moders$short_term, moders$long_term)) %>%
mutate_at(vars(V1, V2), ~str_replace_all(., "p_value", "Personality")) %>%
mutate_at(vars(V1, V2), ~str_replace_all(., ":", " : ")) %>%
mutate_at(vars(V1, V2), str_to_title) %>%
mutate_at(vars(V1, V2), ~str_replace_all(., ":", "x")) %>%
mutate_at(vars(V1, V2), ~str_remove_all(., "[()]")) %>%
mutate_at(vars(V1, V2), ~ifelse(. == "Observation", "Sigma", .)) %>%
select(-group, -estimate, -conf.low, -conf.high) %>%
# unite(V1, Trait, V1, sep = "_", remove = T) %>%
pivot_wider(names_from = "V1", values_from = "est") %>%
rename(" " = V2) %>%
mutate(Trait = factor(Trait, traits$short_name, traits$long_name))
rs <- d2 %>% group_by(Trait) %>% tally() %>%
mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
cap <- if(mod == "none") "Heterogeneity, Variance, and Correlations Among Random Effects for Overall Effects of Personality-Crystallized Domain Associations" else sprintf("Heterogeneity, Variance, and Correlations Among Random Effects for Overall %s Moderation of Personality-Crystallized Domain Associations", moder)
cap <- sprintf("<strong>Table SX</strong><br><em>%s</em>", cap)
tab <- d2 %>%
select(-Trait) %>%
kable(.
, "html"
, align = c("r", rep("c", ncol(d2)-2))
, caption = cap
) %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
footnote("Diagnonal elements represent variances, while off diagnal elements represent correlations. Interval estimates are 95% bootstrapped CI's.")
for (i in 1:nrow(rs)) {
tab <- tab %>%
kableExtra::group_rows(rs$Trait[i], rs$start[i], rs$end[i])
}
save_kable(tab, file = sprintf("%s/results/2b_ipd_mlm/%s/tables/heterogeneity/%s-%s-%s.html", local_path, type, out, mod, cov))
return(tab)
}
nested_ipd2b_het_tab <- nested_ipd2b_het %>%
filter(!Moderator %in% stdyModers$short_name) %>%
unnest(het) %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Outcome, Moderator, Covariate), ip2b_hetero_tab_fun))4.2.2.7.1.4 All Model Terms
ipd2b_mod_tab <- nested_ipd2b_reg %>%
select(-n, -rx) %>%
unnest(fx) %>%
# 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))) %>%
# 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) %>%
arrange(type, Outcome, Trait, Moderator, Covariate) %>%
pivot_wider(names_from = "Trait", values_from = "est")
ipd2b_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(" ", paste0("<strong>", traits$long_name, "</strong>"))
# 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") "1B Pooled Analysis of Individual Participant Data with Cluster Corrected Standard Errors: All Model Estimates of Fixed Effect Personality-Crystallized Domain Associations" else sprintf("1B Pooled Analysis of Individual Participant Data with Cluster Corrected Standard Errors: All Model Estimates of Fixed Effect %s Moderation of Personality-Crystallized Domain Associations", md)
# kable the table
tab <- d %>%
arrange(term) %>%
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)
# save the resulting html table
save_kable(tab, file = sprintf("%s/results/2b_ipd_mlm/%s/tables/all terms/%s-%s-%s.html"
, local_path, type, o, md, cv))
return(tab) # return the html table
}
ipd2b_mod_tab <- ipd2b_mod_tab %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
mutate(tab = pmap(list(data, type, Outcome, Moderator, Covariate), ipd2b_mod_tab_fun))4.2.2.7.2 Figures
4.2.2.7.2.1 Overall Forest
ipd2b_fx_plot_fun <- function(df, mod, type, cov){
m <- mapvalues(mod, moders$short_name, moders$long_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", m)}
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"}
# if(length(unique(df$term)) > 1) df <- df %>% full_join(crossing(Trait = unique(df$Trait), Outcome = unique(df$Outcome), term = unique(df$term)))
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 = term, y = estimate)) +
scale_y_continuous(limits = lim, breaks = brk, labels = lab) +
scale_size_manual(values = c(1.8, 1.3)) +
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)
, width = .01
, 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 2B: Pooled Regression Using Random Effects"
) +
guides(color = "none", size = "none") +
facet_grid(Outcome~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(p, file = sprintf("%s/results/2b_ipd_mlm/%s/figures/overall forest/%s_%s_fixed.png", local_path, type, mod, cov), width = local_path*2, height = 1.25*ht + .75*ht2)
rm(p)
gc()
return(T)
}
nested_ipd2b_reg %>%
select(-rx, -n) %>%
unnest(fx) %>%
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)),
Outcome = forcats::fct_rev(Outcome)) %>%
group_by(type, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
# filter(Moderator == "country") %>%
mutate(pmap(list(data, Moderator, type, Covariate), ipd2b_fx_plot_fun))4.2.2.7.2.2 Study-Specific Forest
ipd2b_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"))) %>%
full_join(tibble(study = " ", estimate = NA, n = NA)) %>%
arrange(estimate)
stds <- df$study[!df$study %in% c("Overall", " ")]
df <- df %>%
mutate(study = factor(study, rev(c(" ", stds, "Overall")))
# , 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 == "Overall", "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(data = d2 %>% filter(est == "est"), aes(label = value), hjust = .5, size = 3.5) +
geom_text(data = d2 %>% filter(est == "n"), 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)
, 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)
gc()
save(p
, file = sprintf("%s/results/2b_ipd_mlm/%s/figures/study specific forest/rdata/%s_%s_%s_%s.RData", local_path, type, outcome, trait, mod, cov))
return(p)
}
## fixed effects
nested_ipd2b_reg_fp <- nested_ipd2b_reg %>%
filter(Moderator %in% moders$short_name) %>%
select(-rx, -n) %>%
unnest(fx) %>%
mutate(study = "Overall") %>%
## random effects
full_join(
nested_ipd2b_reg %>%
filter(Moderator %in% moders$short_name) %>%
mutate(rx = map2(rx, n, ~(.x) %>% full_join(.y))) %>%
select(-fx, -n) %>%
unnest(rx) #%>%
# mutate(term = ifelse(Moderator != "none", paste(term, mapvalues(Moderator, moders$short_name, moders$short_term, warn_missing = F), sep = ":"), term))
) %>%
## 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() %>%
# filter(Trait == "N" & type == "Frequentist") %>%
mutate(p = pmap(list(data, Outcome, Moderator, type, Covariate, Trait), ipd2b_rx_plot_fun))
ipd2b_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 2B: Pooled Regression Using Random Effects",
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/2b_ipd_mlm/%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/2b_ipd_mlm/%s/figures/study specific forest/%s_%s_%s.pdf", local_path, type, outcome, mod, cov)
, width = 10, height = 10)
return(T)
}
nested_ipd2b_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), ipd2b_rx_plot_comb_fun))4.2.2.7.2.3 Overall Simple Effects
loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
path <- sprintf("%s/results/2b_ipd_mlm/%s/%s/%s", local_path, type, folder, fileName)
load(path)
get(ls()[grepl(obj, ls())])
}
## load in "fixed" effects
## first get file names
nested_ipd2b_simp <- tibble(type = c("Frequentist", "Bayesian")) %>%
mutate(file = map(type, ~list.files(sprintf("%s/results/2b_ipd_mlm/%s/predicted", local_path, .)))) %>%
unnest(file) %>%
separate(file, c("Outcome", "Trait", "Moderator", "Covariate"), sep = "_", remove = F) %>%
## read in the files
mutate(Covariate = str_remove(Covariate, ".RData"),
pred.fx = map2(file, type, ~loadRData(.x, .y, "pred.fx", "predicted")),
pred.rx = map2(file, type, ~loadRData(.x, .y, "pred.rx", "predicted"))) %>%
select(-file)
nested_ipd2b_simpsimp_eff_fun <- function(df, outcome, mod, type, cov){
print(paste(outcome, mod))
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, c(moders$short_name, stdyModers$short_name), c(moders$long_name, stdyModers$long_name), warn_missing = F)
d <- round(max(abs(min(df$pred)), abs(max(df$pred))), 3)
# mini <- if(d > 2) .05 else 0-(d+(d/5))
# maxi <- if(d > 2) 2.05 else 0+d+(d/5)
# lim <- c(mini, maxi)
# brk <- if(d > 2) c(0, 1, 2) else{round(c(0-d-(d/10), 0, 0+d+(d/10)),2)}
# lab <- if(d > 2){c("0", "1", "2")} else{str_remove(c(round(0-d-(d/10),2), 0, round(0+d+(d/10),2)), "^0")}
titl <- if(mod == "none"){o} else {sprintf("%s: Personality x %s Simple Effects", o, m)}
# colnames(df)[colnames(df) == mod] <- "mod_value"
df <- df %>% unclass %>% data.frame
df$mod_value <- df[,mod]
df <- df %>% select(-all_of(mod)) %>% as_tibble
if(class(df$mod_value) %in% c("factor", "character")){df <- df %>% mutate(mod_fac = factor(mod_value))}
else {
if(mod == "age") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("-10 yrs", "M", "+10 yrs")))
else if(mod == "baseYear") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("1990", "200)0", "2010")))
else if(mod == "baseAge") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("50", "60", "70")))
else if(mod == "predInt") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-5, 0, 5), labels = c("-5 yrs", "5 yrs", "+5 yrs")))
else if(mod == "education") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-5, 0, 5), labels = c("-5 yrs", "12 years", "+5 yrs")))
else df <- df %>% mutate(mod_fac = factor(mod_value, levels = unique(mod_value), labels = c("-1 SD", "M", "+1 SD")))
}
lt <- c("dotted", "solid", "dashed")[1:length(unique(df$mod_fac))]
mini <- floor(min(df$pred)); maxi <- 10
df %>%
mutate(Trait = factor(Trait, levels = traits$short_name, labels = traits$long_name),
lower = ifelse(lower < mini, mini, lower),
upper = ifelse(upper > maxi, 10, upper)) %>%
ggplot(aes(x = p_value
, y = pred
, group = mod_fac)) +
scale_y_continuous(limits = c(mini,maxi)
, breaks = seq(mini, maxi, by = 2)
, labels = seq(mini, maxi, by = 2)) +
scale_linetype_manual(values = lt) +
geom_ribbon(aes(ymin = lower
, ymax = upper
, fill = mod_fac)
, alpha = .25) +
geom_line(aes(linetype = mod_fac)) +
labs(x = "Personality (POMP)"
, y = paste(o, "(POMP)")
, title = titl
, linetype = m
, fill = m
, subtitle = "Method 2B: Pooled Regression Using Random Effects") +
facet_wrap(~Trait, nrow = 3) +
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/2b_ipd_mlm/%s/figures/overall simple effects/%s_%s_%s.png", local_path, type, outcome, mod, cov), width = 6, height = 6)
}
ipd2b_se_plot <- nested_ipd2b_simp %>%
select(-pred.rx) %>%
group_by(type, Outcome, Moderator, Covariate) %>%
nest() %>%
ungroup() %>%
# filter(Moderator== "education") %>%
mutate(data = map(data, ~(.) %>% unnest(pred.fx)),
plot = pmap(list(data, Outcome, Moderator, type, Covariate), simp_eff_fun))4.2.2.7.2.4 Study-Specific Simple Effects
ipd2b_std_se_plot_fun <- function(df, outcome, trait, mod, cov, type){
print(paste(outcome, mod))
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)
d <- round(max(abs(min(df$pred)), abs(max(df$pred))), 3)
# mini <- if(d > 2) .05 else 0-(d+(d/5))
# maxi <- if(d > 2) 2.05 else 0+d+(d/5)
# lim <- c(mini, maxi)
# brk <- if(d > 2) c(0, 1, 2) else{round(c(0-d-(d/10), 0, 0+d+(d/10)),2)}
# lab <- if(d > 2){c("0", "1", "2")} else{str_remove(c(round(0-d-(d/10),2), 0, round(0+d+(d/10),2)), "^0")}
titl <- if(mod == "none"){sprintf("%s: %s", o, trt)} else {sprintf("%s: %s x %s Simple Effects", o, trt, m)}
# colnames(df)[colnames(df) == mod] <- "mod_value"
df <- df %>% mutate(study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin")))
df <- df %>% unclass %>% data.frame
df$mod_value <- df[,mod]
df <- df %>% select(-all_of(mod)) %>% as_tibble
if(class(df$mod_value) == "factor"){df <- df %>% mutate(mod_fac = factor(mod_value))}
if(class(df$mod_value) %in% c("factor", "character")){df <- df %>% mutate(mod_fac = factor(mod_value))}
else {
if(mod == "age") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("-10 yrs", "M", "+10 yrs")))
else if(mod == "baseYear") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("1990", "200)0", "2010")))
else if(mod == "baseAge") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-10, 0, 10), labels = c("50", "60", "70")))
else if(mod == "predInt") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-5, 0, 5), labels = c("-5 yrs", "5 yrs", "+5 yrs")))
else if(mod == "education") df <- df %>% mutate(mod_fac = factor(mod_value, levels = c(-5, 0, 5), labels = c("-5 yrs", "12 years", "+5 yrs")))
else df <- df %>% mutate(mod_fac = factor(mod_value, levels = unique(mod_value), labels = c("-1 SD", "M", "+1 SD")))
}
std <- unique(df$study)
cols <- (stdcolors %>% filter(studies %in% std))$colors
lt <- (stdcolors %>% filter(studies %in% std))$lt
ht <- length(unique(df$mod_fac))
mini <- floor(min(df$pred)); maxi <- 10
df %>%
mutate(study = factor(study, levels = stdcolors$studies),
lower = ifelse(lower < mini, mini, lower),
upper = ifelse(upper > maxi, 10, upper),
gr = ifelse(study == "Overall", "Overall", "study")) %>%
group_by(study, mod_fac, p_value, gr) %>%
summarize_at(vars(pred, lower, upper), mean) %>%
ungroup() %>%
ggplot(aes(x = p_value
, y = pred
, group = study)) +
scale_y_continuous(limits = c(mini,maxi)
, breaks = seq(mini, maxi, by = 2)
, labels = seq(mini, maxi, by = 2)) +
scale_linetype_manual(values = lt) +
scale_color_manual(values = cols) +
scale_fill_manual(values = cols) +
scale_size_manual(values = c(2,.8)) +
# geom_ribbon(aes(ymin = lower
# , ymax = upper
# , fill = study)
# , alpha = .25) +
geom_line(aes(linetype = study, color = study, size = gr)) +
labs(x = "Personality (POMP)"
, y = paste(o, "(POMP)")
, title = titl
, linetype = "Study"
, color = "Study"
, fill = "Study"
, subtitle = "Method 2B: Pooled Regression Using Random Effects") +
guides(size = "none") +
facet_wrap(~mod_fac, nrow = 1) +
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/2b_ipd_mlm/%s/figures/study specific simple effects/%s_%s_%s_%s.png", local_path, type, outcome, trt, mod, cov), width = 3*ht, height = 5)
}
nested_ipd2b_simp %>%
filter(!Moderator %in% stdyModers$short_name) %>%
# filter(map_lgl(pred.rx, ~!is.null(.)))
mutate(pred.fx = map(pred.fx, ~(.) %>% mutate(study = "Overall")),
comb.fx = map2(pred.fx, pred.rx, full_join)) %>%
select(-pred.fx, -pred.rx) %>%
# filter(Moderator %in% c("age", "education")) %>%
mutate(pmap(list(comb.fx, Outcome, Trait, Moderator, Covariate, type), ipd2b_std_se_plot_fun))load("~/Documents/projects/data synthesis/crystallized/results/2b_ipd_mlm/Frequentist/models/crystallized_C_scale_all.RData")
fixef(m)## (Intercept) p_value age genderFemale education
## 6.737893882 0.070935261 0.001625964 0.103872146 0.255833361
## scaleBFI-S scaleEysenck scaleIPIP NEO scaleMIDI scaleTDA-40
## 1.445698251 1.649520347 -0.586613200 -2.027132359 -1.593313674
## p_value:scaleBFI-S p_value:scaleEysenck p_value:scaleIPIP NEO p_value:scaleMIDI p_value:scaleTDA-40
## -0.066244543 -0.213937506 -0.115922216 0.012948027 -0.043546168cntrm <- rbind(
c(0,1,rep(0,13)) # NEO FFI
, c(0,1,rep(0,8),1,0,0,0,0) # BFI-S
, c(0,1,rep(0,9),1,0,0,0) # Eysenck
, c(0,1,rep(0,10),1,0,0) # IPIP NEO
, c(0,1,rep(0,11),1,0) # MIDI
, c(0,1,rep(0,12),1) # TDA-40
); rownames(cntrm) <- c("NEO-FFI", "BFI-S", "Eysenck", "IPIP NEO", "MIDI", "TDA-40")
(multcomp::glht(m, 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 = sprintf("b = %.2f, 95%% CI [%.2f, %.2f]", Estimate, lwr, upr))## Estimate lwr upr term est
## 1 0.070935261 0.01391155 0.12795897 NEO-FFI b = 0.07, 95% CI [0.01, 0.13]
## 2 0.004690717 -0.09200863 0.10139006 BFI-S b = 0.00, 95% CI [-0.09, 0.10]
## 3 -0.143002245 -0.24210542 -0.04389907 Eysenck b = -0.14, 95% CI [-0.24, -0.04]
## 4 -0.044986955 -0.13309875 0.04312484 IPIP NEO b = -0.04, 95% CI [-0.13, 0.04]
## 5 0.083883288 0.04542256 0.12234401 MIDI b = 0.08, 95% CI [0.05, 0.12]
## 6 0.027389093 -0.01026313 0.06504132 TDA-40 b = 0.03, 95% CI [-0.01, 0.07]load("~/Documents/projects/data synthesis/crystallized/results/2b_ipd_mlm/Frequentist/models/crystallized_A_baseYear_all.RData")
fixef(m)## (Intercept) p_value age genderFemale education baseYear p_value:baseYear
## 6.631036181 -0.019513289 0.002262725 0.103624853 0.274346670 -0.118874391 0.008497705cntrm <- c(
"p_value = 0" # 2000
, "p_value - 10*p_value:baseYear = 0" # 1990
, "p_value + 10*p_value:baseYear = 0" # 2010
); names(cntrm) <- c(2000, 1990, 2010)
(multcomp::glht(m, cntrm) %>% # multcomp hypothesis function
confint(., calpha = multcomp::univariate_calpha()))$confint %>%
data.frame() %>%
data.frame() %>%
rownames_to_column("cntr") %>%
mutate(term = names(cntrm)) %>%
select(-cntr) %>%
mutate(est = sprintf("b = %.2f, 95%% CI [%.2f, %.2f]", Estimate, lwr, upr))## Estimate lwr upr term est
## 1 -0.01951329 -0.04827826 0.009251684 2000 b = -0.02, 95% CI [-0.05, 0.01]
## 2 -0.10449034 -0.17448848 -0.034492203 1990 b = -0.10, 95% CI [-0.17, -0.03]
## 3 0.06546376 0.03595811 0.094969409 2010 b = 0.07, 95% CI [0.04, 0.09]4.2.2.8 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 fully pooled, one-stage multilevel models with individual participant data in 11 studies. Multilevel models allowed us to estimate overall, sample-specific, and heterogeneity estimates as well as both participant-level and sample-level moderators of Big Five personality characteristic-crystallized abilities associations. Heterogeneity estimates were estimated for both personality-cognitive ability associations (\(\tau_{11}^2\)) as well as for participant-level moderators (\(\tau_{22}^2\); see online materials). Fully adjusted prospective Big Five personality characteristic-crystallized abilities associations and participant- and sample-level moderators of these associations can be found in the fourth section of Table S4.
|
E
|
A
|
C
|
N
|
O
|
|
|---|---|---|---|---|---|
| Term | b [CI] | b [CI] | b [CI] | b [CI] | b [CI] |
| None | |||||
| Personality |
0.01 [-0.02, 0.05] |
0.00 [-0.03, 0.03] |
0.01 [-0.04, 0.07] |
-0.06 [-0.09, -0.03] |
0.16 [0.07, 0.25] |
| Age | |||||
| Age |
0.001 [-0.03, 0.03] |
0.000 [-0.09, 0.09] |
-0.001 [-0.01, 0.01] |
0.00 [-0.00, 0.01] |
0.001 [-0.09, 0.09] |
| Gender | |||||
| Gender (Male v Female) |
0.05 [0.01, 0.09] |
0.02 [-0.02, 0.06] |
0.01 [-0.03, 0.06] |
-0.02 [-0.05, 0.01] |
-0.000 [-0.06, 0.06] |
| Education | |||||
| Education (Years) |
-0.01 [-0.07, 0.06] |
-0.01 [-0.02, 0.01] |
0.01 [-0.02, 0.04] |
0.00 [-0.01, 0.01] |
-0.01 [-0.12, 0.09] |
| Continent | |||||
| Continent (North America v Europe) |
0.07 [0.00, 0.14] |
-0.12 [-0.26, 0.03] |
-0.12 [-0.24, 0.00] |
-0.03 [-0.06, 0.01] |
0.14 [-0.08, 0.36] |
| Continent (North America v Australia) |
0.03 [-0.05, 0.10] |
-0.01 [-0.18, 0.15] |
-0.02 [-0.16, 0.11] |
-0.11 [-0.14, -0.08] |
0.11 [-0.19, 0.40] |
| Country | |||||
| Country (United States v Germany) |
0.08 [-0.02, 0.18] |
-0.01 [-0.10, 0.08] |
-0.05 [-0.19, 0.09] |
-0.02 [-0.09, 0.05] |
0.12 [-0.15, 0.38] |
| Country (United States v Sweden) |
0.07 [-0.03, 0.16] |
-0.23 [-0.35, -0.12] |
-0.20 [-0.34, -0.06] |
-0.05 [-0.11, 0.01] |
0.22 [-0.13, 0.58] |
| Country (United States v The Netherlands) |
-0.02 [-0.06, 0.03] |
||||
| Country (United States v Australia) |
0.02 [-0.06, 0.11] |
-0.01 [-0.05, 0.03] |
-0.02 [-0.13, 0.08] |
-0.11 [-0.14, -0.08] |
0.11 [-0.21, 0.43] |
| Personality Scale | |||||
| Scale (NEO-FFI v DPQ) |
0.01 [-0.09, 0.11] |
||||
| Scale (NEO-FFI v Eysenck) |
0.09 [-0.01, 0.18] |
-0.24 [-9.45, 8.97] |
-0.21 [-0.33, -0.10] |
-0.02 [-0.11, 0.08] |
0.14 [-0.50, 0.77] |
| Scale (NEO-FFI v MIDI) |
-0.00 [-0.08, 0.08] |
-0.01 [-9.22, 9.20] |
0.01 [-0.06, 0.08] |
0.03 [-0.06, 0.12] |
-0.07 [-0.70, 0.55] |
| Scale (NEO-FFI v BFI-S) |
0.12 [0.01, 0.23] |
-0.01 [-9.22, 9.20] |
-0.07 [-0.18, 0.05] |
0.03 [-0.09, 0.14] |
-0.12 [-0.75, 0.51] |
| Scale (NEO-FFI v IPIP NEO) |
0.11 [-0.01, 0.22] |
0.01 [-9.20, 9.22] |
-0.12 [-0.22, -0.01] |
-0.00 [-0.12, 0.12] |
0.01 [-0.61, 0.64] |
| Scale (NEO-FFI v TDA-40) |
0.04 [-0.04, 0.13] |
-0.01 [-9.22, 9.20] |
-0.04 [-0.11, 0.02] |
-0.08 [-0.17, 0.01] |
0.02 [-0.60, 0.65] |
| Baseline Age | |||||
| Study Baseline Age |
-0.001 [-0.00, 0.00] |
0.00 [-0.00, 0.01] |
0.00 [-0.00, 0.01] |
0.001 [-0.001, 0.00] |
-0.000 [-0.01, 0.01] |
| Baseline Year | |||||
| Study Baseline Year |
-0.000 [-0.01, 0.00] |
0.01 [0.00, 0.01] |
0.01 [-0.00, 0.01] |
0.001 [-0.00, 0.01] |
-0.01 [-0.02, 0.00] |
| Prediction Interval | |||||
| Prediction Interval |
-0.00 [-0.01, 0.01] |
-0.01 [-0.02, -0.00] |
-0.01 [-0.02, 0.01] |
-0.00 [-0.01, 0.01] |
0.01 [-0.02, 0.04] |
Of the 20 key terms (five main effects and 15 moderators) tested at the participant-level, two (10%) were significant. Openness to experience was associated with later crystallized abilities across samples, and gender moderated the relationship between Extraversion and cognitive ability (b = 0.05, 95% CI [0.01, 0.09]). For the latter, we examined simple effects plots for gender moderating of extraversion and crystallized abilities. The black line in Figure S5 displays the fully adjusted association across genders. In this case, the interaction indicates that although the association was null for both males (b = -0.02, 95% CI [-0.04, 0.01]) and females (b = 0.03, 95% CI [-0.01, 0.08]), the associations did differ from one another.
Figure 4.3: Figure S5. 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.
Next, we examined the sample-specific effects to see how consistent associations were across samples. Figure S6 presents the forest plot of fully adjusted sample-specific and overall prospective associations between the Big Five personality characteristics and crystallized abilities. As shown in the forest plot, openness showed the most consistent association with crystallized ability, with five of seven samples being positively and significantly associated with crystallized abilities. Other prospective associations, such as those with extraversion, agreeableness, and conscientiousness, were less consistent, with at least one significant association in different directions across samples. For example, although there was no overall association between conscientiousness and crystallized abilities, SATSA had a significant negative association, and HILDA and HRS had positive associations.
Figure 4.4: Figure S6. 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.
Finally, we can also examine the consistency of simple effects across samples for each moderator. The thinner dashed and colored lines in Figure S5 display the fully adjusted predicted crystallized abilities levels at different levels of extraversion across samples for males and females. As noted above, the overall trend suggests a positive association for females and null association for males, such that women who were higher in extraversion tended to score better on crystallized / knowledge domain tasks and that these associations were quite consistent across samples. Women show a slight positive relationship in all samples but HRS and ROS. For men, the associations were also quite consistent, with all the samples showing significant negative associations but GSOEP, MAP, and OCTO-Twin, which showed positive associations.
We also examined sample-level moderators of the associations, which are displayed in Table S4 (fully adjusted) above. Across the 67 tested associations, 10 (14.93%) were significant. Some of these differences were a function of the personality scale, with for example, estimates of the association between Conscientiousness and crystallized / knowledge domain cognitive ability differing across scales. The NEO-FFI (b = 0.07, 95% CI [0.01, 0.13]) and MIDI (b = 0.08, 95% CI [0.05, 0.12] ) demonstrated positive associations; the IPIP NEO (b = -0.04, 95% CI [-0.13, 0.04]), BFI-S (b = 0.00, 95% CI [-0.09, 0.10]), and TDA-40 (b = 0.03, 95% CI [-0.01, 0.07]) demonstrated null associations; and the measure used in SATSA and OCTO-Twin (b = -0.14, 95% CI [-0.24, -0.04]) had a negative association. In addition, for Agreeableness, samples with earlier baseline years had more negative effects than samples with later baseline effects. For example, samples with baseline years around 1990 demonstrated a negative association (b = -0.10, 95% CI [-0.17, -0.03]), samples with baseline years around 2000 demonstrated a null association (b = -0.02, 95% CI [-0.05, 0.01]), and samples with baseline years around 2010 demonstrated a positive association (b = 0.07, 95% CI [0.04, 0.09]).