Sensitivity analysis: control parcel definition
Dani Cosme
2022-02-04
This document provides the code for the analyses reported in:
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.
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
domain x age
# 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")
}domain x target x age
# 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
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") %>%
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"))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.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")))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(-.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)))
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, .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
domain x age
# 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")
}domain x target x age
# 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
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") %>%
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"))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.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")))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(-.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)))
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, .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))