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 but estimating only linear relationships with age (i.e., removing the quadratic terms).

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

load MRI data

  • exclude participants with >20% volumes with motion artifacts or low quality FX data
  • standardize betas within parcel (across participants and time)
  • recode standardized betas +/- 3 SDs from the mean as NA (0.8% of all observations)
# 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)

# 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 = read_csv('../data/fxParcellations.csv') %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other'))) %>%
  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 = parcellations %>%
  mutate(beta_std = ifelse(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3, NA, beta_std))

load demographic data and merge

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

# merge data
merged = parcellations_ex %>%
  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

# subset data for modeling
neuro_model_data = merged %>%
  filter(!is.na(beta)) %>%
  select(subjectID, wave, age, age_c, age_c2, target, domain, parcellation, label, beta, beta_std)

neuro_model_data_ex = neuro_model_data  %>%
  na.omit()

# dummy code target and domain
neuro_model_data_dummy = neuro_model_data_ex %>%
  mutate(target = ifelse(target == "self", .5, -.5),
         domain = ifelse(domain == "social", .5, -.5))

run domain x target x age model

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

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

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_linear.RDS")) {
  model_target = readRDS("../data/model_target_linear.RDS")
} else {
  model_target = lmer(model_target_equation, data = neuro_model_data_dummy, 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_linear.RDS")
}

summarize 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.003 [-0.055, 0.050] 0.027 -0.097 381.610 .923
Age 0.003 [-0.002, 0.009] 0.003 1.201 144.269 .232
Target -0.006 [-0.020, 0.008] 0.007 -0.884 73.600 .379
Domain -0.004 [-0.027, 0.018] 0.011 -0.385 92.478 .701
Label (self) -0.025 [-0.232, 0.182] 0.105 -0.237 352.012 .813
Label (social) 0.214 [0.072, 0.355] 0.072 2.972 351.997 .003
Age x Target 0.002 [-0.003, 0.007] 0.002 0.941 69.740 .350
Age x Domain 0.000 [-0.007, 0.008] 0.004 0.068 59.451 .946
Target x Domain -0.007 [-0.036, 0.022] 0.015 -0.499 59.151 .619
Age x Label (self) 0.007 [-0.008, 0.022] 0.008 0.961 352.388 .337
Age x Label (social) 0.018 [0.007, 0.028] 0.005 3.344 351.943 .001
Target x Label (self) 0.184 [0.146, 0.223] 0.019 9.485 356.977 < .001
Target x Label (social) -0.023 [-0.049, 0.003] 0.013 -1.710 355.121 .088
Domain x Label (self) 0.109 [0.055, 0.162] 0.027 4.014 358.719 < .001
Domain x Label (social) 0.154 [0.118, 0.191] 0.018 8.338 357.759 < .001
Age x Target x Domain 0.002 [-0.008, 0.013] 0.005 0.439 54.551 .663
Age x Target x Label (self) 0.024 [0.011, 0.037] 0.007 3.543 6460.934 < .001
Age x Target x Label (social) 0.007 [-0.002, 0.016] 0.005 1.582 6419.379 .114
Age x Domain x Label (self) -0.005 [-0.018, 0.009] 0.007 -0.666 1519.658 .506
Age x Domain x Label (social) -0.003 [-0.013, 0.006] 0.005 -0.716 1510.489 .474
Target x Domain x Label (self) 0.034 [-0.035, 0.103] 0.035 0.967 1982.783 .334
Target x Domain x Label (social) 0.088 [0.041, 0.136] 0.024 3.685 1970.327 < .001
Age x Target x Domain x Label (self) -0.004 [-0.030, 0.022] 0.013 -0.280 18154.597 .779
Age x Target x Domain x Label (social) -0.014 [-0.032, 0.004] 0.009 -1.535 18042.926 .125

simple slopes

Estimate simple slopes to test interactions at specific levels

run models

# 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")

# 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")

make table

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.010 1675.62 2.61 0.009
self social > academic social linear -0.009 0.007 828.95 1.34 0.179

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 | parcellation)")

neuro_plot_data = with(neuro_model_data_dummy, 
                    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")))

hypothesis 0

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

hypothesis 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(-.25, .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(-.25, .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)))

hypothesis 2

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

hypothesis 3

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, .5)) +
  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, .5)) +
  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, .5)) +
  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))