Supplementary material analyses

This document provides the code for the analyses reported in:

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

This script reproduces the analyses reported in supplementary material.

load packages

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

define color palettes and labels

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

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

label_df = expand.grid(label = c("social", "self", "other"),
              target = c("self", "other"),
              domain = c("social", "academic"),
              pdss = 3,
              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') 

# mri exclusions
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 cleaned SPPC and task data

# demographics
race_ethnicity = readxl::read_xlsx("../data/Ethnicity_w1.xlsx") %>%
  extract(SID, "subjectID", "L([0-9]{3})-.*") %>%
  mutate(subjectID = sprintf("s%s", subjectID),
         hispanic_latinx = ifelse(Ethnicity == 1, "yes", "no"),
         race = recode(Race, `1` = "White", `2` = "Multiracial", `3` = "Pacific Islander",
                       `4` = "Native American", `5` = "Black or African American")) %>%
  select(subjectID, hispanic_latinx, race)

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) %>%
  full_join(., race_ethnicity)

# self-report
sppc = readRDS("../data/sppc.RDS") %>%
  filter(!domain == "self_worth")

sppc_self_worth = readRDS("../data/sppc.RDS") %>%
  filter(domain == "self_worth") %>%
  rename("self_worth" = competence) %>%
  select(-domain, -importance)

# task
task = readRDS("../data/task.RDS") %>%
  ungroup() %>%
  mutate(percent_stdgm = (task_percent - mean(task_percent, na.rm = TRUE)) / 10, #center at grand mean, in units of 10%
         percent_std50 = (task_percent - 50) / 10,
         percent_std = (task_percent - 50) / 10) %>% #center at 50, in units of 10%
  group_by(target, domain) %>%
  mutate(percent_stdc = (task_percent - mean(task_percent, na.rm = TRUE)) / 10) #center at condition mean, in units of 10%

# merge data
merged = parcellations_ex %>%
  full_join(., task, by = c("subjectID", "wave", "domain", "target")) %>%
  full_join(., sppc, by = c("subjectID", "wave", "domain")) %>%
  full_join(., sppc_self_worth, by = c("subjectID", "wave")) %>%
  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_pdss = merged %>%
  filter(!is.na(beta)) %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, target, domain, parcellation, label, beta, beta_std)

neuro_model_data_ex_pdss = neuro_model_data_pdss  %>%
  na.omit()

comp_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, domain, competence) %>%
  unique() %>%
  na.omit()

imp_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, domain, importance) %>%
  unique() %>%
  na.omit()

self_worth_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, self_worth) %>%
  unique() %>%
  na.omit()

status_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, target, domain, task_percent, contains("percent")) %>%
  unique() %>%
  na.omit()

# dummy code target and domain
neuro_model_data_dummy_pdss = neuro_model_data_ex_pdss %>%
  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

merged %>% 
  select(subjectID, hispanic_latinx, race) %>%
  unique() %>%
  group_by(hispanic_latinx, race) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = (n/total) * 100) %>%
  select(-n, -total) %>%
  spread(hispanic_latinx, percent) %>%
  kable(format = 'pandoc', digits = 1)

race	no	yes
Black or African American	3.3	NA
Multiracial	11.1	5.6
Native American	4.4	1.1
Pacific Islander	4.4	1.1
White	43.3	25.6

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 %>% 
  select(subjectID, inclusion, hispanic_latinx, race) %>%
  unique() %>%
  group_by(inclusion, hispanic_latinx, race) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = (n/total) * 100) %>%
  select(-n, -total) %>%
  spread(hispanic_latinx, percent) %>%
  kable(format = 'pandoc', digits = 1)

inclusion	race	no	yes
completed MRI	Black or African American	0.8	NA
completed MRI	Multiracial	7.6	2.5
completed MRI	Native American	3.4	0.8
completed MRI	Pacific Islander	2.5	0.8
completed MRI	White	27.7	16.0
didn’t complete MRI	Black or African American	2.5	NA
didn’t complete MRI	Multiracial	1.7	NA
didn’t complete MRI	Native American	NA	0.8
didn’t complete MRI	Pacific Islander	0.8	NA
didn’t complete MRI	White	8.4	3.4
excluded from MRI	Multiracial	0.8	2.5
excluded from MRI	Native American	0.8	NA
excluded from MRI	Pacific Islander	0.8	0.8
excluded from MRI	White	11.8	2.5

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

age and puberty correlations

cor_fun = function(data) pmap(var.names, ~ cor.test(data[[.x]], data[[.y]])) %>% 
  map_df(broom::tidy) %>% 
  cbind(var.names, .)

age_pds = merged %>%
  select(subjectID, wave, gender, age, pdss) %>%
  gather(dev, value, age, pdss) %>%
  unite(wave_dev, dev, wave, sep = " ", remove = TRUE) %>%
  rownames_to_column() %>%
  spread(wave_dev, value) %>%
  group_by(subjectID) %>%
  fill(everything(), .direction = "up") %>%
  fill(everything(), .direction = "down") %>%
  select(-rowname) %>%
  unique() %>%
  ungroup()

var.names = tidystringdist::tidy_comb_all(names(select(age_pds, -c(subjectID, gender)))) %>%
  filter(!(grepl("age", V1) & grepl("age", V2)))

age_pds %>%
  nest(-c(gender)) %>%
  mutate(
    test = map(data, cor_fun)
  ) %>% 
  unnest(test, .drop = TRUE) %>%
  mutate(r = sprintf("%.02f", estimate), #sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high)
         r = ifelse(p.value < .05 & p.value >= .01, paste0(r,"*"),
             ifelse(p.value < .01 & p.value >= .001, paste0(r,"**"),
             ifelse(p.value < .001, paste0(r,"***"), r))),
         r = gsub("0\\.", "\\.", r)) %>%
  select(gender, V1, V2, r) %>%
  rename("Sex" = gender) %>%
  unique() %>%
  mutate(V1 = gsub("age", "Age", V1),
         V2 = gsub("age", "Age", V2),
         V1 = gsub("pdss", "PDS-S", V1),
         V2 = gsub("pdss", "PDS-S", V2),
         V1 = gsub("t", "", V1),
         V2 = gsub("t", "", V2)) %>%
  spread(V2, r) %>%
  arrange(Sex) %>%
  filter(Sex == "F") %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "–-", .))) %>%
  kable(format = "pandoc")

Sex	V1	PDS-S 1	PDS-S 2	PDS-S 3
F	Age 1	.50***	.20	.01
F	Age 2	-.12	.40*	.30
F	Age 3	-.11	.02	.24
F	PDS-S 1	–-	.42*	-.08
F	PDS-S 2	–-	–-	.25

age_pds %>%
  nest(-c(gender)) %>%
  mutate(
    test = map(data, cor_fun)
  ) %>% 
  unnest(test, .drop = TRUE) %>%
  mutate(r = sprintf("%.02f", estimate), #sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high)
         r = ifelse(p.value < .05 & p.value >= .01, paste0(r,"*"),
             ifelse(p.value < .01 & p.value >= .001, paste0(r,"**"),
             ifelse(p.value < .001, paste0(r,"***"), r))),
         r = gsub("0\\.", "\\.", r)) %>%
  select(gender, V1, V2, r) %>%
  rename("Sex" = gender) %>%
  unique() %>%
  mutate(V1 = gsub("age", "Age", V1),
         V2 = gsub("age", "Age", V2),
         V1 = gsub("pdss", "PDS-S", V1),
         V2 = gsub("pdss", "PDS-S", V2),
         V1 = gsub("t", "", V1),
         V2 = gsub("t", "", V2)) %>%
  spread(V1, r) %>%
  arrange(Sex) %>%
  filter(Sex == "M") %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "–-", .))) %>%
  kable(format = "pandoc")

Sex	V2	Age 1	Age 2	Age 3	PDS-S 1	PDS-S 2
M	PDS-S 1	-.32*	-.37	-.26	–-	–-
M	PDS-S 2	.22	.41*	.50*	.16	–-
M	PDS-S 3	-.24	.37	.41	.30	.62**

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 = "\npubertal stage (PSD score)", 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 = "\npubertal stage (PSD score)", 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"))

neural puberty models

visualize raw

Visualize the developmental trajectories using the raw data

hypothesis 0

domain_parc_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, 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_y_continuous(breaks = seq(-2, 2, 1)) + 
  coord_cartesian(ylim = c(-2, 2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', 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_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, 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_y_continuous(breaks = seq(-2, 2, 1)) + 
  coord_cartesian(ylim = c(-2, 2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', 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_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, 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_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.6, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.6, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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"))

(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_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, 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_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.5, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(beta_std_avf = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, beta_std_avf) %>%
  unique() %>%
  spread(target, beta_std_avf) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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(-.5, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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"))

(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_pdss %>%
  distinct(parcellation, target, domain, pdss, beta_std, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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")) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', 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_pdss %>%
  group_by(subjectID, pdss, label, target, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_pdss %>%
  group_by(subjectID, pdss, label, domain, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, target, beta_std_avg) %>%
  unique() %>%
  spread(target, beta_std_avg) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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

Estimate full models

domain x puberty

model_domain_equation = formula(beta_std ~ 1 + pdss_c*domain*label +
                                  pdss_c2*domain*label +
                                  (1 + pdss_c*domain || subjectID) +
                                  (1 + pdss_c*domain + pdss_c2*domain || parcellation))

model_domain_formula = lFormula(model_domain_equation, data = neuro_model_data_dummy_pdss)

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

if (file.exists("../data/model_domain_pdss.RDS")) {
  model_domain = readRDS("../data/model_domain_pdss.RDS")
} else {
  model_domain = lmer(model_domain_equation, data = neuro_model_data_dummy_pdss, 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_pdss.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_pdss
## Control: 
## lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa",  
##     calc.derivs = FALSE)
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  438496.7  439001.4 -219198.4  438396.7    178707 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9960 -0.6001 -0.0084  0.5907  3.9156 
## 
## Random effects:
##  Groups       Name           Variance  Std.Dev. Corr                         
##  parcellation (Intercept)    0.1915074 0.43762                               
##               pdss_c         0.0019204 0.04382   0.32                        
##               domain         0.0058018 0.07617   0.20 -0.17                  
##               pdss_c2        0.0004336 0.02082  -0.08  0.26  0.28            
##               pdss_c:domain  0.0002954 0.01719   0.87  0.57  0.18  0.39      
##               domain:pdss_c2 0.0001758 0.01326   0.08  0.79  0.37  0.49  0.37
##  subjectID    (Intercept)    0.0020536 0.04532                               
##               pdss_c         0.0007719 0.02778   0.14                        
##               domain         0.0048433 0.06959   0.38  0.46                  
##               pdss_c:domain  0.0029625 0.05443  -0.62 -0.05 -0.49            
##  Residual                    0.6687116 0.81775                               
## Number of obs: 178757, groups:  parcellation, 352; subjectID, 71
## 
## Fixed effects:
##                                 Estimate    Std. Error            df t value
## (Intercept)                  -0.00864707    0.02674359  387.34171896  -0.323
## pdss_c                        0.00604326    0.00533745  104.27381117   1.132
## domain                        0.00685629    0.01313967  117.55647744   0.522
## labelself                    -0.01026887    0.10459828  352.12973718  -0.098
## labelsocial                   0.20896198    0.07146475  352.03292870   2.924
## pdss_c2                       0.00004446    0.00286432  553.09450523   0.016
## pdss_c:domain                 0.00789460    0.00886178   66.13809819   0.891
## pdss_c:labelself              0.00387981    0.01270059  353.20284903   0.305
## pdss_c:labelsocial            0.01506527    0.00867452  352.63124023   1.737
## domain:labelself              0.12529333    0.03347653  561.68895019   3.743
## domain:labelsocial            0.14395617    0.02281197  555.63602467   6.311
## domain:pdss_c2               -0.01128767    0.00509669  591.57331074  -2.215
## labelself:pdss_c2            -0.01231471    0.00821216  355.97912813  -1.500
## labelsocial:pdss_c2          -0.00096279    0.00559962  353.05486283  -0.172
## pdss_c:domain:labelself       0.02131053    0.01519558 2372.06515845   1.402
## pdss_c:domain:labelsocial     0.00706300    0.01037100 2361.46413438   0.681
## domain:labelself:pdss_c2     -0.00257839    0.01350342 2736.08663383  -0.191
## domain:labelsocial:pdss_c2    0.00459873    0.00919836 2703.31787444   0.500
##                                  Pr(>|t|)    
## (Intercept)                      0.746618    
## pdss_c                           0.260132    
## domain                           0.602791    
## labelself                        0.921850    
## labelsocial                      0.003680 ** 
## pdss_c2                          0.987623    
## pdss_c:domain                    0.376234    
## pdss_c:labelself                 0.760179    
## pdss_c:labelsocial               0.083309 .  
## domain:labelself                 0.000201 ***
## domain:labelsocial         0.000000000569 ***
## domain:pdss_c2                   0.027160 *  
## labelself:pdss_c2                0.134612    
## labelsocial:pdss_c2              0.863585    
## pdss_c:domain:labelself          0.160922    
## pdss_c:domain:labelsocial        0.495917    
## domain:labelself:pdss_c2         0.848584    
## domain:labelsocial:pdss_c2       0.617150    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

domain x target x puberty

model_target_equation = formula(beta_std ~ 1 + pdss_c*target*domain*label +
                                  pdss_c2*target*domain*label +
                                  (1 + pdss_c*target*domain || subjectID) +
                                  (1 + pdss_c*target*domain + pdss_c2*target*domain || parcellation))

model_target_formula = lFormula(model_target_equation, data = neuro_model_data_dummy_pdss)

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

if (file.exists("../data/model_target_pdss.RDS")) {
  model_target = readRDS("../data/model_target_pdss.RDS")
} else {
  model_target = lmer(model_target_equation, data = neuro_model_data_dummy_pdss, 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_pdss.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_pdss
## Control: 
## lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa",  
##     calc.derivs = FALSE)
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  438286.9  439811.1 -218992.5  437984.9    178606 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0624 -0.5995 -0.0084  0.5900  3.9871 
## 
## Random effects:
##  Groups       Name                  Variance    Std.Dev. Corr             
##  parcellation (Intercept)           0.191534957 0.437647                  
##               pdss_c                0.001939251 0.044037  0.32            
##               target                0.002101585 0.045843  0.53 -0.29      
##               domain                0.005697644 0.075483  0.20 -0.16  0.31
##               pdss_c2               0.000449304 0.021197 -0.08  0.25 -0.74
##               pdss_c:target         0.000009699 0.003114 -0.42  0.62 -0.47
##               pdss_c:domain         0.000282659 0.016812  0.90  0.58  0.17
##               target:domain         0.000847892 0.029119 -0.50  0.46 -0.76
##               target:pdss_c2        0.000159009 0.012610 -0.53  0.46 -0.63
##               domain:pdss_c2        0.000203644 0.014270  0.05  0.70 -0.54
##               pdss_c:target:domain  0.000202909 0.014245 -0.02  0.13  0.23
##               target:domain:pdss_c2 0.000225911 0.015030  0.18 -0.02  0.06
##  subjectID    (Intercept)           0.002058021 0.045365                  
##               pdss_c                0.000794142 0.028181  0.14            
##               target                0.001079353 0.032854 -0.10  0.10      
##               domain                0.004642281 0.068134  0.36  0.40 -0.02
##               pdss_c:target         0.000631602 0.025132 -0.04 -0.26 -0.33
##               pdss_c:domain         0.003204006 0.056604 -0.60  0.01  0.25
##               target:domain         0.007155311 0.084589  0.11  0.27  0.16
##               pdss_c:target:domain  0.003097992 0.055660  0.23 -0.12 -0.49
##  Residual                           0.666315705 0.816282                  
##                                                 
##                                                 
##                                                 
##                                                 
##                                                 
##   0.21                                          
##  -0.04  0.29                                    
##   0.17  0.30 -0.09                              
##   0.17  0.67  0.84 -0.11                        
##  -0.77  0.05  0.63 -0.30  0.46                  
##   0.39  0.68  0.60  0.38  0.80  0.13            
##  -0.81 -0.72 -0.02 -0.18 -0.44  0.55 -0.56      
##   0.96  0.43  0.05  0.24  0.34 -0.66  0.59 -0.90
##                                                 
##                                                 
##                                                 
##                                                 
##  -0.49                                          
##  -0.51  0.24                                    
##  -0.38  0.31  0.52                              
##  -0.16  0.42 -0.60 -0.22                        
##                                                 
## Number of obs: 178757, groups:  parcellation, 352; subjectID, 71
## 
## Fixed effects:
##                                        Estimate    Std. Error            df
## (Intercept)                          -0.0079434     0.0267377   387.1349345
## pdss_c                                0.0070218     0.0053534   105.4100607
## target                               -0.0100045     0.0091956   156.8891860
## domain                                0.0050395     0.0128061   120.9093238
## labelself                            -0.0100628     0.1046023   352.1414269
## labelsocial                           0.2090091     0.0714675   352.0446218
## pdss_c2                               0.0003423     0.0028500   536.2910191
## pdss_c:target                         0.0048547     0.0054353    48.6340142
## pdss_c:domain                         0.0080526     0.0088379    60.6677770
## target:domain                        -0.0028657     0.0189577   143.7290923
## pdss_c:labelself                      0.0040014     0.0127346   354.5810794
## pdss_c:labelsocial                    0.0151135     0.0086977   354.0043656
## target:labelself                      0.2138397     0.0301712   934.3619416
## target:labelsocial                   -0.0187090     0.0205468   922.0088785
## domain:labelself                      0.1258877     0.0333471   559.7949984
## domain:labelsocial                    0.1439768     0.0227233   553.7118900
## target:pdss_c2                       -0.0022266     0.0041311   250.6126935
## domain:pdss_c2                       -0.0093440     0.0049267   488.0487025
## labelself:pdss_c2                    -0.0124303     0.0082565   364.8314190
## labelsocial:pdss_c2                  -0.0010014     0.0056300   361.8618792
## pdss_c:target:domain                 -0.0142544     0.0114306    59.9882324
## pdss_c:target:labelself               0.0265231     0.0146329 42548.1841208
## pdss_c:target:labelsocial             0.0083570     0.0099859 42393.8654519
## pdss_c:domain:labelself               0.0214696     0.0151471  2436.3188028
## pdss_c:domain:labelsocial             0.0070003     0.0103377  2425.2295300
## target:domain:labelself              -0.1285312     0.0567210 10051.4166844
## target:domain:labelsocial             0.0550301     0.0386103  9909.8525590
## target:domain:pdss_c2                 0.0018253     0.0082326   254.5687140
## target:labelself:pdss_c2             -0.0215736     0.0134457  3117.6038304
## target:labelsocial:pdss_c2            0.0016450     0.0091588  3079.8050178
## domain:labelself:pdss_c2             -0.0026889     0.0135386  2283.9185674
## domain:labelsocial:pdss_c2            0.0046446     0.0092224  2256.5412924
## pdss_c:target:domain:labelself       -0.0413546     0.0294226 11161.6901215
## pdss_c:target:domain:labelsocial     -0.0250807     0.0200792 11109.7223056
## target:domain:labelself:pdss_c2       0.1010612     0.0264604  8145.7108099
## target:domain:labelsocial:pdss_c2     0.0199685     0.0180221  8046.8092407
##                                   t value         Pr(>|t|)    
## (Intercept)                        -0.297         0.766561    
## pdss_c                              1.312         0.192492    
## target                             -1.088         0.278280    
## domain                              0.394         0.694625    
## labelself                          -0.096         0.923416    
## labelsocial                         2.925         0.003673 ** 
## pdss_c2                             0.120         0.904435    
## pdss_c:target                       0.893         0.376166    
## pdss_c:domain                       0.911         0.365826    
## target:domain                      -0.151         0.880061    
## pdss_c:labelself                    0.314         0.753545    
## pdss_c:labelsocial                  1.738         0.083144 .  
## target:labelself                    7.088 0.00000000000269 ***
## target:labelsocial                 -0.911         0.362767    
## domain:labelself                    3.775         0.000177 ***
## domain:labelsocial                  6.336 0.00000000048848 ***
## target:pdss_c2                     -0.539         0.590373    
## domain:pdss_c2                     -1.897         0.058471 .  
## labelself:pdss_c2                  -1.506         0.133056    
## labelsocial:pdss_c2                -0.178         0.858928    
## pdss_c:target:domain               -1.247         0.217232    
## pdss_c:target:labelself             1.813         0.069905 .  
## pdss_c:target:labelsocial           0.837         0.402666    
## pdss_c:domain:labelself             1.417         0.156491    
## pdss_c:domain:labelsocial           0.677         0.498367    
## target:domain:labelself            -2.266         0.023471 *  
## target:domain:labelsocial           1.425         0.154110    
## target:domain:pdss_c2               0.222         0.824711    
## target:labelself:pdss_c2           -1.604         0.108706    
## target:labelsocial:pdss_c2          0.180         0.857474    
## domain:labelself:pdss_c2           -0.199         0.842585    
## domain:labelsocial:pdss_c2          0.504         0.614580    
## pdss_c:target:domain:labelself     -1.406         0.159889    
## pdss_c:target:domain:labelsocial   -1.249         0.211659    
## target:domain:labelself:pdss_c2     3.819         0.000135 ***
## target:domain:labelsocial:pdss_c2   1.108         0.267895    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

compare

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	438496.7	–	–	–
model_target	151	438286.9	411.78	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 (puberty 3, 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("pdss_c", "Puberty", 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 (puberty 3, label (control))	-0.008 [-0.060, 0.044]	0.027	-0.297	387.135	.767
Puberty	0.007 [-0.003, 0.018]	0.005	1.312	105.410	.192
Target	-0.010 [-0.028, 0.008]	0.009	-1.088	156.889	.278
Domain	0.005 [-0.020, 0.030]	0.013	0.394	120.909	.695
Label (self)	-0.010 [-0.215, 0.195]	0.105	-0.096	352.141	.923
Label (social)	0.209 [0.069, 0.349]	0.071	2.925	352.045	.004
Puberty2	0.000 [-0.005, 0.006]	0.003	0.120	536.291	.904
Puberty x Target	0.005 [-0.006, 0.016]	0.005	0.893	48.634	.376
Puberty x Domain	0.008 [-0.009, 0.025]	0.009	0.911	60.668	.366
Target x Domain	-0.003 [-0.040, 0.034]	0.019	-0.151	143.729	.880
Puberty x Label (self)	0.004 [-0.021, 0.029]	0.013	0.314	354.581	.754
Puberty x Label (social)	0.015 [-0.002, 0.032]	0.009	1.738	354.004	.083
Target x Label (self)	0.214 [0.155, 0.273]	0.030	7.088	934.362	< .001
Target x Label (social)	-0.019 [-0.059, 0.022]	0.021	-0.911	922.009	.363
Domain x Label (self)	0.126 [0.061, 0.191]	0.033	3.775	559.795	< .001
Domain x Label (social)	0.144 [0.099, 0.189]	0.023	6.336	553.712	< .001
Target x Puberty2	-0.002 [-0.010, 0.006]	0.004	-0.539	250.613	.590
Domain x Puberty2	-0.009 [-0.019, 0.000]	0.005	-1.897	488.049	.058
Label (self) x Puberty2	-0.012 [-0.029, 0.004]	0.008	-1.506	364.831	.133
Label (social) x Puberty2	-0.001 [-0.012, 0.010]	0.006	-0.178	361.862	.859
Puberty x Target x Domain	-0.014 [-0.037, 0.008]	0.011	-1.247	59.988	.217
Puberty x Target x Label (self)	0.027 [-0.002, 0.055]	0.015	1.813	42548.184	.070
Puberty x Target x Label (social)	0.008 [-0.011, 0.028]	0.010	0.837	42393.865	.403
Puberty x Domain x Label (self)	0.021 [-0.008, 0.051]	0.015	1.417	2436.319	.156
Puberty x Domain x Label (social)	0.007 [-0.013, 0.027]	0.010	0.677	2425.230	.498
Target x Domain x Label (self)	-0.129 [-0.240, -0.017]	0.057	-2.266	10051.417	.023
Target x Domain x Label (social)	0.055 [-0.021, 0.131]	0.039	1.425	9909.853	.154
Target x Domain x Puberty2	0.002 [-0.014, 0.018]	0.008	0.222	254.569	.825
Target x Label (self) x Puberty2	-0.022 [-0.048, 0.005]	0.013	-1.604	3117.604	.109
Target x Label (social) x Puberty2	0.002 [-0.016, 0.020]	0.009	0.180	3079.805	.857
Domain x Label (self) x Puberty2	-0.003 [-0.029, 0.024]	0.014	-0.199	2283.919	.843
Domain x Label (social) x Puberty2	0.005 [-0.013, 0.023]	0.009	0.504	2256.541	.615
Puberty x Target x Domain x Label (self)	-0.041 [-0.099, 0.016]	0.029	-1.406	11161.690	.160
Puberty x Target x Domain x Label (social)	-0.025 [-0.064, 0.014]	0.020	-1.249	11109.722	.212
Target x Domain x Label (self) x Puberty2	0.101 [0.049, 0.153]	0.026	3.819	8145.711	< .001
Target x Domain x Label (social) x Puberty2	0.020 [-0.015, 0.055]	0.018	1.108	8046.809	.268

visualize fitted models

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

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

neuro_plot_data_puberty = with(neuro_model_data_dummy_pdss, 
                    expand.grid(target = unique(target), 
                                domain = unique(domain),
                                parcellation = unique(parcellation),
                                pdss = unique(pdss),
                                stringsAsFactors = F)) %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other')),
         pdss_c = pdss - 3, 
         pdss_c2 = pdss_c^2, 
         subjectID = NA)

neuro_plot_data_puberty$expected = predict(model_target, newdata = neuro_plot_data_puberty, re.form = reForm)
neuro_plot_data_puberty$expected_mean = predict(model_target, newdata = neuro_plot_data_puberty, re.form = NA)

neuro_plot_data_puberty = neuro_plot_data_puberty %>%
  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_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, 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_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, 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_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = domain)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  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) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .4)) +
  scale_fill_manual(values = "lightgrey") +
  scale_color_manual(values = pal_social_academic) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_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 = FALSE, color = pal_social_academic[2], size = 1.5) + 
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  scale_fill_manual(values = "lightgrey") +
  coord_cartesian(ylim = c(-.25, .4)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, target, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = target)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  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) +
  scale_y_continuous(breaks = seq(-.1, .3, .1)) +
  coord_cartesian(ylim = c(-.15, .3)) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = "lightgrey") +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_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 = FALSE, color = pal_self_other[2], size = 1.5) + 
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.1, .3, .1)) +
  scale_fill_manual(values = "lightgrey") +
  coord_cartesian(ylim = c(-.15, .3)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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) + 
  #geom_point() + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  group_by(subjectID, pdss, label, target, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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, .45)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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_puberty %>%
  group_by(subjectID, pdss, label, domain, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, 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, .45)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", 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))

self/other evaluation task analyses

visualize raw

age

(status_plot_raw_age = status_model_data %>%
  filter(!is.na(task_percent)) %>%
  ggplot(aes(age, task_percent, color = domain, fill = domain, linetype = target)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_x_continuous(limits = c(min(status_model_data$age), 18), breaks = c(10, 12, 14, 16, 18)) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PSD score)", y = "social status or academic proficiency\n(% endorsed)\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(status_plot_raw = status_model_data %>%
  filter(!is.na(task_percent)) %>%
  ggplot(aes(pdss, task_percent, color = domain, fill = domain, linetype = target)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PDS score)", y = "social status or academic proficiency\n(% endorsed)\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(task_percent ~ subjectID, data = status_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)

## [1] 0.18

age

`status model 1` = lmer(task_percent ~ 1 + target*domain + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 2 age` = lmer(task_percent ~ 1 + target*domain + age_c + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 3 age` = lmer(task_percent ~ 1 + target*domain + age_c + age_c2 + 
                         (1 + target * domain | subjectID), 
                         data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 4 age` = lmer(task_percent ~ 1 + target*domain*age_c + 
                            (1 + target * domain | subjectID), 
                           data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 5 age` = lmer(task_percent ~ 1 + target*domain*age_c + 
                            target*domain*age_c2 + 
                            (1 + target * domain | subjectID), 
                            data = status_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`status model 2 puberty` = lmer(task_percent ~ 1 + target*domain + pdss_c + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 3 puberty` = lmer(task_percent ~ 1 + target*domain + pdss_c + pdss_c2 + 
                         (1 + target * domain | subjectID), 
                         data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 4 puberty` = lmer(task_percent ~ 1 + target*domain*pdss_c + 
                            (1 + target * domain | subjectID), 
                           data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 5 puberty` = lmer(task_percent ~ 1 + target*domain*pdss_c + target*domain*pdss_c2 + 
                            (1 + target * domain | subjectID), 
                            data = status_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`status model 1`, `status model 2 age`, `status model 3 age`, `status model 4 age`, `status model 5 age`) %>%
  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
status model 1	15	5237.97	–	–	–
status model 2 age	16	5239.91	0.05	1	.815
status model 3 age	17	5241.69	0.22	1	.637
status model 4 age	19	5233.46	12.23	2	.002
status model 5 age	23	5231.07	10.39	4	.034

puberty

anova(`status model 1`, `status model 2 puberty`, `status model 3 puberty`, `status model 4 puberty`, `status model 5 puberty`) %>%
  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
status model 1	15	5237.97	–	–	–
status model 2 puberty	16	5238.96	1	1	.317
status model 3 puberty	17	5238.77	2.2	1	.138
status model 4 puberty	19	5232.81	9.95	2	.007
status model 5 puberty	23	5237.02	3.8	4	.434

summarize best fitting models

age

`status model 5 age` %>%
  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", p)),
         term = gsub("\\(Intercept\\)", "Intercept (academic, other, age 13)", term),
         term = gsub("domainsocial", "Domain (social)", term),
         term = gsub("targetself", "Target (self)", 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("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  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 (academic, other, age 13)	57.87 [51.28, 64.47]	3.36	17.21	174.40	< .001
Target (self)	4.82 [-4.34, 13.97]	4.67	1.03	182.28	.304
Domain (social)	2.29 [-6.42, 11.01]	4.45	0.52	188.62	.607
Age	-0.56 [-1.65, 0.53]	0.55	-1.01	365.43	.312
Age2	0.14 [-0.42, 0.69]	0.28	0.48	347.30	.631
Target (self) x Domain (social)	18.51 [6.41, 30.61]	6.17	3.00	194.63	.003
Target (self) x Age	-0.68 [-2.22, 0.86]	0.79	-0.87	366.83	.386
Domain (social) x Age	0.55 [-0.97, 2.07]	0.78	0.70	371.77	.481
Target (self) x Age2	0.43 [-0.35, 1.21]	0.40	1.07	348.69	.283
Domain (social) x Age2	-0.48 [-1.26, 0.29]	0.40	-1.22	352.91	.222
Target (self) x Domain (social) x Age	2.09 [-0.04, 4.22]	1.09	1.92	372.86	.056
Target (self) x Domain (social) x Age2	-0.64 [-1.73, 0.46]	0.56	-1.14	356.85	.255

print(VarCorr(`status model 5 age`), comp = c("Variance","Std.Dev."), digits = 5)

##  Groups    Name                    Variance Std.Dev. Corr                
##  subjectID (Intercept)              377.05  19.418                       
##            targetself               689.50  26.258   -0.815              
##            domainsocial             565.08  23.771   -0.866  0.798       
##            targetself:domainsocial 1040.78  32.261    0.772 -0.936 -0.905
##  Residual                           239.13  15.464

puberty

`status model 4 puberty` %>%
  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", p)),
         term = gsub("\\(Intercept\\)", "Intercept (academic, other, stage 3)", term),
         term = gsub("domainsocial", "Domain (social)", term),
         term = gsub("targetself", "Target (self)", term),
         term = gsub("pdss_c", "Puberty", 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("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  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 (academic, other, stage 3)	59.23 [54.19, 64.28]	2.57	23.01	78.80	< .001
Target (self)	8.48 [1.51, 15.45]	3.56	2.38	80.81	.02
Domain (social)	-1.54 [-8.03, 4.94]	3.31	-0.47	80.29	.642
Puberty	-1.68 [-3.91, 0.54]	1.13	-1.48	380.92	.139
Target (self) x Domain (social)	12.33 [3.43, 21.24]	4.54	2.71	80.80	.008
Target (self) x Puberty	-1.14 [-4.28, 2.00]	1.60	-0.71	379.97	.478
Domain (social) x Puberty	1.37 [-1.72, 4.46]	1.57	0.87	382.90	.385
Target (self) x Domain (social) x Puberty	3.53 [-0.78, 7.83]	2.20	1.60	383.96	.109

print(VarCorr(`status model 4 puberty`), comp = c("Variance","Std.Dev."), digits = 5)

##  Groups    Name                    Variance Std.Dev. Corr                
##  subjectID (Intercept)              367.36  19.167                       
##            targetself               687.68  26.224   -0.825              
##            domainsocial             561.36  23.693   -0.869  0.796       
##            targetself:domainsocial 1030.63  32.103    0.797 -0.941 -0.908
##  Residual                           244.39  15.633

visualize fitted

age

model_to_plot = `status model 5 age`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
status_model_data$expected_mean = predict(model_to_plot, newdata = status_model_data, re.form = NA)

(status_plot_age = status_model_data %>%
  filter(!is.na(target)) %>%
  ggplot(aes(age, expected_mean, color = domain, linetype = target)) +
    geom_line(aes(age, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_point(aes(age, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE, size = 1.25) +
    scale_x_continuous(limits = c(min(status_model_data$age), 18), breaks = c(10, 12, 14, 16, 18)) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    coord_cartesian(ylim = c(0, 100)) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PSD score)", y = "social status or academic proficiency\n") +
    dcbw +
    theme(legend.position = "none"))

puberty

model_to_plot = `status model 4 puberty`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
status_model_data$expected_mean = predict(model_to_plot, newdata = status_model_data, re.form = NA)

(status_plot_pdss = status_model_data %>%
  filter(!is.na(target)) %>%
  ggplot(aes(pdss, expected_mean, color = domain, linetype = target)) +
    geom_line(aes(pdss, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_point(aes(pdss, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE, size = 1.25) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    coord_cartesian(ylim = c(0, 100)) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PDS score)", y = "social status or academic proficiency\n") +
    dcbw +
    theme(legend.position = c(.8, .14),
          legend.spacing.x = unit(.1, 'cm'),
          legend.box = "horizontal"))

combined

cowplot::plot_grid(status_plot_age, status_plot_pdss, labels = c('A', 'B'), ncol = 2)

SPPC competence analysis

visualize raw

age

(comp_plot_raw_age = comp_model_data %>%
  ggplot(aes(age, competence, color = domain, fill = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\nage (years)", y = "mean competence rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(comp_plot_raw = comp_model_data %>%
  ggplot(aes(pdss, competence, color = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "mean competence rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(competence ~ subjectID, data = comp_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)

## [1] 0.36

age

`comp model 1` = lmer(competence ~ 1 + domain + 
                              (1 | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 2 age` = lmer(competence ~ 1 + domain + age_c + 
                             (1 + age_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 3 age` = lmer(competence ~ 1 + domain + age_c + age_c2 + 
                             (1 + age_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 4 age` = lmer(competence ~ 1 + domain*age_c + 
                              (1 + age_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 5 age` = lmer(competence ~ 1 + domain*age_c + domain*age_c2 + 
                              (1 + age_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`comp model 2 puberty` = lmer(competence ~ 1 + domain + pdss_c + 
                             (1 + pdss_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 3 puberty` = lmer(competence ~ 1 + domain + pdss_c + pdss_c2 + 
                             (1 + pdss_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 4 puberty` = lmer(competence ~ 1 + domain*pdss_c + 
                              (1 + pdss_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 5 puberty` = lmer(competence ~ 1 + domain*pdss_c + domain*pdss_c2 + 
                              (1 + pdss_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`comp model 1`, `comp model 2 age`, `comp model 3 age`, `comp model 4 age`, `comp model 5 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "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	model df	AIC	X2	X2 df	p
comp model 1	4	491.88	–	–	–
comp model 2 age	7	490.33	7.55	3	.056
comp model 3 age	8	488.18	4.14	1	.042
comp model 4 age	8	492.32	0	0	–
comp model 5 age	10	492.15	4.18	2	.124

puberty

anova(`comp model 1`, `comp model 2 puberty`, `comp model 3 puberty`, `comp model 4 puberty`, `comp model 5 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "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	model df	AIC	X2	X2 df	p
comp model 1	4	491.88	–	–	–
comp model 2 puberty	7	493.55	4.32	3	.228
comp model 3 puberty	8	494.34	1.22	1	.27
comp model 4 puberty	8	495.54	0	0	–
comp model 5 puberty	10	497.65	1.89	2	.389

summarize best fitting model

age

# summarize best fitting model
`comp model 3 age` %>%
  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", p)),
         term = gsub("\\(Intercept\\)", "intercept (academic, age 13)", term),
         term = gsub("domainsocial", "domain (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("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--"),
         SE = sprintf("%.3f", SE)) %>%
  mutate_if(is.numeric, round, 2) %>%
  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 (academic, age 13)	3.26 [3.16, 3.36]	0.052	62.97	227.38	< .001
domain (social)	-0.13 [-0.22, -0.05]	0.044	-2.99	223.05	.003
age	0.02 [0.00, 0.04]	0.011	1.96	65.85	.054
age2	-0.01 [-0.02, -0.00]	0.005	-2.03	317.50	.043

# summarize random effects
print(VarCorr(`comp model 3 age`), comp = c("Variance","Std.Dev."), digits = 4)

##  Groups    Name        Variance Std.Dev. Corr
##  subjectID (Intercept) 0.040830 0.20206      
##            age_c       0.001835 0.04284  0.28
##  Residual              0.175631 0.41908

visualize fitted

age

model_to_plot = `comp model 3 age`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
comp_model_data$expected = predict(model_to_plot, newdata = comp_model_data, re.form = reForm)
comp_model_data$expected_mean = predict(model_to_plot, newdata = comp_model_data, re.form = NA)

(comp_plot_age = comp_model_data %>%
  ggplot(aes(age, expected_mean, color = domain)) +
    geom_line(aes(age, competence, group = interaction(subjectID, domain)), alpha = .1) +
    geom_point(aes(age, competence, group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\nage (years)", y = "SPPC competence rating\n") +
    dcbw +
    theme(legend.position = "none"))

SPPC importance analysis

visualize raw

age

(imp_plot_raw_age = imp_model_data %>%
  filter(!is.na(importance)) %>%
  ggplot(aes(age, importance, color = domain, fill = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    scale_x_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
    labs(x = "\nage (years)", y = "mean importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(imp_plot_raw = imp_model_data %>%
  filter(!is.na(importance)) %>%
  ggplot(aes(pdss, importance, color = domain)) +
    geom_jitter(alpha = .1, width = .05, height = .05) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "mean importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(importance ~ subjectID, data = imp_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)

## [1] 0.24

age

`imp model 1` = lmer(importance ~ 1 + domain +
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 2 age` = lmer(importance ~ 1 + domain + age_c + 
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 3 age` = lmer(importance ~ 1 + domain*age_c + 
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`imp model 2 puberty` = lmer(importance ~ 1 + domain + pdss_c + 
                              (1 + domain | subjectID), 
                              data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 3 puberty` = lmer(importance ~ 1 + domain*pdss_c + 
                              (1 + domain | subjectID), 
                              data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`imp model 1`, `imp model 2 age`, `imp model 3 age`) %>%
  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
imp model 1	6	356.28	–	–	–
imp model 2 age	7	356.65	1.62	1	.203
imp model 3 age	8	356.69	1.97	1	.161

puberty

anova(`imp model 1`, `imp model 2 puberty`, `imp model 3 puberty`) %>%
  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
imp model 1	6	356.28	–	–	–
imp model 2 puberty	7	357.77	0.51	1	.476
imp model 3 puberty	8	359.27	0.5	1	.482

summarize best fiting model

`imp model 1` %>%
  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", p)),
         term = gsub("\\(Intercept\\)", "intercept (academic, age 13)", term),
         term = gsub("domainsocial", "domain (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("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  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 (academic, age 13)	3.35 [3.23, 3.48]	0.06	52.25	55.06	< .001
domain (social)	-1.01 [-1.21, -0.81]	0.10	-10.00	55.63	< .001

# summarize random effects
print(VarCorr(`imp model 1`), comp = c("Variance","Std.Dev."), digits = 4)

##  Groups    Name         Variance Std.Dev. Corr 
##  subjectID (Intercept)  0.04644  0.2155        
##            domainsocial 0.19984  0.4470   -0.24
##  Residual               0.29299  0.5413

visualize fitted

imp_model_data_na = imp_model_data %>%
  filter(!is.na(importance))

model_to_plot = `imp model 1`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
imp_model_data_na$expected_mean = predict(model_to_plot, newdata = filter(imp_model_data_na, !is.na(importance)), re.form = NA)

(imp_plot_age = imp_model_data_na %>%
  ggplot(aes(age, expected_mean, color = domain)) +
    geom_line(aes(age, importance, group = interaction(subjectID, domain)), alpha = .1) +
    geom_point(aes(age, importance, group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_x_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\nage (years)", y = "SPPC importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .14),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

SPPC competence importance plot

cowplot::plot_grid(comp_plot_age, imp_plot_age, ncol = 2, labels = c("A", "B"))

SPPC global self-worth

visualize raw

age

(sw_plot_raw_age = self_worth_model_data %>%
  ggplot(aes(age, self_worth)) +
    geom_point(alpha = .1, color = pal_self_other[2]) +
    geom_line(aes(group = interaction(subjectID)), alpha = .2, color = pal_self_other[2]) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .4,
                color = pal_self_other[2], fill = pal_self_other[2]) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\nage (years)", y = "mean global self-worth rating\n") +
    dcbw)

puberty

(sw_plot_raw = self_worth_model_data %>%
  ggplot(aes(pdss, self_worth, color = domain)) +
    geom_point(alpha = .1, color = pal_self_other[2]) +
    geom_line(aes(group = interaction(subjectID)), alpha = .2, color = pal_self_other[2]) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .4,
                color = pal_self_other[2], fill = pal_self_other[2]) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "mean self_worth rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(self_worth ~ subjectID, data = self_worth_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)

## [1] 0.7

age

`sw model 1` = lmer(self_worth ~ 1 + (1 | subjectID), 
                              data = self_worth_model_data, control = lmerControl(optimizer = "bobyqa"))
`sw model 2 age` = lmer(self_worth ~ 1 + age_c + 
                             (1 + age_c | subjectID), 
                             data = self_worth_model_data, control = lmerControl(optimizer = "bobyqa"))
`sw model 3 age` = lmer(self_worth ~ 1 + age_c + age_c2 + 
                             (1 + age_c | subjectID), 
                             data = self_worth_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`sw model 2 puberty` = lmer(self_worth ~ 1 + pdss_c + 
                             (1 + pdss_c | subjectID), 
                             data = self_worth_model_data, control = lmerControl(optimizer = "bobyqa"))
`sw model 3 puberty` = lmer(self_worth ~ 1 + pdss_c + pdss_c2 +
                             (1 + pdss_c | subjectID), 
                             data = self_worth_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`sw model 1`, `sw model 2 age`, `sw model 3 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "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	model df	AIC	X2	X2 df	p
sw model 1	3	210.90	–	–	–
sw model 2 age	6	209.63	7.28	3	.064
sw model 3 age	7	209.33	2.29	1	.13

puberty

anova(`sw model 1`, `sw model 2 puberty`, `sw model 3 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "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	model df	AIC	X2	X2 df	p
sw model 1	3	210.90	–	–	–
sw model 2 puberty	6	208.18	8.72	3	.033
sw model 3 puberty	7	209.56	0.62	1	.431

summarize best fitting model

puberty

# summarize best fitting model
`sw model 2 puberty` %>%
  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", p)),
         term = gsub("\\(Intercept\\)", "Intercept (PDS 3)", term),
         term = gsub("pdss_c", "puberty", 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("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--"),
         SE = sprintf("%.3f", SE)) %>%
  mutate_if(is.numeric, round, 2) %>%
  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 (PDS 3)	3.23 [3.15, 3.31]	0.041	78.35	85.14	< .001
puberty	-0.04 [-0.09, 0.00]	0.023	-1.93	62.45	.058

# summarize random effects
print(VarCorr(`sw model 2 puberty`), comp = c("Variance","Std.Dev."), digits = 4)

##  Groups    Name        Variance Std.Dev. Corr
##  subjectID (Intercept) 0.087621 0.29601      
##            pdss_c      0.007484 0.08651  0.35
##  Residual              0.098641 0.31407

visualize fitted

age

model_to_plot = `sw model 2 puberty`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
self_worth_model_data$expected = predict(model_to_plot, newdata = self_worth_model_data, re.form = reForm)
self_worth_model_data$expected_mean = predict(model_to_plot, newdata = self_worth_model_data, re.form = NA)

(sw_plot_pdss = self_worth_model_data %>%
  ggplot(aes(pdss, expected_mean)) +
    geom_line(aes(pdss, self_worth, group = subjectID), alpha = .1, color = pal_self_other[2]) +
    geom_point(aes(pdss, self_worth, group = subjectID), alpha = .4, color = pal_self_other[2]) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .4,
                color = pal_self_other[2], fill = pal_self_other[2]) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "SPPC global self-worth rating\n") +
    dcbw +
    theme(legend.position = "none"))

correlations

Repeated measures correlations between social and academic competence and importance for self-worth

competence

data_comp_corr = comp_model_data %>%
  select(-expected_mean, -expected) %>%
  spread(domain, competence) %>%
  mutate(subjectID = as.factor(subjectID))

rmcorr::rmcorr(subjectID, social, academic, dataset = data_comp_corr)

## 
## Repeated measures correlation
## 
## r
## 0.3100034
## 
## degrees of freedom
## 95
## 
## p-value
## 0.002000875
## 
## 95% confidence interval
## 0.1156885 0.4814649

importance

data_imp_corr = imp_model_data %>%
  spread(domain, importance) %>%
  mutate(subjectID = as.factor(subjectID))

rmcorr::rmcorr(subjectID, social, academic, dataset = data_imp_corr)

## 
## Repeated measures correlation
## 
## r
## 0.1159001
## 
## degrees of freedom
## 35
## 
## p-value
## 0.4945415
## 
## 95% confidence interval
## -0.2260783 0.4324442