• load packages
  • define color palettes and labels
  • prep data
  • sensorimotor regions
    • run models
    • compare models
    • summarize best fitting model
    • simple slopes
    • visualize fitted models
  • frontoparietal regions
    • run models
    • compare models
    • summarize best fitting model
    • simple slopes
    • visualize fitted models

This document provides the code for the analyses reported in:

Cosme et al. (Preprint) Testing the adolescent social reorientation model using hierarchical growth curve modeling with parcellated fMRI data

This script reproduces the analyses reported in the main manuscript using different sets of parcels categorized as control regions. Specifically, parcels from the original set of control regions that overlapped with the Yeo et al. (2011) 7-network sensorimotor or frontoparietal maps were categorized as control regions.

load packages

library(tidyverse)
library(knitr)
library(lme4)
library(lmerTest)

define color palettes and labels

pal_self_other = c("#FFA90A", "#247BA0")
pal_social_academic = c("#63647E", "#F25F5C")
pal_wave = c("#693668", "#A74482", "#F84AA7")
pal_label = c("#47A8BD", "#DBC057", "#FF3366")
pal_gender = c("#70c1b3","#247BA0")

parcel_labeller = labeller(label = c('social' = 'social parcels', 'other' = 'control parcels', 'self' = 'self parcels'),
                           domain = c('social' = 'social domain', 'academic' = 'academic domain'),
                           wave = c("t1" = "wave 1", "t2" = "wave 2", "t3" = "wave 3"))

label_df = expand.grid(label = c("social", "self", "other"),
              target = c("self", "other"),
              domain = c("social", "academic"),
              age = 13,
              expected_avg = 1,
              expected_diff = 1)

dcbw = theme_classic() +
  theme(text = element_text(size = 14, family = "Futura Medium", color = "black"),
        panel.background = element_blank(),
        plot.background = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 14),
        legend.background = element_rect(fill = NA, color = NA),
        axis.line = element_line(color = "black"),
        axis.text = element_text(color = "black"),
        panel.grid.minor = element_blank())

prep data

# demographics
demo = read.csv("../data/SFIC_age.pds.gender.csv") %>%
  rename("subjectID" = SID,
         "wave" = wavenum,
         "gender" = Gender) %>%
  mutate(subjectID = sprintf("s%03d", subjectID),
         wave = paste0("t", wave),
         age_c = age - 13,
         age_c2 = age_c^2,
         pdss_c = pdss - 3,
         pdss_c2 = pdss_c^2)

# load MRI data and merge
# define parcels
self_parcels = c(5, 17, 28, 47, 52, 66, 116, 147, 152, 156, 184, 198, 207, 225, 249, 292, 309, 354, 380)
social_parcels = c(18, 23, 29, 49, 54, 59, 62, 63, 67, 76, 111, 114, 139, 143, 146, 150, 178, 179, 189, 203, 212, 224, 229, 236, 238, 239, 245, 250, 259, 266, 271, 301, 305, 310, 322, 324, 328, 331, 333, 342, 343, 350, 374, 391)
control_yeo_sm_parcels = c(1, 8, 12, 14, 21, 24, 35, 37, 42, 51, 53, 64, 75, 79, 82, 87, 88, 89, 91, 106, 107, 112, 120, 121, 129, 131, 132, 137, 140, 148, 151, 158, 164, 166, 168, 170, 182, 183, 188, 191, 193, 204, 210, 222, 226, 228, 231, 242, 243, 256, 260, 264, 265, 280, 285, 287, 288, 291, 299, 320, 325, 327, 338, 339, 340, 341, 353, 357, 365, 366, 369, 370, 372, 376, 386, 387, 390)
control_yeo_fp_parcels = c(3, 4, 9, 14, 19, 22, 26, 30, 34, 36, 41, 42, 45, 50, 55, 60, 61, 68, 70, 72, 80, 81, 83, 84, 86, 92, 93, 94, 95, 97, 105, 108, 118, 122, 123, 124, 126, 128, 130, 133, 135, 137, 141, 148, 153, 154, 155, 161, 165, 169, 175, 176, 186, 188, 192, 197, 199, 200, 201, 205, 208, 211, 213, 214, 220, 231, 234, 237, 242, 244, 246, 253, 254, 257, 262, 267, 268, 273, 286, 289, 299, 306, 308, 311, 313, 314, 315, 317, 318, 319, 325, 332, 336, 341, 342, 346, 353, 357, 366, 369, 383, 384, 387, 388, 393, 395, 397)

# mri exclusions
mri_exclusions = c('s002_t1', 's004_t1', 's008_t1', 's011_t1', 's017_t1', 
                   's026_t1', 's033_t2', 's034_t1', 's041_t1', 's044_t1', 
                   's047_t1', 's051_t1', 's054_t1', 's057_t1', 's059_t1', 
                   's061_t1', 's063_t1', 's070_t2', 's074_t1', 's074_t2', 
                   's078_t1', 's084_t1', 's090_t2', 's090_t3', 's094_t1', 
                   's094_t2', 's096_t1') 

# load and tidy parcel data
parcellations_sm = read_csv('../data/fxParcellations.csv') %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 
                 ifelse(parcellation %in% control_yeo_sm_parcels, 'other', 'remove')))) %>%
  filter(!label == "remove") %>%
  mutate(wave = paste0("t", as.numeric(c(`10` = 1, `13` = 2, `16` = 3)[as.character(age)]))) %>%
  select(-age) %>%
  unite(sub_wave, c(subjectID, wave), remove = FALSE) %>%
  group_by(parcellation) %>%
  mutate(inclusion = ifelse(sub_wave %in% mri_exclusions, "excluded from MRI", "completed MRI"),
         beta = ifelse(sub_wave %in% mri_exclusions, NA, beta),
         sd = ifelse(sub_wave %in% mri_exclusions, NA, sd),
         beta_std = scale(beta, center = FALSE, scale = TRUE),
         mean_beta_std = mean(beta_std, na.rm = TRUE)) %>%
    select(-sub_wave) %>%
  ungroup()

parcellations_fp = read_csv('../data/fxParcellations.csv') %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 
                 ifelse(parcellation %in% control_yeo_fp_parcels, 'other', 'remove')))) %>%
  filter(!label == "remove") %>%
  mutate(wave = paste0("t", as.numeric(c(`10` = 1, `13` = 2, `16` = 3)[as.character(age)]))) %>%
  select(-age) %>%
  unite(sub_wave, c(subjectID, wave), remove = FALSE) %>%
  group_by(parcellation) %>%
  mutate(inclusion = ifelse(sub_wave %in% mri_exclusions, "excluded from MRI", "completed MRI"),
         beta = ifelse(sub_wave %in% mri_exclusions, NA, beta),
         sd = ifelse(sub_wave %in% mri_exclusions, NA, sd),
         beta_std = scale(beta, center = FALSE, scale = TRUE),
         mean_beta_std = mean(beta_std, na.rm = TRUE)) %>%
    select(-sub_wave) %>%
  ungroup()

# exclude parameter estimates 3 SD from the mean
parcellations_ex_sm = parcellations_sm %>%
  mutate(beta_std = ifelse(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3, NA, beta_std))

parcellations_ex_fp = parcellations_fp %>%
  mutate(beta_std = ifelse(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3, NA, beta_std))

# merge data
merged_sm = parcellations_sm %>%
  full_join(., demo, by = c("subjectID", "wave")) %>%
  mutate(inclusion = ifelse(is.na(inclusion), "didn't complete MRI", inclusion)) %>%
  filter(!(subjectID == "s086" & wave == "t3")) %>% #no MRI, task, or self-report data was collected
  filter(!is.na(beta)) %>%
  select(subjectID, wave, age, age_c, age_c2, target, domain, parcellation, label, beta, beta_std) %>%
  na.omit() %>%
  mutate(target = ifelse(target == "self", .5, -.5),
         domain = ifelse(domain == "social", .5, -.5))

merged_fp = parcellations_fp %>%
  full_join(., demo, by = c("subjectID", "wave")) %>%
  mutate(inclusion = ifelse(is.na(inclusion), "didn't complete MRI", inclusion)) %>%
  filter(!(subjectID == "s086" & wave == "t3")) %>% #no MRI, task, or self-report data was collected
  filter(!is.na(beta)) %>%
  select(subjectID, wave, age, age_c, age_c2, target, domain, parcellation, label, beta, beta_std) %>%
  na.omit() %>%
  mutate(target = ifelse(target == "self", .5, -.5),
         domain = ifelse(domain == "social", .5, -.5))

sensorimotor regions

run models

# specify model
model_domain_equation = formula(beta_std ~ 1 + age_c*domain*label +
                                  age_c2*domain*label +
                                  (1 + age_c*domain | subjectID) +
                                  (1 + age_c*domain + age_c2*domain | parcellation))

model_domain_formula = lFormula(model_domain_equation, data = merged_sm)

# calculate max number of iterations
model_domain_numFx = length(dimnames(model_domain_formula$X)[[2]])

model_domain_numRx = sum(as.numeric(lapply(model_domain_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_domain_maxfun = 10*(model_domain_numFx + model_domain_numRx + 1)^2

# run or load the model
if (file.exists("../data/model_domain_yeo_sm.RDS")) {
  model_domain = readRDS("../data/model_domain_yeo_sm.RDS")
} else {
  model_domain = lmer(model_domain_equation, data = merged_sm, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_domain, "../data/model_domain_yeo_sm.RDS")
}
# specify model
model_target_equation = formula(beta_std ~ 1 + age_c*target*domain*label +
                                  age_c2*target*domain*label +
                                  (1 + age_c*target*domain | subjectID) +
                                  (1 + age_c*target*domain + age_c2*target*domain | parcellation))

# calculate max number of iterations
model_target_formula = lFormula(model_target_equation, data = merged_sm)

model_target_numFx = length(dimnames(model_target_formula$X)[[2]])

model_target_numRx = sum(as.numeric(lapply(model_target_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_target_maxfun = 10*(model_target_numFx + model_target_numRx + 1)^2

# run or load the model
if (file.exists("../data/model_target_yeo_sm.RDS")) {
  model_target = readRDS("../data/model_target_yeo_sm.RDS")
} else {
  model_target = lmer(model_target_equation, data = merged_sm, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_target, "../data/model_target_yeo_sm.RDS")
}

compare models

anova(model_domain, model_target) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
model_domain 50 179086.6 – – –
model_target 151 178695.2 593.34 101 < .001

summarize best fitting model

model_target %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept (age 13, label (control))", term),
         term = gsub("target", "Target", term),
         term = gsub("domain", "Domain", term),
         term = gsub("labelself", "Label (self)", term),
         term = gsub("labelsocial", "Label (social)", term),
         term = gsub("age_c", "Age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.3f [%.3f, %.3f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 3) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
Intercept (age 13, label (control)) -0.089 [-0.207, 0.029] 0.060 -1.488 154.354 .139
Age -0.009 [-0.023, 0.005] 0.007 -1.274 110.500 .205
Target 0.085 [0.051, 0.118] 0.017 4.963 185.756 < .001
Domain -0.030 [-0.074, 0.014] 0.022 -1.366 160.579 .174
Label (self) 0.054 [-0.205, 0.313] 0.131 0.411 140.035 .682
Label (social) 0.330 [0.139, 0.521] 0.097 3.411 140.021 .001
Age2 0.006 [0.004, 0.008] 0.001 5.174 275.060 < .001
Age x Target 0.008 [-0.001, 0.018] 0.005 1.679 83.733 .097
Age x Domain -0.016 [-0.030, -0.002] 0.007 -2.270 73.731 .026
Target x Domain -0.043 [-0.107, 0.022] 0.033 -1.305 230.162 .193
Age x Label (self) 0.020 [0.002, 0.039] 0.009 2.193 139.855 .030
Age x Label (social) 0.032 [0.018, 0.046] 0.007 4.637 139.721 < .001
Target x Label (self) 0.180 [0.118, 0.242] 0.032 5.694 304.359 < .001
Target x Label (social) -0.126 [-0.172, -0.080] 0.023 -5.406 302.249 < .001
Domain x Label (self) 0.095 [0.021, 0.169] 0.038 2.531 186.386 .012
Domain x Label (social) 0.147 [0.093, 0.202] 0.028 5.321 185.475 < .001
Target x Age2 -0.010 [-0.013, -0.006] 0.002 -5.439 994.816 < .001
Domain x Age2 -0.004 [-0.008, -0.001] 0.002 -2.297 1727.074 .022
Label (self) x Age2 -0.004 [-0.008, 0.000] 0.002 -1.853 145.026 .066
Label (social) x Age2 -0.010 [-0.013, -0.006] 0.002 -6.028 144.205 < .001
Age x Target x Domain -0.001 [-0.021, 0.020] 0.010 -0.058 72.710 .954
Age x Target x Label (self) 0.021 [0.007, 0.035] 0.007 2.972 1628.706 .003
Age x Target x Label (social) 0.001 [-0.009, 0.011] 0.005 0.245 1617.949 .806
Age x Domain x Label (self) 0.006 [-0.009, 0.021] 0.008 0.759 375.885 .448
Age x Domain x Label (social) 0.007 [-0.004, 0.019] 0.006 1.323 373.762 .187
Target x Domain x Label (self) 0.020 [-0.098, 0.137] 0.060 0.327 1108.322 .743
Target x Domain x Label (social) 0.140 [0.053, 0.226] 0.044 3.170 1099.879 .002
Target x Domain x Age2 0.003 [-0.004, 0.010] 0.004 0.794 2534.928 .427
Target x Label (self) x Age2 -0.003 [-0.011, 0.004] 0.004 -0.955 1003.759 .340
Target x Label (social) x Age2 0.012 [0.006, 0.017] 0.003 4.314 995.825 < .001
Domain x Label (self) x Age2 0.010 [0.003, 0.017] 0.004 2.828 1581.638 .005
Domain x Label (social) x Age2 0.009 [0.004, 0.014] 0.003 3.482 1568.709 .001
Age x Target x Domain x Label (self) -0.001 [-0.029, 0.026] 0.014 -0.082 2761.297 .935
Age x Target x Domain x Label (social) -0.009 [-0.029, 0.011] 0.010 -0.866 2742.990 .387
Target x Domain x Label (self) x Age2 0.004 [-0.010, 0.017] 0.007 0.507 8601.358 .612
Target x Domain x Label (social) x Age2 -0.006 [-0.016, 0.004] 0.005 -1.199 8535.720 .231 
# summarize random effects
print(VarCorr(model_target), comp = c("Variance","Std.Dev."), digits = 3)
##  Groups       Name                 Variance   Std.Dev. Corr                   
##  parcellation (Intercept)          0.25816928 0.50810                         
##               age_c                0.00114872 0.03389   0.55                  
##               target               0.00229747 0.04793  -0.35 -0.23            
##               domain               0.00863484 0.09292   0.07 -0.07 -0.46      
##               age_c2               0.00002430 0.00493  -0.50 -0.51  0.39 -0.43
##               age_c:target         0.00002846 0.00533  -0.10  0.63  0.45 -0.11
##               age_c:domain         0.00017380 0.01318   0.78  0.84  0.02 -0.32
##               target:domain        0.00288012 0.05367   0.03  0.23 -0.09  0.82
##               target:age_c2        0.00001420 0.00377   0.94  0.61 -0.55  0.35
##               domain:age_c2        0.00000568 0.00238   0.54  0.40  0.21  0.21
##               age_c:target:domain  0.00006673 0.00817  -0.25 -0.31  0.54 -0.92
##               target:domain:age_c2 0.00000504 0.00225  -0.90 -0.56  0.15 -0.15
##  subjectID    (Intercept)          0.00961680 0.09807                         
##               age_c                0.00164330 0.04054   0.30                  
##               target               0.00500456 0.07074  -0.19  0.29            
##               domain               0.01167156 0.10804   0.08  0.21  0.14      
##               age_c:target         0.00071070 0.02666   0.04 -0.19 -0.13 -0.41
##               age_c:domain         0.00215177 0.04639  -0.13  0.26  0.30 -0.12
##               target:domain        0.01877223 0.13701   0.14 -0.10  0.29 -0.28
##               age_c:target:domain  0.00370953 0.06091   0.14 -0.02 -0.09  0.02
##  Residual                          0.61199067 0.78230                         
##                                           
##                                           
##                                           
##                                           
##                                           
##                                           
##  -0.27                                    
##  -0.45  0.46                              
##  -0.62  0.45 -0.02                        
##  -0.69 -0.04  0.68  0.30                  
##  -0.76  0.46  0.60  0.57  0.57            
##   0.64 -0.14  0.00 -0.87 -0.55 -0.36      
##   0.73 -0.17 -0.79 -0.31 -0.89 -0.85  0.35
##                                           
##                                           
##                                           
##                                           
##                                           
##   0.23                                    
##   0.55  0.45                              
##   0.30 -0.10  0.04                        
##

simple slopes

Estimate simple slopes to test interactions at specific levels

# self social > academic
self_social = emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c2", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "quadratic"))

# social self > other
social_self = emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c2", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "quadratic"))
social_self %>%
  bind_rows(self_social) %>%
  select(contrast, parcel, age_effect, estimate, SE, df, t.ratio, p.value) %>%
  rename("b" = estimate,
         "t" = t.ratio,
         "p" = p.value) %>%
  mutate(b = round(b, 3) * -1, #flip signs for it's .5 - (-.5)
         SE = round(SE, 3),
         df = round(df, 2),
         t = abs(round(t, 2)),
         p = round(p, 3)) %>%
  kable(format = "pandoc")
contrast parcel age_effect b SE df t p
social self > other self linear 0.028 0.011 275.32 2.57 0.011
social self > other self quadratic -0.010 0.005 3314.59 2.14 0.033
self social > academic social linear -0.013 0.009 88.56 1.42 0.159
self social > academic social quadratic 0.003 0.003 4486.70 0.99 0.324

visualize fitted models

Visualize the developmental trajectory using the fitted values from the domain x target x age model

reForm = as.formula("~(1 + age_c*target*domain + age_c2*target*domain | parcellation)")

neuro_plot_data = with(merged_sm, 
                    expand.grid(target = unique(target), 
                                domain = unique(domain),
                                parcellation = unique(parcellation),
                                age = unique(age),
                                stringsAsFactors = F)) %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other')),
         age_c = age - 13, 
         age_c2 = age_c^2, 
         subjectID = NA)
neuro_plot_data$expected = predict(model_target, newdata = neuro_plot_data, re.form = reForm)
neuro_plot_data$expected_mean = predict(model_target, newdata = neuro_plot_data, re.form = NA)

neuro_plot_data = neuro_plot_data %>%
  mutate(target = factor(target, levels = c(-.5, .5), labels = c("other", "self")),
         domain = factor(domain, levels = c(-.5, .5), labels = c("academic", "social")))
domain_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(name = "", values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n", color = "") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_fitted = cowplot::plot_grid(domain_parc_plot, target_parc_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

domain_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), color = NA, alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.3, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_social_academic[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.3, .45)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h1_fitted = cowplot::plot_grid(domain_plot, soc_acad_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

target_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, target, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), color = NA, alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .35)) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_self_other[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .35)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h2_fitted = cowplot::plot_grid(target_plot, self_other_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

int_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_mean, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = 1, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .55)) +
  scale_color_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '', linetype = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .55)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))


self_other_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .55)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_fitted = cowplot::plot_grid(int_plot, soc_acad_int_plot, self_other_int_plot,
                       labels = c('A', 'B', 'C'), ncol = 3))

frontoparietal regions

run models

Estimate full models

# specify model
model_domain_equation = formula(beta_std ~ 1 + age_c*domain*label +
                                  age_c2*domain*label +
                                  (1 + age_c*domain | subjectID) +
                                  (1 + age_c*domain + age_c2*domain | parcellation))

model_domain_formula = lFormula(model_domain_equation, data = merged_fp)

# calculate max number of iterations
model_domain_numFx = length(dimnames(model_domain_formula$X)[[2]])

model_domain_numRx = sum(as.numeric(lapply(model_domain_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_domain_maxfun = 10*(model_domain_numFx + model_domain_numRx + 1)^2

# run or load the model
if (file.exists("../data/model_domain_yeo_fp.RDS")) {
  model_domain = readRDS("../data/model_domain_yeo_fp.RDS")
} else {
  model_domain = lmer(model_domain_equation, data = merged_fp, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_domain, "../data/model_domain_yeo_fp.RDS")
}
# specify model
model_target_equation = formula(beta_std ~ 1 + age_c*target*domain*label +
                                  age_c2*target*domain*label +
                                  (1 + age_c*target*domain | subjectID) +
                                  (1 + age_c*target*domain + age_c2*target*domain | parcellation))

# calculate max number of iterations
model_target_formula = lFormula(model_target_equation, data = merged_fp)

model_target_numFx = length(dimnames(model_target_formula$X)[[2]])

model_target_numRx = sum(as.numeric(lapply(model_target_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_target_maxfun = 10*(model_target_numFx + model_target_numRx + 1)^2

# run or load the model
if (file.exists("../data/model_target_yeo_fp.RDS")) {
  model_target = readRDS("../data/model_target_yeo_fp.RDS")
} else {
  model_target = lmer(model_target_equation, data = merged_fp, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_target, "../data/model_target_yeo_fp.RDS")
}

compare models

anova(model_domain, model_target) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
model_domain 50 221943.0 – – –
model_target 151 221323.4 821.64 101 < .001

summarize best fitting model

model_target %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept (age 13, label (control))", term),
         term = gsub("target", "Target", term),
         term = gsub("domain", "Domain", term),
         term = gsub("labelself", "Label (self)", term),
         term = gsub("labelsocial", "Label (social)", term),
         term = gsub("age_c", "Age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.3f [%.3f, %.3f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 3) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
Intercept (age 13, label (control)) 0.006 [-0.083, 0.094] 0.045 0.128 198.660 .898
Age 0.020 [0.010, 0.030] 0.005 3.893 92.502 < .001
Target -0.006 [-0.042, 0.031] 0.018 -0.305 161.983 .761
Domain 0.064 [0.007, 0.122] 0.029 2.223 110.094 .028
Label (self) -0.032 [-0.249, 0.185] 0.110 -0.290 169.041 .772
Label (social) 0.243 [0.087, 0.399] 0.079 3.076 169.018 .002
Age2 -0.000 [-0.002, 0.002] 0.001 -0.308 329.568 .759
Age x Target 0.001 [-0.010, 0.012] 0.005 0.152 60.362 .880
Age x Domain -0.010 [-0.026, 0.006] 0.008 -1.259 62.238 .213
Target x Domain 0.059 [-0.007, 0.124] 0.033 1.770 142.093 .079
Age x Label (self) -0.009 [-0.023, 0.004] 0.007 -1.385 169.146 .168
Age x Label (social) 0.002 [-0.008, 0.011] 0.005 0.362 168.846 .718
Target x Label (self) 0.275 [0.207, 0.342] 0.034 8.029 235.979 < .001
Target x Label (social) -0.032 [-0.080, 0.017] 0.025 -1.284 234.664 .200
Domain x Label (self) 0.001 [-0.079, 0.082] 0.041 0.034 202.414 .973
Domain x Label (social) 0.052 [-0.006, 0.110] 0.029 1.761 201.624 .080
Target x Age2 -0.006 [-0.009, -0.003] 0.002 -3.845 1268.643 < .001
Domain x Age2 -0.004 [-0.008, -0.001] 0.002 -2.414 2874.494 .016
Label (self) x Age2 0.000 [-0.004, 0.004] 0.002 0.030 174.605 .976
Label (social) x Age2 -0.005 [-0.009, -0.002] 0.002 -3.492 173.663 .001
Age x Target x Domain -0.000 [-0.022, 0.022] 0.011 -0.003 60.729 .997
Age x Target x Label (self) 0.025 [0.011, 0.039] 0.007 3.548 2188.276 < .001
Age x Target x Label (social) 0.006 [-0.004, 0.016] 0.005 1.127 2174.365 .260
Age x Domain x Label (self) -0.003 [-0.018, 0.011] 0.007 -0.442 783.921 .659
Age x Domain x Label (social) -0.002 [-0.013, 0.008] 0.005 -0.470 779.237 .638
Target x Domain x Label (self) -0.074 [-0.191, 0.042] 0.059 -1.247 2396.209 .213
Target x Domain x Label (social) 0.047 [-0.037, 0.131] 0.043 1.102 2378.915 .271
Target x Domain x Age2 -0.006 [-0.013, 0.000] 0.003 -1.947 1898.791 .052
Target x Label (self) x Age2 -0.009 [-0.016, -0.002] 0.004 -2.594 1819.207 .010
Target x Label (social) x Age2 0.006 [0.001, 0.011] 0.003 2.165 1804.803 .031
Domain x Label (self) x Age2 0.009 [0.002, 0.016] 0.004 2.480 4710.718 .013
Domain x Label (social) x Age2 0.008 [0.003, 0.013] 0.003 3.175 4673.747 .002
Age x Target x Domain x Label (self) -0.003 [-0.030, 0.025] 0.014 -0.186 3924.226 .853
Age x Target x Domain x Label (social) -0.010 [-0.030, 0.010] 0.010 -0.979 3899.379 .328
Target x Domain x Label (self) x Age2 0.011 [-0.003, 0.025] 0.007 1.525 8776.891 .127
Target x Domain x Label (social) x Age2 0.001 [-0.009, 0.011] 0.005 0.108 8711.047 .914 
# summarize random effects
print(VarCorr(model_target), comp = c("Variance","Std.Dev."), digits = 3)
##  Groups       Name                 Variance   Std.Dev. Corr                   
##  parcellation (Intercept)          0.19097710 0.43701                         
##               age_c                0.00054244 0.02329   0.37                  
##               target               0.00508366 0.07130   0.00 -0.08            
##               domain               0.01303988 0.11419   0.22  0.07 -0.26      
##               age_c2               0.00002621 0.00512  -0.38 -0.35 -0.02 -0.60
##               age_c:target         0.00002627 0.00513  -0.40  0.14  0.80 -0.49
##               age_c:domain         0.00009194 0.00959   0.79  0.62  0.11 -0.30
##               target:domain        0.00169283 0.04114  -0.09  0.64 -0.32  0.57
##               target:age_c2        0.00000677 0.00260   0.20  0.00 -0.95  0.20
##               domain:age_c2        0.00000188 0.00137  -0.14 -0.20  0.19 -0.59
##               age_c:target:domain  0.00006148 0.00784   0.62  0.05  0.56 -0.47
##               target:domain:age_c2 0.00000521 0.00228  -0.62 -0.42  0.33  0.29
##  subjectID    (Intercept)          0.01115553 0.10562                         
##               age_c                0.00096577 0.03108   0.18                  
##               target               0.00866331 0.09308   0.22  0.28            
##               domain               0.03659842 0.19131   0.14  0.25  0.22      
##               age_c:target         0.00117666 0.03430  -0.24 -0.59 -0.17 -0.09
##               age_c:domain         0.00302648 0.05501  -0.30  0.18 -0.08 -0.41
##               target:domain        0.03169049 0.17802  -0.29 -0.22 -0.29 -0.43
##               age_c:target:domain  0.00522723 0.07230  -0.30 -0.40  0.05  0.10
##  Residual                          0.65432972 0.80891                         
##                                           
##                                           
##                                           
##                                           
##                                           
##                                           
##   0.04                                    
##  -0.02 -0.02                              
##  -0.64 -0.02 -0.16                        
##  -0.06 -0.85  0.01  0.16                  
##  -0.16  0.39 -0.01 -0.23 -0.02            
##   0.19  0.21  0.76 -0.66 -0.36  0.18      
##   0.25  0.27 -0.70  0.02 -0.53 -0.45 -0.38
##                                           
##                                           
##                                           
##                                           
##                                           
##  -0.05                                    
##   0.54  0.43                              
##   0.15 -0.23 -0.15                        
##

simple slopes

Estimate simple slopes to test interactions at specific levels

# self social > academic
self_social = emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c2", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "quadratic"))

# social self > other
social_self = emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c2", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "quadratic"))
social_self %>%
  bind_rows(self_social) %>%
  select(contrast, parcel, age_effect, estimate, SE, df, t.ratio, p.value) %>%
  rename("b" = estimate,
         "t" = t.ratio,
         "p" = p.value) %>%
  mutate(b = round(b, 3) * -1, #flip signs for it's .5 - (-.5)
         SE = round(SE, 3),
         df = round(df, 2),
         t = abs(round(t, 2)),
         p = round(p, 3)) %>%
  kable(format = "pandoc")
contrast parcel age_effect b SE df t p
social self > other self linear 0.025 0.012 208.05 2.09 0.037
social self > other self quadratic -0.013 0.005 4848.62 2.78 0.005
self social > academic social linear -0.018 0.010 103.02 1.73 0.087
self social > academic social quadratic 0.001 0.003 6762.88 0.32 0.749

visualize fitted models

Visualize the developmental trajectory using the fitted values from the domain x target x age model

reForm = as.formula("~(1 + age_c*target*domain + age_c2*target*domain | parcellation)")

neuro_plot_data = with(merged_fp, 
                    expand.grid(target = unique(target), 
                                domain = unique(domain),
                                parcellation = unique(parcellation),
                                age = unique(age),
                                stringsAsFactors = F)) %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other')),
         age_c = age - 13, 
         age_c2 = age_c^2, 
         subjectID = NA)
neuro_plot_data$expected = predict(model_target, newdata = neuro_plot_data, re.form = reForm)
neuro_plot_data$expected_mean = predict(model_target, newdata = neuro_plot_data, re.form = NA)

neuro_plot_data = neuro_plot_data %>%
  mutate(target = factor(target, levels = c(-.5, .5), labels = c("other", "self")),
         domain = factor(domain, levels = c(-.5, .5), labels = c("academic", "social")))
domain_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(name = "", values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n", color = "") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_fitted = cowplot::plot_grid(domain_parc_plot, target_parc_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

domain_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), color = NA, alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.3, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_social_academic[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.3, .45)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h1_fitted = cowplot::plot_grid(domain_plot, soc_acad_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

target_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, target, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), color = NA, alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .35)) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_self_other[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .35)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h2_fitted = cowplot::plot_grid(target_plot, self_other_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

int_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_mean, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = 1, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .55)) +
  scale_color_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '', linetype = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .55)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))


self_other_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .55)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_fitted = cowplot::plot_grid(int_plot, soc_acad_int_plot, self_other_int_plot,
                       labels = c('A', 'B', 'C'), ncol = 3))