Chapter 6 Method 3: One-Stage Individual Participant Analyses Reported Together

6.1 Analytic Plan

In the present study, we estimate associations between the Big Five personality traits and crystallized cognitive abilities using a series of separate regression models for each sample (also known as a coordinated data analysis when multiple analysts are used). The procedure is as follows:

Models were run separately for each sample, personality characteristic, outcome, covariate, and moderator combination. To do so, we wrote a series of functions in the R programming language (see online materials) to (1) set up and run the model and extract model coefficients, (2) extract simple-effects predictions for moderator models (i.e. predicted values across levels of the moderator values), and (3) extract effect size metrics for the later meta-analysis. The basic form of the model is as follows:

\(Y_{ij}=b_0+b_1\ast predictor_{ij}+\epsilon_{ij}\) \(\epsilon_{ij}\sim\mathcal{N}(0, \sigma^2)\),

where \(b_1\) represents the effect of personality predicting the outcome separately in each sample.
The modeling function used the base R lm() function to run the models and the tidy() function from the broom package to extract model coefficients and confidence intervals (CI). Inferences will be made based on the 95% confidence intervals. Effect sizes and their standard errors were extracted from the model summary(). Simple-effects predictions were calculated by providing the full range of personality levels (0-10) and average levels of the covariates from the data used to estimate the model as the "newdata" argument in the base R predict() function.

6.2 Step 1: Combine Data

We again need to combine data. However, rather than combining data across studies, for the two-stage approach, we’ll be combining data within studies in order to run separate analyses for each before combining via meta-analytic tools.

loadRData <- function(fileName, type){
#loads an RData file, and returns it
    path <- sprintf("%s/data/clean/%s_cleaned.RData", local_path, fileName)
    load(path)
    get(ls()[grepl(type, ls())])
}

ipd4_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))

ipd4_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 rows

6.2.1 Harmonize Data

ipd4_reg_data <- ipd4_reg_data %>%
  group_by(study, Trait, Outcome) %>%
  mutate_at(vars(p_value, o_value, SRhealth), 
            ~((. - min(., na.rm = T))/(max(., na.rm = T) - min(., na.rm = T))*10)) %>%
  mutate(gender = factor(gender, levels = c(0,1), labels = c("Male", "Female")),
         education = education - 12,
         age = age - mean(age, na.rm = T))
ipd4_reg_data

## # A tibble: 119,597 × 13
## # Groups:   study, Trait, Outcome [40]
##    study Trait p_value p_year SID   Outcome      o_value o_year education gender SRhealth yearBrth    age
##    <chr> <chr>   <dbl>  <dbl> <chr> <chr>          <dbl>  <dbl>     <dbl> <fct>     <dbl>    <dbl>  <dbl>
##  1 BASE  E         7     1990 10004 crystallized   2.17    2000         5 Male        7.5     1918 -6.23 
##  2 BASE  E         7     1990 10010 crystallized   3.08    1995        -1 Male        2.5     1911  0.766
##  3 BASE  E         3     1990 10033 crystallized   0.909   1997        -4 Female     10       1910  1.77 
##  4 BASE  E         4.5   1990 10034 crystallized   6.54    1995        -2 Male       10       1914 -2.23 
##  5 BASE  E         3     1990 10111 crystallized   6.54    1995         1 Female      7.5     1901 10.8  
##  6 BASE  E         5     1990 10115 crystallized  10       1995        -1 Male        2.5     1918 -6.23 
##  7 BASE  E         4.5   1990 10116 crystallized   3.85    1995        -2 Male        7.5     1917 -5.23 
##  8 BASE  E         4.5   1990 10145 crystallized   6.36    1997        -2 Female     10       1897 14.8  
##  9 BASE  E         5.5   1990 10175 crystallized   2.31    1995        -4 Female      5       1911  0.766
## 10 BASE  E         5.5   1990 10188 crystallized   5.77    2004        -2 Female      2.5     1903  8.77 
## # ℹ 119,587 more rows

6.2.2 Save Data Files

save_fun <- function(d, trait, outcome, study){
  save(d, file = sprintf("%s/data/two_stage/%s_%s_%s.RData", local_path, trait, outcome, study))
}

ipd4_reg_data %>%
  group_by(study, Trait, Outcome) %>%
  nest() %>%
  ungroup() %>%
  mutate(data = pmap(list(data, Trait, Outcome, study), save_fun))

## # A tibble: 40 × 4
##    study Trait Outcome      data  
##    <chr> <chr> <chr>        <list>
##  1 BASE  E     crystallized <NULL>
##  2 BASE  N     crystallized <NULL>
##  3 BASE  O     crystallized <NULL>
##  4 EAS   A     crystallized <NULL>
##  5 EAS   C     crystallized <NULL>
##  6 EAS   E     crystallized <NULL>
##  7 EAS   N     crystallized <NULL>
##  8 EAS   O     crystallized <NULL>
##  9 GSOEP A     crystallized <NULL>
## 10 GSOEP C     crystallized <NULL>
## # ℹ 30 more rows

6.3 Step 2: Run Models for Each Study

6.3.1 Functions

6.3.1.1 Model Function

ipd4_study_mod_fun <- function(trait, outcome, type, mod, study, cov){
  ## load the data
  load(sprintf("%s/data/two_stage/%s_%s_%s.RData", local_path, trait, outcome, study))
  
  ## compiled Bayesian model to speed up processing and avoid crashing
  if(type == "Bayesian") load(sprintf("%s/results/4_ca_reptog/bayes_sample_mod.RData", local_path))
  
  ## model formula 
  if (cov == "all") cv <- c("age", "gender", "education")
  if (!cov %in% c("all", "none")) cv <- cov
  rhs <- "p_value"
  rhs <- if(cov != "none") c(rhs, cv) else rhs
  if(mod != "none"){rhs <- c(rhs, paste("p_value", mod, sep = "*"))}
  rhs <- paste(rhs, collapse = " + ")
  f <- paste("o_value ~ ", rhs, collapse = "")
  
  ## run the models & save
  m <- if(type == "Frequentist"){do.call("lm", list(formula = f, data = quote(d)))} else {update(m, formula = f, newdata = d, warmup = 1000, iter = 2000)}
  save(m, file = sprintf("%s/results/4_ca_reptog/%s/models/%s_%s_%s_%s_%s.RData"
                         , local_path, type, outcome, trait, mod, cov, study))
  
  ## extract model terms and confidence intervals & save
  fx <- tidy(m, conf.int = T) %>%
    select(term, estimate, conf.low, conf.high)
  save(fx, file = sprintf("%s/results/4_ca_reptog/%s/summary/%s_%s_%s_%s_%s.RData"
                          , local_path, type, outcome, trait, mod, cov, study))
  
  ## calculate simple effects  
  if(mod != "none"){
    pred.rx <- ipd4_study_simpeff_fun(m, mod, type)
    save(pred.rx, file = sprintf("%s/results/4_ca_reptog/%s/predicted/%s_%s_%s_%s_%s.RData"
                                 , local_path, type, outcome, trait, mod, cov, study))
  }
  
  ## clean up the local function environment
  rm(list = c("d", "f", "m", "fx", "es", "rhs"))
  gc()
}

6.3.1.2 Study-Specific Simple Effects Function

ipd4_study_simpeff_fun <- function(m, moder, type){
  d <- if(type == "Bayesian") m$data else m$model
  d <- d %>% select(-o_value, -p_value)
  cols <- colnames(d)
  md_cl <- class(d[,moder])
  if(any(sapply(d, class) == "numeric")){
    msd <- d %>%
      select_if(is.numeric) %>%
      pivot_longer(everything()
                   , names_to = "item"
                   , values_to = "value") %>%
      group_by(item) %>%
      summarize_at(vars(value), lst(mean, sd)) %>%
      ungroup()
  }
  if(any(sapply(d, class) == "factor")){
    fct_lev <- d %>% 
      select_if(is.factor) %>%
      summarize_all(~list(unique(.)))
  }
  d <- d %>% select(-one_of(moder))
  md_levs <- if(md_cl == "numeric"){
    if(moder %in% c("age", "baseAge", "baseYear")) {
      c(-10, 0, 10)
      } else if (moder %in% c("predInt", "education")) {
      c(-5, 0, 5) 
      } else {
        with(msd, c(mean[item == moder] - sd[item == moder], mean[item == moder], mean[item == moder] + sd[item == moder]))
      }
  } else { 
    unique(fct_lev[,moder][[1]])
  }
  
  mod_frame <- crossing(
    p_value = seq(0,10,.5)
    , modvalue = md_levs
    ) %>% setNames(c("p_value", moder))
  
  if(ncol(d) > 0){
    if(any(sapply(d, class) == "numeric")){
      mod_frame <- tibble(mod_frame, d %>% select_if(is.numeric) %>% summarize_all(mean))
    } 
    if(any(sapply(d, class) == "factor")){
      mod_frame <- tibble(mod_frame, d %>% select_if(is.factor) %>% summarize_all(~levels(.)[1]))
    }
  }
  
  pred.fx <- if(type == "Bayesian"){
    bind_cols(
      mod_frame, 
      fitted(m, newdata = mod_frame) %>% data.frame
    ) %>%
      select(one_of(colnames(m$data)), pred = Estimate, lower = Q2.5, upper = Q97.5)
  } else {
    bind_cols(
      mod_frame, 
      predict(m, newdata = mod_frame, interval = "confidence") %>% data.frame 
    ) %>%
      select(one_of(colnames(m$model)), pred = fit, lower = lwr, upper = upr)
  }
  
  rm(list = c("m", "mod_frame", "d", "md_levs"))
  gc()
  return(pred.fx)
}

6.3.2 Run Models

load(sprintf("%s/data/two_stage/N_crystallized_HILDA.RData", local_path))

# clean data & keep only needed columns and a subset of the used variables
d <- d %>%
  filter(row_number() %in% sample(1:nrow(.), 100, replace = F))

# 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(o_value ~ p_value + age + gender + education)
m <- brm(formula = f
            , data = d
            , prior = Prior
            , iter = Iter
            , warmup = Warmup)

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: WARNING: There aren't enough warmup iterations to fit the
## Chain 1:          three stages of adaptation as currently configured.
## Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
## Chain 1:          the given number of warmup iterations:
## Chain 1:            init_buffer = 3
## Chain 1:            adapt_window = 16
## Chain 1:            term_buffer = 2
## Chain 1: 
## Chain 1: Iteration:  1 / 30 [  3%]  (Warmup)
## Chain 1: Iteration:  3 / 30 [ 10%]  (Warmup)
## Chain 1: Iteration:  6 / 30 [ 20%]  (Warmup)
## Chain 1: Iteration:  9 / 30 [ 30%]  (Warmup)
## Chain 1: Iteration: 12 / 30 [ 40%]  (Warmup)
## Chain 1: Iteration: 15 / 30 [ 50%]  (Warmup)
## Chain 1: Iteration: 18 / 30 [ 60%]  (Warmup)
## Chain 1: Iteration: 21 / 30 [ 70%]  (Warmup)
## Chain 1: Iteration: 22 / 30 [ 73%]  (Sampling)
## Chain 1: Iteration: 24 / 30 [ 80%]  (Sampling)
## Chain 1: Iteration: 27 / 30 [ 90%]  (Sampling)
## Chain 1: Iteration: 30 / 30 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.002 seconds (Warm-up)
## Chain 1:                0 seconds (Sampling)
## Chain 1:                0.002 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: WARNING: There aren't enough warmup iterations to fit the
## Chain 2:          three stages of adaptation as currently configured.
## Chain 2:          Reducing each adaptation stage to 15%/75%/10% of
## Chain 2:          the given number of warmup iterations:
## Chain 2:            init_buffer = 3
## Chain 2:            adapt_window = 16
## Chain 2:            term_buffer = 2
## Chain 2: 
## Chain 2: Iteration:  1 / 30 [  3%]  (Warmup)
## Chain 2: Iteration:  3 / 30 [ 10%]  (Warmup)
## Chain 2: Iteration:  6 / 30 [ 20%]  (Warmup)
## Chain 2: Iteration:  9 / 30 [ 30%]  (Warmup)
## Chain 2: Iteration: 12 / 30 [ 40%]  (Warmup)
## Chain 2: Iteration: 15 / 30 [ 50%]  (Warmup)
## Chain 2: Iteration: 18 / 30 [ 60%]  (Warmup)
## Chain 2: Iteration: 21 / 30 [ 70%]  (Warmup)
## Chain 2: Iteration: 22 / 30 [ 73%]  (Sampling)
## Chain 2: Iteration: 24 / 30 [ 80%]  (Sampling)
## Chain 2: Iteration: 27 / 30 [ 90%]  (Sampling)
## Chain 2: Iteration: 30 / 30 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.004 seconds (Warm-up)
## Chain 2:                0 seconds (Sampling)
## Chain 2:                0.004 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 3e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: WARNING: There aren't enough warmup iterations to fit the
## Chain 3:          three stages of adaptation as currently configured.
## Chain 3:          Reducing each adaptation stage to 15%/75%/10% of
## Chain 3:          the given number of warmup iterations:
## Chain 3:            init_buffer = 3
## Chain 3:            adapt_window = 16
## Chain 3:            term_buffer = 2
## Chain 3: 
## Chain 3: Iteration:  1 / 30 [  3%]  (Warmup)
## Chain 3: Iteration:  3 / 30 [ 10%]  (Warmup)
## Chain 3: Iteration:  6 / 30 [ 20%]  (Warmup)
## Chain 3: Iteration:  9 / 30 [ 30%]  (Warmup)
## Chain 3: Iteration: 12 / 30 [ 40%]  (Warmup)
## Chain 3: Iteration: 15 / 30 [ 50%]  (Warmup)
## Chain 3: Iteration: 18 / 30 [ 60%]  (Warmup)
## Chain 3: Iteration: 21 / 30 [ 70%]  (Warmup)
## Chain 3: Iteration: 22 / 30 [ 73%]  (Sampling)
## Chain 3: Iteration: 24 / 30 [ 80%]  (Sampling)
## Chain 3: Iteration: 27 / 30 [ 90%]  (Sampling)
## Chain 3: Iteration: 30 / 30 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.003 seconds (Warm-up)
## Chain 3:                0 seconds (Sampling)
## Chain 3:                0.003 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: WARNING: There aren't enough warmup iterations to fit the
## Chain 4:          three stages of adaptation as currently configured.
## Chain 4:          Reducing each adaptation stage to 15%/75%/10% of
## Chain 4:          the given number of warmup iterations:
## Chain 4:            init_buffer = 3
## Chain 4:            adapt_window = 16
## Chain 4:            term_buffer = 2
## Chain 4: 
## Chain 4: Iteration:  1 / 30 [  3%]  (Warmup)
## Chain 4: Iteration:  3 / 30 [ 10%]  (Warmup)
## Chain 4: Iteration:  6 / 30 [ 20%]  (Warmup)
## Chain 4: Iteration:  9 / 30 [ 30%]  (Warmup)
## Chain 4: Iteration: 12 / 30 [ 40%]  (Warmup)
## Chain 4: Iteration: 15 / 30 [ 50%]  (Warmup)
## Chain 4: Iteration: 18 / 30 [ 60%]  (Warmup)
## Chain 4: Iteration: 21 / 30 [ 70%]  (Warmup)
## Chain 4: Iteration: 22 / 30 [ 73%]  (Sampling)
## Chain 4: Iteration: 24 / 30 [ 80%]  (Sampling)
## Chain 4: Iteration: 27 / 30 [ 90%]  (Sampling)
## Chain 4: Iteration: 30 / 30 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.002 seconds (Warm-up)
## Chain 4:                0 seconds (Sampling)
## Chain 4:                0.002 seconds (Total)
## Chain 4:

save(m, file = sprintf("%s/results/4_ca_reptog/bayes_sample_mod.RData", local_path))
rm(list = c("d", "Prior", "Iter", "Warmup", "treedepth", "f", "m"))

# done <- tibble(file = list.files(sprintf("%s/results/4_ca_reptog/Bayesian/studyModels", local_path))) %>%
#   separate(file, c("Outcome", "Trait", "Moderator", "Covariate", "study"), sep = "_") %>%
#   mutate(study = str_remove_all(study, ".RData"),
#          done = "done")
plan(multisession(workers = 12L))
nested_ipd_reg <- tibble(files = list.files(sprintf("%s/data/two_stage", local_path))) %>%
  separate(files, c("Trait", "Outcome", "study"), sep = "_") %>%
  filter(!is.na(study)) %>%
  mutate(study = str_remove(study, ".RData")) %>%
  left_join(
    crossing(
      study = studies
      , Trait = traits$short_name
      , Outcome = outcomes$short_name
      , type = c("Frequentist", "Bayesian")
      , Moderator = c("age", "gender", "education")
      , Covariate = c("none", "all")
    ) %>% 
      full_join(
        crossing(
          study = studies
          , Trait = traits$short_name
          , Outcome = outcomes$short_name
          , type = c("Frequentist", "Bayesian")
          , Moderator = "none"
          , Covariate = c("none", "age", "gender", "education", "all")
          )
        )
    ) %>%
  filter(!is.na(Covariate)) %>%
  # full_join(done) %>% filter(is.na(done)) %>%
  # filter(type != "Frequentist") %>%
  filter(Trait == "N") %>%
  mutate(run = 
           # pmap(list(Trait, Outcome, type, Moderator, study, Covariate)
           future_pmap(list(Trait, Outcome, type, Moderator, study, Covariate)
                    , ipd4_study_mod_fun
                    , .progress = T
             , .options = furrr_options(
                                    globals = c("ipd4_study_simpeff_fun"
                                                , "read_path"
                                                , "local_path"
                                                , "local_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()

6.3.3 Compile Results

loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
    path <- sprintf("%s/results/4_ca_reptog/%s/%s/%s", local_path, type, folder, fileName)
    # print(path)
    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 <- nrow(d)
  return(n)
}

## load in effect size data 
## first get file names
nested_ipd4_ca <- tibble(type = c("Frequentist", "Bayesian")) %>%
  mutate(file = map(type, ~list.files(sprintf("%s/results/4_ca_reptog/%s/summary", local_path, .)))) %>%
  unnest(file) %>%
  separate(file, c("Outcome", "Trait", "Moderator", "Covariate", "study"), sep = "_", remove = F) %>% 
  filter(!is.na(study)) %>%
  ## read in the files
  mutate(study = str_remove(study, ".RData"),
         fx = map2(file, type, ~loadRData(.x, .y, "fx", "summary")),
         n = map2_dbl(file, type, n_fun)) %>%
  select(-file) %>%
  unnest(fx)

This results in a nested data frame, with columns:

• studyEff = standardized study-specific effects from stage 1 regressions

nested_ipd4_ca

## # A tibble: 4,120 × 11
##    type        Outcome      Trait Moderator Covariate study term         estimate conf.low conf.high     n
##    <chr>       <chr>        <chr> <chr>     <chr>     <chr> <chr>           <dbl>    <dbl>     <dbl> <dbl>
##  1 Frequentist crystallized A     age       all       EAS   (Intercept)   5.67      4.99     6.35      788
##  2 Frequentist crystallized A     age       all       EAS   p_value       0.0812   -0.0121   0.174     788
##  3 Frequentist crystallized A     age       all       EAS   age           0.0136   -0.110    0.138     788
##  4 Frequentist crystallized A     age       all       EAS   genderFemale -0.445    -0.774   -0.116     788
##  5 Frequentist crystallized A     age       all       EAS   education     0.233     0.185    0.281     788
##  6 Frequentist crystallized A     age       all       EAS   p_value:age  -0.0103   -0.0276   0.00702   788
##  7 Frequentist crystallized A     age       all       GSOEP (Intercept)   8.08      7.67     8.48      647
##  8 Frequentist crystallized A     age       all       GSOEP p_value       0.0196   -0.0343   0.0736    647
##  9 Frequentist crystallized A     age       all       GSOEP age           0.00909  -0.0209   0.0391    647
## 10 Frequentist crystallized A     age       all       GSOEP genderFemale  0.0439   -0.140    0.227     647
## # ℹ 4,110 more rows

6.3.3.0.1 Tables

6.3.3.0.1.1 Study-Specific

We’ll make a table of the results, separately for each moderator. To do this efficiently, we’ll make a function that creates those tables across all combinations. Before calling that function, we’ll do some reformatting and reshaping to get the data ready.

ipd4_res_tab <- nested_ipd4_ca %>%
  filter((Moderator == "none" & term == "p_value") | (Moderator != "none" & grepl("p_value:", term))) %>%
  mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns")) %>% # significance marker
  mutate_at(vars(estimate:conf.high), 
            ~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f",.))) %>%
  mutate(est = sprintf("%s<br>[%s, %s]", estimate, conf.low, conf.high),
         est = ifelse(sig == "sig", sprintf("<strong>%s</strong>", est), est),
         study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin")),
         study = factor(study, c(studies_long, "Meta-Analytic")),
         Outcome = factor(Outcome, outcomes$short_name, outcomes$long_name)) %>%
  select(type:Covariate, term, study, est) %>%
  pivot_wider(names_from = "Trait", values_from = "est") %>%
  select(type:study, E, A, C, N, O) %>%
  arrange(type, Outcome, Moderator, Covariate, study)
ipd4_res_tab

## # A tibble: 242 × 11
##    type     Outcome              Moderator Covariate term        study     E                         A      C     N     O    
##    <chr>    <fct>                <chr>     <chr>     <chr>       <fct>     <chr>                     <chr>  <chr> <chr> <chr>
##  1 Bayesian Crystallized Ability age       all       p_value:age BASE      0.002<br>[-0.03, 0.04]    <NA>   <NA>  0.01… -0.0…
##  2 Bayesian Crystallized Ability age       all       p_value:age EAS       -0.005<br>[-0.02, 0.01]   -0.01… <str… <str… -0.0…
##  3 Bayesian Crystallized Ability age       all       p_value:age GSOEP     -0.000<br>[-0.004, 0.003] -0.00… 0.00… 0.00… -0.0…
##  4 Bayesian Crystallized Ability age       all       p_value:age HILDA     0.001<br>[-0.000, 0.003]  0.000… <str… <str… <str…
##  5 Bayesian Crystallized Ability age       all       p_value:age HRS       -0.002<br>[-0.005, 0.001] -0.00… <str… 0.00… <str…
##  6 Bayesian Crystallized Ability age       all       p_value:age LASA      <NA>                      <NA>   <NA>  -0.0… <NA> 
##  7 Bayesian Crystallized Ability age       all       p_value:age MAP       -0.002<br>[-0.009, 0.005] <NA>   0.00… -0.0… <NA> 
##  8 Bayesian Crystallized Ability age       all       p_value:age MARS      <NA>                      <NA>   <NA>  -0.0… <NA> 
##  9 Bayesian Crystallized Ability age       all       p_value:age OCTO-Twin -0.01<br>[-0.06, 0.03]    <NA>   <NA>  0.03… <NA> 
## 10 Bayesian Crystallized Ability age       all       p_value:age ROS       0.004<br>[-0.006, 0.01]   <stro… 0.00… <str… <str…
## # ℹ 232 more rows

ipd4_std_table_fun <- function(d, type, moder, cov){
  cv <- mapvalues(cov, covars$short_name, covars$long_name, warn_missing = F)
  if(!grepl("djust", cv)) cv <- paste(cv, "Adjusted")
  md <- mapvalues(moder, c(moders$short_name, stdyModers$short_name),
                  c(moders$long_name, stdyModers$long_name), warn_missing = F)
  rs <- d %>% group_by(Outcome) %>% tally() %>% 
    mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
  cs <- rep(1,6); names(cs) <- c(" ", paste0("<strong>", traits$short_name, "</strong>"))
  cap <- if(moder == "none") {
      sprintf("<strong>Table X.</strong><br><em>Method 4 Coordinated Analyses Reported Together: Sample Estimates of %s Personality-Cognitive Domain Relationships</em>", cv) 
    } else {
      sprintf("<strong>Table X.</strong><br><em>Method 4 Coordinated Analyses Reported Together: Sample %s Moderation of %s Personality-Cognitive Domain Relationships</em>", md, cv)
    }
  tab <- d %>%
    select(-Outcome) %>%
    kable(., "html"
          , escape = F
          , booktabs = T
          , col.names = c("Study", rep("b [CI]", 5))
          , align = c("r", rep("c", 5))
          , 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/4_ca_reptog/%s/tables/study specific/%s_%s.html",
                                 local_path, type, moder, cov))
  return(tab)
}

ipd4_std_tab <- ipd4_res_tab %>%
  filter(!Moderator %in% stdyModers$short_name) %>%
  select(-term) %>%
  group_by(type, Moderator, Covariate) %>%
  nest() %>%
  ungroup() %>%
  mutate(tab = pmap(list(data, type, Moderator, Covariate), ipd4_std_table_fun))

6.3.3.0.1.2 Meta-Analytic

Results are not meta-analyzed, so there are no meta-analytic results.

6.3.3.0.1.3 All Model Terms

ipd4_mod_tab <- nested_ipd4_ca %>%
  # 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, moders$short_name, moders$long_name),
         Covariate = factor(Covariate, covars$short_name, str_wrap(covars$long_name, 15)),
         term = str_replace_all(term, "metamod", paste0("p_value:", Moderator))) %>%
  # format values as text, combine estimates and CI's, bold significance
  mutate_at(vars(estimate, conf.low, conf.high), 
            ~ifelse(abs(.) < .01, sprintf("%.3f", .), sprintf("%.2f", .))) %>%
  mutate(est = sprintf("%s<br>[%s, %s]", estimate, conf.low, conf.high),
         est = ifelse(sig == "sig", sprintf("<strong>%s</strong>", est), est)) %>%
  # mutate(est = sprintf("%s [%s, %s]", estimate, conf.low, conf.high),
  # est = ifelse(sig == "sig", sprintf("\\textbf{%s}", est), est)) %>%
  # final reshaping, remove extra columns, arrange values, and change to wide format
  select(-estimate, -conf.low, -conf.high, -sig, -n) %>%
  arrange(type, Outcome, Trait, Moderator, Covariate) %>%
  pivot_wider(names_from = "Trait", values_from = "est") 

ipd4_mod_tab_fun <- function(d, type, out, moder, cov){
  md <- mapvalues(moder, c(moders$long_name, stdyModers$long_name), c(moders$short_name, stdyModers$short_name), warn_missing = F)
  o <- mapvalues(out, outcomes$long_name, outcomes$short_name, warn_missing = F)
  cv <- mapvalues(cov, covars$long_name, covars$short_name, warn_missing = F)
  cs <- rep(1,6)
  names(cs) <- c(" ", traits$long_name)
  # cln <- if(length(unique(d$term2)) == 1) c("Covariate", rep("\\textit{b} [CI]", 5)) else c(" ", "Term", rep("\\textit{b} [CI]", 5))
  cln <- c("Term", rep("<em>b</em> [CI]", 5))
  al <- c("r", rep("c", 5))
  # caption 
  cap <- if(md == "none") "4 Coordinated Analyses Reported Together: All Model Estimates of Fixed Effect Personality-Crystallized Domain Associations" else sprintf("4 Coordinated Analyses Reported Together: All Model Estimates of Fixed Effect %s Moderation of Personality-Crystallized Domain Associations", md)
  
  d <- d %>% arrange(study, term)
  
  rs <- d %>% group_by(study) %>% tally() %>% 
    mutate(end = cumsum(n), start = lag(end) + 1, start = ifelse(is.na(start), 1, start))
  
  # kable the table
  tab <- d %>%
    select(-study) %>%
    kable(., "html"
    # kable(., "latex"
          , booktabs = T
          , escape = F
          , col.names = cln
          , align = al
          , caption = cap
    ) %>% 
    kable_classic(full_width = F, html_font = "Times New Roman") %>%
    # kable_styling(full_width = F, font_size = 7) %>%
    add_header_above(cs)
  
  for(i in 1:nrow(rs)){
    tab <- tab %>% kableExtra::group_rows(rs$study[i], rs$start[i], rs$end[i])
  }
  # save the resulting html table
  save_kable(tab, file = sprintf("%s/results/4_ca_reptog/%s/tables/all terms/%s-%s-%s.html"
                                 , local_path, type, o, md, cv))
  return(tab) # return the html table
}

ipd4_mod_tab2 <- ipd4_res_tab %>%
  group_by(type, Outcome, Moderator, Covariate) %>%
  nest() %>%
  ungroup() %>%
  mutate(tab = pmap(list(data, type, Outcome, Moderator, Covariate), ipd4_mod_tab_fun)) 

6.3.4 Figures

6.3.4.1 Overall Forest

Results are not meta-analyzes, so there is no overall estimate.

6.3.4.2 Study-Specific Forest

ipd4_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))
  df <- df %>% arrange(estimate)
  stds <- df$study[df$study != " "]
  df <- df %>%
    mutate(study = factor(study, rev(c(" ", stds)))
           # , 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 = "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) + .5, 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) + .75, y = "est", hjust = .5, vjust = 0) +
      annotate("text", label = "N", x = length(stds) + .75, y = "n", hjust = .5, vjust = 0) +
      # geom_vline(aes(xintercept = .5)) +
      geom_vline(aes(xintercept = length(stds) + .5)) +
      coord_flip() +
      theme_classic() +
      theme(plot.title = element_text(face = "bold", hjust = 0)
            , axis.text = element_blank()
            , axis.ticks = element_blank()
            , axis.title = element_blank()
            , axis.line.y = element_blank())
  
  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
      )

  p3 <- cowplot::plot_grid(p1, p2
                     , rel_widths = c(.4, .6)
                     , align = "h"
                     )
  # p <- cowplot::plot_grid(ttl, subttl, p3, rel_heights = c(.05, .05, .9), nrow = 3)
  p <- cowplot::plot_grid(ttl, p3, rel_heights = c(.05, .95), nrow = 2)
  gc()
  save(p
       , file = sprintf("%s/results/4_ca_reptog/%s/figures/study specific forest/rdata/%s_%s_%s_%s.RData", local_path, type, outcome, trait, mod, cov))
  return(p)
}

## fixed effects
nested_ipd4_reg_fp <- nested_ipd4_ca %>%
  filter(Moderator %in% moders$short_name) %>%
  ## filter key terms
  filter((Moderator == "none" & term == "p_value")|
         (Moderator != "none" & grepl("^p_value:", term) & !grepl("study", term))) %>%
  ## significance
  mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns")
         , study = mapvalues(study, studies_long, studies_sp, warn_missing = F)) %>%
  ## grouping for plotting
  group_by(Outcome, Moderator, type, Covariate, Trait) %>%
  nest() %>%
  ungroup() %>%
  mutate(p = pmap(list(data, Outcome, Moderator, type, Covariate, Trait), ipd4_rx_plot_fun))

ipd4_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 4: Coordinated Analyses Reported Together",
          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/4_ca_reptog/%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/4_ca_reptog/%s/figures/study specific forest/%s_%s_%s.pdf", local_path, type, outcome, mod, cov)
         , width = 10, height = 10)
  return(T)
}

nested_ipd4_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), ipd4_rx_plot_comb_fun))

ipd4_rx_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$estimate)), abs(max(df$estimate))), 3)
  lim <- if(mod == "none"){c(-.25, .25)} else{c(0-d-(d/2.5), 0+d+(d/2.5))}
  brk <- if(mod == "none"){seq(-.2,.2,.2)} else{round(c(0-d-(d/5), 0, 0+d+(d/5)),2)}
  lab <- if(mod == "none"){c("-.2", "0", ".2")} else{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"){o} else {sprintf("%s: Personality x %s,", o, m)}
  titl <- if(!cov %in% c("none", "all")) paste(titl, cv, "Adjusted", collapse = " ") else paste(titl, cv, 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),
           study = factor(study, levels = c("Overall", studies_long)),
           Trait = factor(Trait, levels = traits$short_name, labels = traits$long_name),
           type = ifelse(study == "Overall", "fixed", "random")) %>%
    ggplot(aes(x = study, y = estimate)) +
      scale_y_continuous(limits = lim, breaks = brk, labels = lab) +
      scale_size_manual(values = c(2.5, 1.5)) +
      scale_shape_manual(values = shapes) + 
      scale_color_manual(values = c("blue", "black")) +
      scale_linetype_manual(values = lt) +
      geom_hline(aes(yintercept = 0), linetype = "dashed") +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high, linetype = term)
                    , width = 0
                    , position = position_dodge(width = .9)) + 
      geom_point(aes(color = type, size = type, shape = term)
                    , position = position_dodge(width = .9)) +
      labs(x = NULL
           , y = "Estimate (POMP)"
           , title = titl
           , subtitle = "Method 2B: Pooled Regression Using Random Effects"
         ) +
      facet_wrap(~Trait, scales = "free_y", nrow = 3) +
      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),
            strip.background = element_rect(fill = "black"),
            strip.text = element_text(face = "bold", color = "white"),
            axis.text = element_text(color = "black"))
  ht <- length(unique(df$study))
  local_path <- length(unique(df$Trait))
  ggsave(p, file = sprintf("%s/results/4_ca_reptog/%s/figures/study specific forest/%s_%s_fixed.png"
                          , local_path, type, mod, cov)
         , width = local_path*2
         , height = ht*.75)
  rm(p)
  gc()
  return(T)
}

## fixed effects
nested_ipd4_ca %>% 
    mutate(term = ifelse(is.na(term), names, term)) %>%
    # rename(names = term)
  ## filter key terms
      filter(Moderator %in% moders$short_name) %>%
      filter((Moderator == "none" & term == "p_value")|
             (Moderator != "none" & grepl("^p_value:", term))) %>%
  ## significance
  mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns")) %>%
  ## grouping for plotting
  group_by(Outcome, Moderator, type, Covariate) %>%
  nest() %>%
  ungroup() %>%
  mutate(pmap(list(data, Outcome, Moderator, type, Covariate), ipd4_rx_plot_fun))

6.3.4.3 Overall Simple Effects Plots

There are no overall effects when coordinated analyses are reported together without random effects meta-analysis.

loadRData <- function(fileName, type, obj, folder){
#loads an RData file, and returns it
    path <- sprintf("%s/results/4_ca_reptog/%s/%s/%s", local_path, type, folder, fileName)
    # print(path)
    load(path)
    get(ls()[grepl(obj, ls())])
}

## load in effect size data 
## first get file names
nested_ipd4_ca <- tibble(type = c("Frequentist", "Bayesian")) %>%
  mutate(file = map(type, ~list.files(sprintf("%s/results/4_ca_reptog/%s/predicted", local_path, .)))) %>%
  unnest(file) %>%
  separate(file, c("Outcome", "Trait", "Moderator", "Covariate", "study"), sep = "_", remove = F) %>% 
  filter(!is.na(study)) %>%
  ## read in the files
  mutate(study = str_remove(study, ".RData"),
         pred.rx = map2(file, type, ~loadRData(.x, .y, "pred.rx", "predicted"))) %>%
  select(-file) %>%
  unnest(pred.rx) %>%
  group_by(type, Trait, Outcome, Moderator, Covariate) %>%
  nest(pred.rx = study:upper) %>%
  ungroup()

6.3.4.4 Study-Specific Simple Effects Plots

ipd4_std_se_plot_fun <- function(d, outcome, trait, mod, cov, type){
  # print(paste(int, mod, cov, random, imp))
  o <- mapvalues(outcome, outcomes$short_name, outcomes$long_name, warn_missing = F)
  trt <- mapvalues(trait, traits$short_name, traits$long_name, warn_missing = F)
  cv <- mapvalues(cov, covars$short_name, covars$long_name, warn_missing = F)
  m <- mapvalues(mod, moders$short_name, moders$long_name, warn_missing = F)
  titl <- if(mod == "none"){sprintf("%s: %s", o, trt)} else {sprintf("%s: %s x %s Simple Effects", o, trt, m)}
  d <- d %>% mutate(study = mapvalues(study, c("BASE-I", "OCTO-TWIN"), c("BASE", "OCTO-Twin")),
                    study = factor(study, levels = stdcolors$studies))
  std <- unique(d$study)
  cols <- (stdcolors %>% filter(studies %in% std))$colors
  lt <- (stdcolors %>% filter(studies %in% std))$lt
  d <- d %>% unclass %>% data.frame
  d$mod_value <- d[,mod]
  d <- d %>% select(-all_of(mod)) %>% as_tibble
  d <- if(class(d$mod_value) %in% c("factor", "character")){
    d %>% mutate(mod_fac = factor(mod_value))
    } else{
      d %>%
      full_join(
        d %>% select(study, mod_value) %>% distinct() %>% arrange(study, mod_value) %>%
          group_by(study) %>% mutate(mod_fac = factor(c("-1 SD", "M", "+1 SD"), levels = c("-1 SD", "M", "+1 SD"))) %>% ungroup()
        )
  }
  ht <- length(unique(d$mod_fac))
  p <- d %>% 
    ggplot(aes(x = p_value, y = pred)) + 
      geom_line(aes(color = study
                    , group = interaction(study, mod_fac)
                    , linetype = study)
                , size = 1) + 
      # geom_ribbon(aes(fill = study
      #                 , group = interaction(study, mod_fac)
      #                 , ymin = lower
      #                 , ymax = upper)
      #             , alpha = .25) + 
      # stat_smooth(aes(weight = n)
      #             , method = "lm"
      #             , formula = y~x
      #             , size = 1.2
      #             , color = "darkslateblue"
      #             ) +
      scale_y_continuous(limits = c(0,10)
                         , breaks = c(0,5,10)
                         , labels = c(0,5,10)) +
      scale_color_manual(values = cols) +
      scale_linetype_manual(values = lt) +
      labs(x = "Personality Score (POMP)"
           , y = "Cognition Score (POMP)"
           , color = "Study"
           , fill = "Study"
           , linetype = "Study"
           , title = titl
           , subtitle = "Method 4: Coordinated Analyses Reported Together"
           )  +
      facet_wrap(~mod_fac) + 
      theme_classic() + 
      theme(legend.position = "bottom"
            , plot.title = element_text(face = "bold", hjust = .5, size = rel(.95))
            , plot.subtitle = element_text(size = rel(1.1), hjust = .5)
            , strip.background = element_rect(fill = "black")
            , strip.text = element_text(face = "bold", color = "white")
            , axis.text = element_text(color = "black"))
  ggsave(file = sprintf("%s/results/4_ca_reptog/%s/figures/study specific simple effects/%s_%s_%s_%s.png", local_path, type, outcome, trt, mod, cov), width = 3*ht, height = 5)
}

nested_ipd4_ca %>%
  mutate(p = pmap(list(pred.rx, Outcome, Trait, Moderator, Covariate, type), ipd4_std_se_plot_fun))

6.4 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 separate regression models of individual participant data of 11 studies. Fully adjusted Big Five personality trait-crystallized ability associations as well as participant- and sample-level moderators of these associations can be found in the forest plot in Figure S8. Of all estimates, openness was most consistently associated with later crystallized abilities, with positive associations in six of seven studies (all except ROS). Neuroticism (-) was somewhat consistently associated with crystallized abilities, demonstrating a negative association in five of the 11 samples. In contrast, extraversion and conscientiousness were less consistently associated with crystallized abilities with significant effects for at least one study in both directions. For extraversion, only two of nine studies were significant: a positive association in GSOEP and a negative association in HRS. Similarly, while there was a positive association between conscientiousness and crystallized abilities in HRS, MAP, and ROS but a negative association in SATSA.

Figure 6.1: Figure S8. Forest plot of fully adjusted prospective associations between Big Five personality characteristics and crystallized abilities across samples for using one-stage pooled integrative data analysis via effects coded contrasts. Overall point estimates (squares) represent the grand-mean estimates of the association across samples, while sample point estimates represent regression terms or linear combinations of regression terms. Error bars capture the 95% CI around the point estimate. Arrows indicate the confidence band was truncated to better visualize the estimates.

Next, we examined participant-level moderators of the association between personality traits and crystallized abilities. Moderator associations in each sample can be seen in Table S6. For age, 12 of the 40 tested moderators were significant. Of these, the most consistent were for conscientiousness (3/7) and openness (3/7). However, for each of these, the overall pattern was less clear because moderator associations were in both directions. For example, for openness (see Figure S9), age negatively moderated the association between openness and crystallized abilities in HILDA and HRS such that associations were weaker for younger relative to older participants. In contrast, age positively moderated the association in ROS such that associations were stronger for older relative to younger participants.

Table 6.1: Table X.
Method 4 Coordinated Analyses Reported Together: Sample Age Moderation of Fully Adjusted Personality-Cognitive Domain Relationships

	E	A	C	N	O
Study	b [CI]	b [CI]	b [CI]	b [CI]	b [CI]
Crystallized Ability
BASE	0.002 [-0.03, 0.04]			0.01 [-0.01, 0.03]	-0.006 [-0.03, 0.02]
EAS	-0.005 [-0.02, 0.01]	-0.01 [-0.03, 0.007]	-0.03 [-0.05, -0.009]	0.02 [0.006, 0.04]	-0.007 [-0.03, 0.01]
GSOEP	-0.000 [-0.004, 0.003]	-0.002 [-0.005, 0.002]	0.004 [-0.000, 0.008]	0.001 [-0.002, 0.004]	-0.001 [-0.004, 0.002]
HILDA	0.001 [0.000, 0.003]	0.000 [-0.001, 0.002]	0.003 [0.002, 0.005]	-0.003 [-0.004, -0.001]	-0.002 [-0.003, -0.000]
HRS	-0.002 [-0.005, 0.001]	-0.001 [-0.005, 0.002]	-0.006 [-0.010, -0.003]	0.002 [-0.001, 0.004]	-0.005 [-0.008, -0.002]
LASA				-0.000 [-0.003, 0.002]
MAP	-0.002 [-0.009, 0.005]		0.003 [-0.006, 0.01]	-0.005 [-0.01, 0.002]
MARS				-0.01 [-0.03, 0.009]
OCTO-Twin	-0.01 [-0.06, 0.03]			0.03 [-0.01, 0.07]
ROS	0.004 [-0.006, 0.01]	0.01 [0.001, 0.03]	0.002 [-0.01, 0.01]	-0.01 [-0.02, -0.001]	0.02 [0.009, 0.04]
SATSA	-0.007 [-0.02, 0.001]	-0.005 [-0.01, 0.005]	0.002 [-0.007, 0.01]	0.009 [0.001, 0.02]	0.006 [-0.007, 0.02]

Figure 6.2: Figure S9. Prospective sample-specific and overall associations between extraversion (in POMP units, 0-10) and crystallized abilities (in POMP units, 0-10) across age (M = 60, +/- 1 SD in each sample). Different colors and line types indicate different samples.

For education, 12 of the 40 tested moderators were significant. Of these, the traits where education most consistently moderated the association with crystallized abilities included extraversion (3/9) and agreeableness (3/6). For agreeableness (see Figure S10), all of the moderator associations were in the same direction (-) suggested that more educated people had weaker associations than less educated people in those samples. In contrast, for extraversion, education negatively moderated the association in HILDA and GSOEP (i.e. the association was weaker for more educated people relative to less educated people) but positive for HRS (i.e. the association was less negative for more educated people relative to less educated people).

Figure 6.3: Figure S10. Prospective sample-specific and overall associations between extraversion (in POMP units, 0-10) and crystallized abilities (in POMP units, 0-10) across education (M = 12 years, +/- 1 SD in each sample). Different colors and line types indicate different samples.

rm(list = ls()[grepl("ipd4", ls())])