Main manuscript analyses
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.
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))sample descriptives
full sample
merged %>%
select(subjectID, wave, age, gender) %>%
unique() %>%
group_by(wave) %>%
summarize(n = n(),
meanAge = mean(age, na.rm = TRUE),
sdAge = sd(age, na.rm = TRUE)) %>%
kable(format = 'pandoc', digits = 2)| wave | n | meanAge | sdAge |
|---|---|---|---|
| t1 | 90 | 10.07 | 0.31 |
| t2 | 57 | 13.04 | 0.31 |
| t3 | 44 | 16.46 | 0.58 |
merged %>%
select(subjectID, wave, age, gender) %>%
unique() %>%
group_by(wave, gender) %>%
summarize(n = n()) %>%
spread(gender, n) %>%
kable(format = 'pandoc', digits = 2)| wave | F | M |
|---|---|---|
| t1 | 45 | 45 |
| t2 | 30 | 27 |
| t3 | 25 | 19 |
broken down by inclusion
age_table = merged %>%
select(subjectID, inclusion, wave, age, pdss, gender) %>%
unique() %>%
group_by(inclusion, wave) %>%
summarize(N = n(),
Mage = mean(age, na.rm = TRUE),
SDage = sd(age, na.rm = TRUE),
Mpdss = mean(pdss, na.rm = TRUE),
SDpdss = sd(pdss, na.rm = TRUE))
gender_table = merged %>%
select(subjectID, inclusion, wave, age, pds, gender) %>%
unique() %>%
group_by(inclusion, wave, gender) %>%
summarize(N = n()) %>%
spread(gender, N) %>%
rename("Female" = F,
"Male" = M)
merged %>%
filter(!inclusion == "didn't complete MRI") %>%
mutate(inclusion = ifelse(grepl("excluded", inclusion), "excluded", "included")) %>%
select(subjectID, wave, gender, age, inclusion) %>%
unique() %>%
filter(inclusion == "included") %>%
group_by(subjectID) %>%
summarize(`number of waves included` = n()) %>%
group_by(`number of waves included`) %>%
summarize(n = n()) %>%
kable(format = 'pandoc')| number of waves included | n |
|---|---|
| 1 | 35 |
| 2 | 17 |
| 3 | 22 |
age_table %>%
left_join(gender_table) %>%
extract(wave, "Wave", "t([1-3]{1})") %>%
arrange(Wave) %>%
select(Wave, inclusion, N, Mage, SDage, Mpdss, SDpdss, Female, Male) %>%
rename("MRI status" = inclusion) %>%
kable(format = 'pandoc', digits = 2)| Wave | MRI status | N | Mage | SDage | Mpdss | SDpdss | Female | Male |
|---|---|---|---|---|---|---|---|---|
| 1 | completed MRI | 57 | 10.08 | 0.32 | 2.16 | 0.74 | 33 | 24 |
| 1 | didn’t complete MRI | 12 | 10.00 | 0.25 | 2.00 | 0.77 | 4 | 8 |
| 1 | excluded from MRI | 21 | 10.07 | 0.33 | 1.83 | 0.58 | 8 | 13 |
| 2 | completed MRI | 44 | 13.06 | 0.33 | 3.57 | 1.05 | 25 | 19 |
| 2 | didn’t complete MRI | 8 | 13.07 | 0.20 | 3.69 | 0.65 | 4 | 4 |
| 2 | excluded from MRI | 5 | 12.79 | 0.28 | 2.62 | 0.85 | 1 | 4 |
| 3 | completed MRI | 34 | 16.33 | 0.46 | 4.75 | 0.39 | 19 | 15 |
| 3 | didn’t complete MRI | 9 | 17.06 | 0.60 | 4.78 | 0.67 | 6 | 3 |
| 3 | excluded from MRI | 1 | 15.72 | NA | 3.50 | NA | NA | 1 |
plot
full sample
merged %>%
select(subjectID, wave, gender, age, inclusion) %>%
group_by(subjectID) %>%
unique() %>%
mutate(n.waves = n()) %>%
ungroup() %>%
arrange(n.waves, gender, age, inclusion) %>%
mutate(order = row_number()) %>%
ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
geom_line(size = .3) +
geom_point(aes(shape = inclusion), size = 2, fill = "white") +
scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
scale_shape_manual(name = "", values = c(19, 17, 21)) +
labs(x = "\nage", y = "participant\n") +
dcbw +
theme(axis.line = element_line(colour = "black"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.75,.92),
legend.spacing.y = unit(.01, 'cm'),
legend.margin = unit(0, "cm"))
MRI only
merged %>%
filter(!inclusion == "didn't complete MRI") %>%
mutate(inclusion = ifelse(grepl("excluded", inclusion), "excluded", "included")) %>%
select(subjectID, wave, gender, age, inclusion) %>%
group_by(subjectID) %>%
unique() %>%
mutate(n.waves = n()) %>%
ungroup() %>%
arrange(n.waves, gender, age, inclusion) %>%
mutate(order = row_number()) %>%
ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
geom_line(size = .3) +
geom_point(aes(shape = inclusion), size = 2, fill = "white") +
scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
scale_shape_manual(name = "", values = c(21, 19)) +
labs(x = "\nage", y = "participant\n") +
dcbw +
theme(axis.line = element_line(colour = "black"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.75,.92),
legend.spacing.y = unit(.01, 'cm'),
legend.margin = unit(0, "cm"))
mri descriptives
# collapsed across wave
neuro_model_data_dummy %>%
mutate(target = ifelse(target == .5, "self", "other"),
domain = ifelse(domain == .5, "social", "academic")) %>%
group_by(domain, label) %>%
summarize(M = mean(beta_std, na.rm = TRUE),
SD = sd(beta_std, na.rm = TRUE)) %>%
gather(var, value, M, SD) %>%
unite(label, label, var, sep = " ") %>%
spread(label, value) %>%
kable(format = "pandoc", digits = 2)| domain | other M | other SD | self M | self SD | social M | social SD |
|---|---|---|---|---|---|---|
| academic | 0.01 | 0.92 | -0.08 | 0.91 | 0.14 | 0.92 |
| social | -0.01 | 0.95 | 0.02 | 0.92 | 0.27 | 0.93 |
neuro_model_data_dummy %>%
mutate(target = ifelse(target == .5, "self", "other"),
domain = ifelse(domain == .5, "social", "academic")) %>%
group_by(target, label) %>%
summarize(M = mean(beta_std, na.rm = TRUE),
SD = sd(beta_std, na.rm = TRUE)) %>%
gather(var, value, M, SD) %>%
unite(label, label, var, sep = " ") %>%
spread(label, value) %>%
kable(format = "pandoc", digits = 2)| target | other M | other SD | self M | self SD | social M | social SD |
|---|---|---|---|---|---|---|
| other | 0.00 | 0.94 | -0.11 | 0.91 | 0.22 | 0.93 |
| self | -0.01 | 0.94 | 0.05 | 0.92 | 0.19 | 0.93 |
# as a function of wave
neuro_model_data_dummy %>%
mutate(target = ifelse(target == .5, "self", "other"),
domain = ifelse(domain == .5, "social", "academic")) %>%
group_by(domain, label, wave) %>%
summarize(M = mean(beta_std, na.rm = TRUE),
SD = sd(beta_std, na.rm = TRUE)) %>%
gather(var, value, M, SD) %>%
unite(label, label, var, sep = " ") %>%
spread(label, value) %>%
kable(format = "pandoc", digits = 2)| domain | wave | other M | other SD | self M | self SD | social M | social SD |
|---|---|---|---|---|---|---|---|
| academic | t1 | 0.00 | 0.92 | -0.13 | 0.95 | 0.06 | 0.91 |
| academic | t2 | 0.00 | 0.96 | -0.04 | 0.92 | 0.19 | 0.95 |
| academic | t3 | 0.03 | 0.89 | -0.04 | 0.82 | 0.21 | 0.88 |
| social | t1 | -0.02 | 0.94 | 0.00 | 0.91 | 0.22 | 0.91 |
| social | t2 | 0.00 | 0.96 | 0.00 | 0.95 | 0.30 | 0.98 |
| social | t3 | -0.01 | 0.96 | 0.05 | 0.89 | 0.33 | 0.90 |
neuro_model_data_dummy %>%
mutate(target = ifelse(target == .5, "self", "other"),
domain = ifelse(domain == .5, "social", "academic")) %>%
group_by(target, label, wave) %>%
summarize(M = mean(beta_std, na.rm = TRUE),
SD = sd(beta_std, na.rm = TRUE)) %>%
gather(var, value, M, SD) %>%
unite(label, label, var, sep = " ") %>%
spread(label, value) %>%
kable(format = "pandoc", digits = 2)| target | wave | other M | other SD | self M | self SD | social M | social SD |
|---|---|---|---|---|---|---|---|
| other | t1 | 0.00 | 0.92 | -0.10 | 0.92 | 0.16 | 0.90 |
| other | t2 | 0.00 | 0.97 | -0.14 | 0.94 | 0.27 | 0.98 |
| other | t3 | 0.02 | 0.92 | -0.10 | 0.85 | 0.27 | 0.88 |
| self | t1 | -0.02 | 0.94 | -0.03 | 0.94 | 0.11 | 0.92 |
| self | t2 | 0.00 | 0.94 | 0.11 | 0.92 | 0.22 | 0.95 |
| self | t3 | 0.01 | 0.93 | 0.11 | 0.86 | 0.27 | 0.90 |
create hypothetical plot (figure 1)
hyp_data = data.frame(age = c(10, 13, 16, 10, 13, 16, 10, 13, 16),
social = c(.5, .9, 1.3, .5, 1.5, .6, .5, 1, 1),
academic = c(0, .1, .2, 0, .2, .1, 0, .2, .1),
trajectory= rep(c("linear\nincrease", "early adolescent\nspecific", "adolescent\nemergent"), each = 3)) %>%
mutate(trajectory = factor(trajectory,
levels = c("linear\nincrease", "early adolescent\nspecific", "adolescent\nemergent"))) %>%
gather(domain, value, social, academic)
domain_plot_hyp = hyp_data %>%
ggplot(aes(x = age, y = value, color = domain)) +
geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) +
facet_grid(~ trajectory) +
scale_color_manual(values = pal_social_academic) +
scale_x_continuous(breaks = c(10, 13, 16)) +
labs(x = "\nage", y = "hypothetical BOLD signal value\n", color = '') +
dcbw +
theme(legend.position = c(.12, .83),
legend.spacing.y = unit(.01, 'cm'),
legend.margin = unit(0, "cm"))
soc_acad_plot_hyp = hyp_data %>%
spread(domain, value) %>%
mutate(diff = social - academic) %>%
ggplot(aes(x = age, y = diff, color = domain)) +
geom_smooth(method = 'lm', formula = y ~ poly(x,2),
se = FALSE, color = pal_social_academic[2], size = 1.5) +
facet_grid(~ trajectory) +
scale_x_continuous(breaks = c(10, 13, 16)) +
labs(x = "\nage", y = "hypothetical BOLD signal value\n", color = '') +
dcbw +
theme(legend.position = "none")
(hyp = cowplot::plot_grid(domain_plot_hyp, soc_acad_plot_hyp,
labels = c('A', 'B'), ncol = 2,
rel_widths = c(1, 1)))
visualize raw
Visualize the developmental trajectories using the raw data
hypothesis 0
domain_parc_plot_raw = neuro_model_data %>%
distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
group_by(subjectID, wave, label, domain, parcellation) %>%
mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
ggplot(aes(x = age, y = beta_std_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(aes(fill = domain), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
scale_color_manual(values = pal_social_academic) +
scale_fill_manual(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 standardized parameter estimate\n', color = '', fill = '') +
theme_minimal(base_size = 12) +
dcbw +
theme(legend.position = c(.85, .15),
legend.spacing.y = unit(.01, 'cm'),
legend.margin = unit(0, "cm"))
target_parc_plot_raw = neuro_model_data %>%
distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
group_by(subjectID, wave, label, target, parcellation) %>%
mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
ggplot(aes(x = age, y = beta_std_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(aes(fill = target), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
scale_color_manual(values = pal_self_other) +
scale_fill_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 standardized parameter estimate\n', color = '', fill = '') +
theme_minimal(base_size = 12) +
dcbw +
theme(legend.position = c(.85, .15),
legend.spacing.y = unit(.01, 'cm'),
legend.margin = unit(0, "cm"))
(h0_raw = cowplot::plot_grid(domain_parc_plot_raw, target_parc_plot_raw,
labels = c('A', 'B'), ncol = 2,
rel_widths = c(1, 1)))
hypothesis 1
domain_plot_raw = neuro_model_data %>%
distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
group_by(subjectID, wave, label, domain, parcellation) %>%
mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
ggplot(aes(x = age, y = beta_std_avg, group = domain, color = domain)) +
geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) +
scale_color_manual(values = pal_social_academic) +
scale_fill_manual(values = pal_social_academic) +
scale_x_continuous(breaks = c(10, 13, 16)) +
scale_y_continuous(breaks = round(seq(-.4, .6, .2), 1)) +
coord_cartesian(ylim = c(-.45, .65)) +
facet_grid(~ label, labeller = parcel_labeller) +
labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
dcbw +
theme(legend.position = "none")
soc_acad_plot_raw = neuro_model_data %>%
group_by(subjectID, age, label, domain, parcellation) %>%
mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
select(subjectID, age, label, domain, beta_std_avg) %>%
unique() %>%
spread(domain, beta_std_avg) %>%
mutate(avg_diff = social - academic) %>%
distinct(parcellation, age, label, .keep_all = T) %>%
ggplot(aes(x = age, y = avg_diff)) +
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 = TRUE, color = pal_social_academic[2], fill = pal_social_academic[2], size = 1.5) +
scale_size_manual(values = c(.03, .1, .1)) +
scale_y_continuous(breaks = round(seq(-.4, .6, .2), 1)) +
coord_cartesian(ylim = c(-.45, .65)) +
scale_x_continuous(breaks = c(10, 13, 16)) +
facet_grid(~ label, labeller = parcel_labeller) +
labs(x = "\nage", y = "mean standardized BOLD signal value\n", color = '', fill = '') +
dcbw +
theme(legend.position = "none")
(h1_raw = cowplot::plot_grid(domain_plot_raw, soc_acad_plot_raw,
labels = c('A', 'B'), ncol = 2,
rel_widths = c(1, 1)))
hypothesis 2
target_plot_raw = neuro_model_data %>%
distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
group_by(subjectID, wave, label, target, parcellation) %>%
mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
ggplot(aes(x = age, y = beta_std_avg, group = target, color = target)) +
geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), alpha = .2) +
scale_color_manual(values = pal_self_other) +
scale_fill_manual(values = pal_self_other) +
scale_x_continuous(breaks = c(10, 13, 16)) +
scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
coord_cartesian(ylim = c(-.4, .5)) +
facet_grid(~ label, labeller = parcel_labeller) +
labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
dcbw +
theme(legend.position = "none")
self_other_plot_raw = neuro_model_data %>%
group_by(subjectID, age, label, target, parcellation) %>%
mutate(beta_std_avf = mean(beta_std, na.rm = TRUE)) %>%
select(subjectID, age, label, target, beta_std_avf) %>%
unique() %>%
spread(target, beta_std_avf) %>%
mutate(avg_diff = self - other) %>%
distinct(parcellation, age, label, .keep_all = T) %>%
ggplot(aes(x = age, y = avg_diff)) +
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 = TRUE, color = pal_self_other[2], fill = pal_self_other[2], size = 1.5) +
scale_size_manual(values = c(.03, .1, .1)) +
scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
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 standardized BOLD signal value\n", color = '', fill = '') +
dcbw +
theme(legend.position = "none")
(h2_raw = cowplot::plot_grid(target_plot_raw, self_other_plot_raw,
labels = c('A', 'B'), ncol = 2,
rel_widths = c(1, 1)))
hypothesis 3
int_plot_raw = neuro_model_data %>%
distinct(parcellation, target, domain, age, beta_std, label, .keep_all = T) %>%
ggplot(aes(x = age, y = beta_std, group = interaction(target, domain), color = domain, linetype = target)) +
geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) +
scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
coord_cartesian(ylim = c(-.4, .45)) +
scale_color_manual(values = pal_social_academic) +
scale_fill_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 standardized parameter estimate\n', color = '', linetype = '', fill = '') +
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_raw = neuro_model_data %>%
group_by(subjectID, age, label, target, domain, parcellation) %>%
mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
select(subjectID, age, label, target, domain, beta_std_avg) %>%
unique() %>%
spread(domain, beta_std_avg) %>%
mutate(avg_diff = social - academic) %>%
distinct(parcellation, age, label, .keep_all = T) %>%
ggplot(aes(x = age, y = avg_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(aes(fill = target), method = 'lm', formula = y ~ poly(x,2),
se = TRUE, size = 1.5) +
scale_color_manual(values = pal_self_other) +
scale_fill_manual(values = pal_self_other) +
scale_size_manual(values = c(.03, .1, .1)) +
scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
coord_cartesian(ylim = c(-.65, .65)) +
scale_x_continuous(breaks = c(10, 13, 16)) +
facet_grid(~ label, labeller = parcel_labeller) +
labs(x = "\nage", y = "mean standardized BOLD signal value\n", color = '', fill = '') +
dcbw +
theme(legend.position = c(.85, .15),
legend.spacing.y = unit(.01, 'cm'),
legend.margin = unit(0, "cm"))
self_other_int_plot_raw = neuro_model_data %>%
group_by(subjectID, age, label, domain, target, parcellation) %>%
mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
select(subjectID, age, label, domain, target, beta_std_avg) %>%
unique() %>%
spread(target, beta_std_avg) %>%
mutate(avg_diff = self - other) %>%
distinct(parcellation, age, label, .keep_all = T) %>%
ggplot(aes(x = age, y = avg_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(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2),
se = TRUE, size = 1.5) +
scale_color_manual(values = pal_social_academic) +
scale_fill_manual(values = pal_social_academic) +
scale_size_manual(values = c(.03, .1, .1)) +
scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
coord_cartesian(ylim = c(-.65, .65)) +
scale_x_continuous(breaks = c(10, 13, 16)) +
facet_grid(~ label, labeller = parcel_labeller) +
labs(x = "\nage", y = "mean standardized BOLD signal value\n", color = '', fill = '') +
dcbw +
theme(legend.position = c(.85, .15),
legend.spacing.y = unit(.01, 'cm'),
legend.margin = unit(0, "cm"))
(h3_raw = cowplot::plot_grid(int_plot_raw, soc_acad_int_plot_raw, self_other_int_plot_raw,
labels = c('A', 'B', 'C'), ncol = 3))
run models
domain x age
# specify model
model_domain_equation = formula(beta_std ~ 1 + age_c*domain*label +
age_c2*domain*label +
(1 + domain + age_c | subjectID) +
(1 + domain + age_c + age_c2 | parcellation))
model_domain_formula = lFormula(model_domain_equation, data = neuro_model_data_dummy)
# 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 = 50*(model_domain_numFx + model_domain_numRx + 1)^2
# run or load the model
if (file.exists("../data/model_domain.RDS")) {
model_domain = readRDS("../data/model_domain.RDS")
} else {
model_domain = lmer(model_domain_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_domain_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
saveRDS(model_domain, "../data/model_domain.RDS")
}
summary(model_domain)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: model_domain_equation
## Data: neuro_model_data_dummy
## Control:
## lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa",
## calc.derivs = FALSE)
##
## AIC BIC logLik deviance df.resid
## 462224.6 462579.8 -231077.3 462154.6 188542
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0506 -0.6008 -0.0084 0.5920 3.8505
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## parcellation (Intercept) 0.20821426 0.456305
## domain 0.00806375 0.089798 0.26
## age_c 0.00088083 0.029679 0.47 0.26
## age_c2 0.00003623 0.006019 -0.37 -0.42 -0.38
## subjectID (Intercept) 0.00215410 0.046412
## domain 0.00499073 0.070645 0.08
## age_c 0.00016868 0.012988 0.22 0.34
## Residual 0.66779327 0.817186
## Number of obs: 188577, groups: parcellation, 352; subjectID, 74
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -0.00470119 0.02776899 383.14722407
## age_c 0.00241423 0.00273527 149.71847282
## domain 0.01381497 0.01246246 193.31491209
## labelself -0.03361084 0.10898061 352.10270654
## labelsocial 0.24213531 0.07446145 352.05411001
## age_c2 0.00002738 0.00064418 574.64433087
## age_c:domain -0.00306884 0.00190350 6273.57253448
## age_c:labelself 0.00709832 0.00777864 352.26517211
## age_c:labelsocial 0.01841752 0.00531320 351.79581518
## domain:labelself 0.06467357 0.03525574 1012.95937367
## domain:labelsocial 0.11817131 0.02405744 1007.58093919
## domain:age_c2 -0.00434010 0.00091709 17310.53722166
## labelself:age_c2 0.00130288 0.00221804 352.85964350
## labelsocial:age_c2 -0.00428726 0.00151355 351.00849091
## age_c:domain:labelself -0.00618361 0.00666326 186980.53100214
## age_c:domain:labelsocial -0.00457075 0.00454456 186978.07006976
## domain:labelself:age_c2 0.00657677 0.00339852 186985.06357544
## domain:labelsocial:age_c2 0.00545640 0.00231693 186980.94510218
## t value Pr(>|t|)
## (Intercept) -0.169 0.865653
## age_c 0.883 0.378850
## domain 1.109 0.269011
## labelself -0.308 0.757952
## labelsocial 3.252 0.001257 **
## age_c2 0.043 0.966114
## age_c:domain -1.612 0.106967
## age_c:labelself 0.913 0.362109
## age_c:labelsocial 3.466 0.000593 ***
## domain:labelself 1.834 0.066886 .
## domain:labelsocial 4.912 0.00000105 ***
## domain:age_c2 -4.732 0.00000224 ***
## labelself:age_c2 0.587 0.557309
## labelsocial:age_c2 -2.833 0.004884 **
## age_c:domain:labelself -0.928 0.353400
## age_c:domain:labelsocial -1.006 0.314531
## domain:labelself:age_c2 1.935 0.052969 .
## domain:labelsocial:age_c2 2.355 0.018523 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1domain 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 = 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.RDS")) {
model_target = readRDS("../data/model_target.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.RDS")
}
summary(model_target)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: model_target_equation
## Data: neuro_model_data_dummy
## Control:
## lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa",
## calc.derivs = FALSE)
##
## AIC BIC logLik deviance df.resid
## 461839.9 463372.1 -230768.9 461537.9 188426
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1138 -0.6003 -0.0075 0.5917 4.0134
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## parcellation (Intercept) 0.208315479 0.456416
## age_c 0.000884423 0.029739 0.46
## target 0.004500598 0.067086 -0.11 -0.35
## domain 0.009993281 0.099966 0.23 0.30 -0.32
## age_c2 0.000038415 0.006198 -0.37 -0.37 0.53 -0.30
## age_c:target 0.000028694 0.005357 -0.30 0.06 0.78 -0.47
## age_c:domain 0.000099945 0.009997 0.75 0.58 -0.02 -0.33
## target:domain 0.002383879 0.048825 0.31 0.71 -0.60 0.84
## target:age_c2 0.000021163 0.004600 0.43 0.43 -0.86 0.19
## domain:age_c2 0.000007047 0.002655 0.22 -0.28 0.00 -0.65
## age_c:target:domain 0.000026276 0.005126 0.65 -0.10 0.38 -0.45
## target:domain:age_c2 0.000011578 0.003403 -0.85 -0.59 0.56 -0.26
## subjectID (Intercept) 0.001986582 0.044571
## age_c 0.000194035 0.013930 0.18
## target 0.001417266 0.037647 -0.16 0.13
## domain 0.004732643 0.068794 0.10 0.31 0.04
## age_c:target 0.000154815 0.012442 -0.09 -0.07 -0.11 -0.39
## age_c:domain 0.000540990 0.023259 -0.50 0.04 0.31 0.18
## target:domain 0.007392230 0.085978 0.18 0.45 0.11 -0.32
## age_c:target:domain 0.000894980 0.029916 0.06 -0.08 -0.43 -0.56
## Residual 0.664221997 0.814998
##
##
##
##
##
##
## 0.48
## -0.37 0.12
## -0.52 -0.43 0.02
## -0.82 -0.72 0.42 0.52
## -0.41 -0.11 0.51 -0.52 0.38
## -0.17 0.10 0.72 -0.46 0.08 0.74
## 0.69 0.56 -0.71 -0.52 -0.83 -0.31 -0.41
##
##
##
##
##
## 0.11
## 0.33 -0.08
## 0.44 -0.74 0.30
##
## Number of obs: 188577, groups: parcellation, 352; subjectID, 74
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -0.0039281 0.0277014 380.6628882
## age_c 0.0032238 0.0027573 146.1193893
## target 0.0038205 0.0095641 188.9704616
## domain 0.0134584 0.0129922 139.1236009
## labelself -0.0334801 0.1090019 352.1218175
## labelsocial 0.2423061 0.0744760 352.0732087
## age_c2 0.0001401 0.0006444 577.2447465
## age_c:target 0.0024352 0.0024481 75.1749279
## age_c:domain 0.0004922 0.0036977 58.0095169
## target:domain -0.0060469 0.0186190 137.3024152
## age_c:labelself 0.0071203 0.0077878 352.4916624
## age_c:labelsocial 0.0184586 0.0053195 352.0225729
## target:labelself 0.2552736 0.0322323 475.0766633
## target:labelsocial -0.0510006 0.0219886 472.0334880
## domain:labelself 0.0661994 0.0367017 447.8687911
## domain:labelsocial 0.1180557 0.0250466 445.6636084
## target:age_c2 -0.0017697 0.0009762 959.6250219
## domain:age_c2 -0.0030908 0.0010288 1434.4132242
## labelself:age_c2 0.0012939 0.0022421 366.8014351
## labelsocial:age_c2 -0.0043157 0.0015300 364.9159721
## age_c:target:domain 0.0028374 0.0054019 54.8431999
## age_c:target:labelself 0.0260003 0.0067656 3904.8268258
## age_c:target:labelsocial 0.0063011 0.0046146 3877.7360937
## age_c:domain:labelself -0.0060762 0.0070547 1437.5574380
## age_c:domain:labelsocial -0.0045981 0.0048125 1428.2933190
## target:domain:labelself -0.0220377 0.0572682 3796.0053311
## target:domain:labelsocial 0.0975762 0.0390508 3765.8994302
## target:domain:age_c2 -0.0001262 0.0018775 1743.8797920
## target:labelself:age_c2 -0.0107295 0.0035603 1444.4065785
## target:labelsocial:age_c2 0.0042813 0.0024277 1432.5625907
## domain:labelself:age_c2 0.0064675 0.0034473 3673.4163646
## domain:labelsocial:age_c2 0.0054693 0.0023503 3641.7187801
## age_c:target:domain:labelself -0.0055670 0.0133465 17643.1850585
## age_c:target:domain:labelsocial -0.0136281 0.0091027 17526.3680023
## target:domain:labelself:age_c2 0.0085129 0.0068268 10786.6270743
## target:domain:labelsocial:age_c2 -0.0013490 0.0046542 10695.0553010
## t value Pr(>|t|)
## (Intercept) -0.142 0.887311
## age_c 1.169 0.244227
## target 0.399 0.690001
## domain 1.036 0.302053
## labelself -0.307 0.758909
## labelsocial 3.253 0.001250 **
## age_c2 0.217 0.827965
## age_c:target 0.995 0.323062
## age_c:domain 0.133 0.894562
## target:domain -0.325 0.745850
## age_c:labelself 0.914 0.361192
## age_c:labelsocial 3.470 0.000585 ***
## target:labelself 7.920 0.000000000000017 ***
## target:labelsocial -2.319 0.020799 *
## domain:labelself 1.804 0.071948 .
## domain:labelsocial 4.713 0.000003261470680 ***
## target:age_c2 -1.813 0.070178 .
## domain:age_c2 -3.004 0.002707 **
## labelself:age_c2 0.577 0.564221
## labelsocial:age_c2 -2.821 0.005054 **
## age_c:target:domain 0.525 0.601525
## age_c:target:labelself 3.843 0.000123 ***
## age_c:target:labelsocial 1.365 0.172177
## age_c:domain:labelself -0.861 0.389222
## age_c:domain:labelsocial -0.955 0.339510
## target:domain:labelself -0.385 0.700396
## target:domain:labelsocial 2.499 0.012507 *
## target:domain:age_c2 -0.067 0.946409
## target:labelself:age_c2 -3.014 0.002627 **
## target:labelsocial:age_c2 1.764 0.078024 .
## domain:labelself:age_c2 1.876 0.060722 .
## domain:labelsocial:age_c2 2.327 0.020019 *
## age_c:target:domain:labelself -0.417 0.676598
## age_c:target:domain:labelsocial -1.497 0.134371
## target:domain:labelself:age_c2 1.247 0.212435
## target:domain:labelsocial:age_c2 -0.290 0.771944
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1compare models
anova(model_domain, model_target) %>%
as.data.frame() %>%
select(npar, AIC, BIC, 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),
BIC = round(BIC, 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 | BIC | X2 | X2 df | p |
|---|---|---|---|---|---|---|
| model_domain | 35 | 462224.6 | 462579.8 | – | – | – |
| model_target | 151 | 461839.9 | 463372.1 | 616.74 | 116 | < .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.004 [-0.058, 0.051] | 0.028 | -0.142 | 380.663 | .887 |
| Age | 0.003 [-0.002, 0.009] | 0.003 | 1.169 | 146.119 | .244 |
| Target | 0.004 [-0.015, 0.023] | 0.010 | 0.399 | 188.970 | .690 |
| Domain | 0.013 [-0.012, 0.039] | 0.013 | 1.036 | 139.124 | .302 |
| Label (self) | -0.033 [-0.248, 0.181] | 0.109 | -0.307 | 352.122 | .759 |
| Label (social) | 0.242 [0.096, 0.389] | 0.074 | 3.253 | 352.073 | .001 |
| Age2 | 0.000 [-0.001, 0.001] | 0.001 | 0.217 | 577.245 | .828 |
| Age x Target | 0.002 [-0.002, 0.007] | 0.002 | 0.995 | 75.175 | .323 |
| Age x Domain | 0.000 [-0.007, 0.008] | 0.004 | 0.133 | 58.010 | .895 |
| Target x Domain | -0.006 [-0.043, 0.031] | 0.019 | -0.325 | 137.302 | .746 |
| Age x Label (self) | 0.007 [-0.008, 0.022] | 0.008 | 0.914 | 352.492 | .361 |
| Age x Label (social) | 0.018 [0.008, 0.029] | 0.005 | 3.470 | 352.023 | .001 |
| Target x Label (self) | 0.255 [0.192, 0.319] | 0.032 | 7.920 | 475.077 | < .001 |
| Target x Label (social) | -0.051 [-0.094, -0.008] | 0.022 | -2.319 | 472.033 | .021 |
| Domain x Label (self) | 0.066 [-0.006, 0.138] | 0.037 | 1.804 | 447.869 | .072 |
| Domain x Label (social) | 0.118 [0.069, 0.167] | 0.025 | 4.713 | 445.664 | < .001 |
| Target x Age2 | -0.002 [-0.004, 0.000] | 0.001 | -1.813 | 959.625 | .070 |
| Domain x Age2 | -0.003 [-0.005, -0.001] | 0.001 | -3.004 | 1434.413 | .003 |
| Label (self) x Age2 | 0.001 [-0.003, 0.006] | 0.002 | 0.577 | 366.801 | .564 |
| Label (social) x Age2 | -0.004 [-0.007, -0.001] | 0.002 | -2.821 | 364.916 | .005 |
| Age x Target x Domain | 0.003 [-0.008, 0.014] | 0.005 | 0.525 | 54.843 | .602 |
| Age x Target x Label (self) | 0.026 [0.013, 0.039] | 0.007 | 3.843 | 3904.827 | < .001 |
| Age x Target x Label (social) | 0.006 [-0.003, 0.015] | 0.005 | 1.365 | 3877.736 | .172 |
| Age x Domain x Label (self) | -0.006 [-0.020, 0.008] | 0.007 | -0.861 | 1437.557 | .389 |
| Age x Domain x Label (social) | -0.005 [-0.014, 0.005] | 0.005 | -0.955 | 1428.293 | .340 |
| Target x Domain x Label (self) | -0.022 [-0.134, 0.090] | 0.057 | -0.385 | 3796.005 | .700 |
| Target x Domain x Label (social) | 0.098 [0.021, 0.174] | 0.039 | 2.499 | 3765.899 | .013 |
| Target x Domain x Age2 | -0.000 [-0.004, 0.004] | 0.002 | -0.067 | 1743.880 | .946 |
| Target x Label (self) x Age2 | -0.011 [-0.018, -0.004] | 0.004 | -3.014 | 1444.407 | .003 |
| Target x Label (social) x Age2 | 0.004 [-0.000, 0.009] | 0.002 | 1.764 | 1432.563 | .078 |
| Domain x Label (self) x Age2 | 0.006 [-0.000, 0.013] | 0.003 | 1.876 | 3673.416 | .061 |
| Domain x Label (social) x Age2 | 0.005 [0.001, 0.010] | 0.002 | 2.327 | 3641.719 | .020 |
| Age x Target x Domain x Label (self) | -0.006 [-0.032, 0.021] | 0.013 | -0.417 | 17643.185 | .677 |
| Age x Target x Domain x Label (social) | -0.014 [-0.031, 0.004] | 0.009 | -1.497 | 17526.368 | .134 |
| Target x Domain x Label (self) x Age2 | 0.009 [-0.005, 0.022] | 0.007 | 1.247 | 10786.627 | .212 |
| Target x Domain x Label (social) x Age2 | -0.001 [-0.010, 0.008] | 0.005 | -0.290 | 10695.055 | .772 |
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.027 | 0.010 | 1664.99 | 2.77 | 0.006 |
| social self > other | self | quadratic | -0.008 | 0.005 | 4879.37 | 1.75 | 0.081 |
| self social > academic | social | linear | -0.010 | 0.007 | 910.55 | 1.43 | 0.153 |
| self social > academic | social | quadratic | 0.002 | 0.003 | 7651.76 | 0.53 | 0.598 |
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(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))