index.qmd

---
title: "Online supplementary materials"
subtitle: "``Behavioral and neural responses to social rejection: Individual differences in developmental trajectories across childhood and adolescence'' by Mulder, Dobbelaar, and Achterberg (2024)"
author: 
  - name: Jeroen D. Mulder
    url: https://orcid.org/0000-0002-5553-0856
    affiliation: Utrecht University
    affiliation_url: https://www.uu.nl/staff/jdmulder
    orcid_id: 0000-0002-5553-0856
format:
  html:
    self-contained: true
    grid:
      margin-width: 350px

toc: true
toc-location: left
reference-location: margin
citation-location: document
bibliography: references.bib
---

```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
library(tidyverse)
library(semTools)
library(haven) # Read in .sav files
library(naniar) # Missing data handling
library(ggplot2) # Data visualization
```

# Introduction
This document contains the online supplementary materials to [Mulder, Dobbelaar, and Achterberg (2024)](https://doi.org/10.1016/j.dcn.2024.101365). It serves as a methodological and statistical rationale, and highlights the R code that was used to fit the models and explore the results. Questions regarding the analyses in the paper can be addressed to the first author, [Jeroen D. Mulder](https://orcid.org/0000-0002-5553-0856). 

First, the study's research aims are briefly summarized. Second, the data used in the study are introduced, in which the focus lies on those characteristics of the data that require due consideration in the subsequent statistical analysis. Third, the statistical models are introduced in detail. Fourth, the models are fitted to the data, and model convergence and fit are assessed. Finally, numerical results are interpreted and visualized. The complete R code (i.e., including code that was used for data manipulation and the creation of visualizations) can be found on [GitHub](https://github.com/JeroenDMulder/social-emotion-regulation). 

# Research aims 
This study has three main aims:

1. To describe developmental trajectories of behavioral (aggressive) and brain responses of negative social feedback throughout childhood and emerging adolescence,
2. to estimate the associations between development in behavioral aggression regulation and development in neural responses to social rejection, and
3. to estimate the relationship between development in social emotion regulation and social well-being in adolescence.

# The data
The data in this study come from the experimental longitudinal twin sample from the Leiden Consortium Individual Development (L-CID). Participants in the sample were invited to participate in fMRI sessions at three measurement occasions with approximately two-year intervals, termed wave 1, 3, and 5, respectively. During the fMRI sessions, the social network aggression task (SNAT) was used to measure behavioral (aggressive) responses to social rejection, and neural responses to social rejection. A sixth wave of data was collected in which participants' social well-being was measured. More detailed information about the sample used in this study is given in this study's preregistration by @achterberg_individual_2022, and in @dobbelaar_development_2022. For details on the Leiden Consortium Individual Development project, the reader is referred to @crone_neural_2020. 

## Behavioural responses
Aggressive behavioral responses to social rejection were measured during each fMRI session using the social network aggression task (SNAT). In short, during the SNAT, a participant would receive feedback from multiple peers on their own personal profile (which they created themselves specifically for the study). The feedback was either positive, neutral, or negative. Participants were told that they could press a button to send a loud noise blast to the peers in response to the feedback. The longer the button press, the longer, and louder, the noise blast would be. The length of the button press (in ms) was recorded, and taken as a proxy for behavioral aggressive response (maximum length of the button press was 3500 ms). In total, there were 20 neutral trials and 20 negative trials per fMRI session, resulting in three sets of 40 repeated behavioral responses for individuals who completed measurements at all three waves. 

Several characteristics of the behavioural response measurements warrant attention. First, the data have **a hierarchical structure**, with repeated feedback *trials* (level 1) nested in measurement *waves* (level 2), which are nested in *individuals* (level 3), which are nested in *families* (level 4). Second, there are **individually-varying times of observation** for the fMRI sessions (see the reported age range per measurement wave in Table 1 of the main paper). Hence, at any given wave, participants might have been at a different stage in their developmental process. Third, @fig-cens shows that a significant amount of behavioral response measurements have been censored from above, resulting from a time limit on the length of the button press.

```{r echo=FALSE, warning=FALSE, message=FALSE}
#| label: fig-cens
#| fig-cap: "Noise blast duration (seconds) per measurement occasion."
#| fig-width: 11
#| fig-height: 3
#| fig-column: page-right
#| warning: false
read.csv("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_SNAT_Be.csv") |> 
  naniar::replace_with_na_all(~ .x == -999) |>
  dplyr::mutate(
    noise1 = T1_noiseblast / 1000, 
    noise2 = T2_noiseblast / 1000,
    noise3 = T3_noiseblast / 1000
  ) |>
  dplyr::select(noise1, noise2, noise3) |>
  tidyr::pivot_longer(
    cols = noise1:noise3,
    names_to = "wave",
    names_prefix = "noise",
    names_transform = list(wave = as.integer),
    values_to = "noise"
  ) |>
  dplyr::mutate(
    censored = ifelse(noise == 3.500, 1, 0),
    wave_factor = factor(   
      wave, 
      levels = c("1", "2", "3"), 
      labels = c("Wave 1", "Wave 2", "Wave 3")
    )
  ) |>
  ggplot2::ggplot(aes(x = noise)) + 
    ggplot2::geom_histogram() +
    ggplot2::scale_x_continuous(breaks = c(0, .5, 1, 1.5, 2, 2.5, 3, 3.5)) + 
    ggplot2::facet_wrap(. ~ wave_factor, ncol = 3) +
    ggplot2::labs(x = "Behavioural response (ms)", y = "Frequency") + 
    ggplot2::theme_light() + 
    ggplot2::theme(text = element_text(size = 12))
```

Fourth, there is a fair amount of **missing data**. For the behavioral response measurements, this exists almost entirely of dropout. The `noise` columns in @fig-missing display the different missing data patterns across time for these measurements. 

```{r missing, echo=FALSE}
#| label: fig-missing
#| fig-cap: "Missing data patterns. The behavioural measurements are denoted by `noise`. `DLPFC`, `INS`, and `MPFC` are the neural responses in the dorsolateral prefrontal cortex, anterior insula, and medial prefrontal cortext, respectively."
#| fig-width: 8
#| fig-height: 4
#| warning: false

df_noise_wide <- read.csv("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_SNAT_Be.csv") |>
  select(SubjectID, T1_noiseblast, T2_noiseblast, T3_noiseblast) |>
  rename(id_orig = SubjectID) |>
  naniar::replace_with_na_all(~ .x == -999) |>
  group_by(id_orig) |>
  summarise(
    noise1 = ifelse(all(is.na(T1_noiseblast)), NA, 1),
    noise2 = ifelse(all(is.na(T2_noiseblast)), NA, 1),
    noise3 = ifelse(all(is.na(T3_noiseblast)), NA, 1)
  )
  
df_INS_wide <- haven::read_spss("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_MRI.sav") |>
  dplyr::select(SubjectID, ends_with("INS_Neg_Neut")) |>
  dplyr::rename(
    id_orig = SubjectID,
    AI1 = T1_INS_Neg_Neut,
    AI2 = T2_INS_Neg_Neut,
    AI3 = T3_INS_Neg_Neut
  )

df_MPFC_wide <- haven::read_spss("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_MRI.sav") |>
  dplyr::select(SubjectID, ends_with("MPFC_Neg_Neut")) |>
  dplyr::rename(
    id_orig = SubjectID,
    MPFC1 = T1_MPFC_Neg_Neut,
    MPFC2 = T2_MPFC_Neg_Neut,
    MPFC3 = T3_MPFC_Neg_Neut
  )

df_wide <- haven::read_spss("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_MRI.sav") |>
  dplyr::select(SubjectID, ends_with("DLPFC_Neg_Neut")) |>
  dplyr::rename(
    id_orig = SubjectID,
    DLPFC1 = T1_DLPFC_Neg_Neut,
    DLPFC2 = T2_DLPFC_Neg_Neut,
    DLPFC3 = T3_DLPFC_Neg_Neut
  ) |>
  dplyr::right_join(df_INS_wide, by = "id_orig") |>
  dplyr::right_join(df_MPFC_wide, by = "id_orig") |>
  dplyr::right_join(df_noise_wide, by = "id_orig") |>
  select(-id_orig)

rm("df_noise_wide", "df_INS_wide", "df_MPFC_wide")

ggmice::plot_pattern(
  data = df_wide,
  square = FALSE,
  caption = FALSE
) +
  ggplot2::theme(text = element_text(size = 12)) +
  ggplot2::scale_x_discrete(
    limits = c("DLPFC1", "INS1", "MPFC1", "noise1",
               "DLPFC2", "INS2", "MPFC2", "noise2",
               "DLPFC3", "INS3", "MPFC3", "noise3")
  )
```

## Neural responses
In addition to behavioral responses after social feedback, neural responses after social feedback were also measured using fMRI scans during the SNAT, specifically in the anterior insula `AI`, the medial prefrontal cortex `MPFC`, and the dorsolateral prefrontal corect `DLPFC`. Like the behavioral response measure, the neural responses also have a hierarchical structure. However, preprocessing of the fMRI data resulted in measurements that are contrasts between average brain response for neutral feedback versus average brain response for negative feedback. Hence, the neural measurements have a three level structure compared to the behavioral measurements (i.e., there is no neural data at the trial level, level 1). Furthermore, there is missing data in the neural measurements, as visualized in @fig-missing. This is a combination of dropout and measurement error (e.g., when children move too much to get an accurate fMRI scan; Table 1 in the main paper contains an overview of reasons for exclusion frequencies of fMRI measurements per wave). Finally, as the behavioral and neural measurements were taken simultaneously, there are **individually-varying times of observation** to also need to be taken into account for these measurements. 

The neural responses (i.e., contrast) across time are plotted in the spaghetti plots below. These plots are useful as an initial exploration of individual development in brain response. The thick black line represent the average linear development in brain response, and colored lines represent individual brain responses across time. From the plot it becomes clear that there are many differences in how individuals develop over time, and that on average development appears to be only minor.

```{r echo=FALSE, warning=FALSE, message=FALSE}
#| label: fig-spaghetti
#| fig-cap: "Observed neural responses across age."
#| fig-width: 11
#| fig-height: 4
#| fig-column: page-right
#| warning: false
df_age <- read.table(
  "C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_Age.txt",
  header = TRUE
) |>
  naniar::replace_with_na_all(~ .x == 999) |>
  dplyr::mutate(
    id_orig = SubjectID,
    id = SubjectID %% 100,
    age1 = Age_wave1,
    age2 = Age_wave3,
    age3 = Age_wave5,
  ) |>
  dplyr::select(id, id_orig, age1, age2, age3) |>
  tidyr::pivot_longer(
    cols = age1:age3,
    names_to = "wave",
    names_prefix = "age",
    names_transform = list(wave = as.integer),
    values_to = "age"
  ) |>
  as.data.frame()

df_INS <- read_spss("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_MRI.sav") |>
  select(SubjectID, ends_with("INS_Neg_Neut")) |>
  dplyr::rename(id_orig = SubjectID,
                AI1 = T1_INS_Neg_Neut,
                AI2 = T2_INS_Neg_Neut,
                AI3 = T3_INS_Neg_Neut) |>
  tidyr::pivot_longer(
    cols = !id_orig,
    names_to = "wave",
    names_prefix = "AI", 
    values_to = "AI",
    names_transform = list(wave = as.integer)
  ) 

df_MPFC <- read_spss("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_MRI.sav") |>
  select(SubjectID, ends_with("MPFC_Neg_Neut")) |>
  dplyr::rename(id_orig = SubjectID,
                MPFC1 = T1_MPFC_Neg_Neut,
                MPFC2 = T2_MPFC_Neg_Neut,
                MPFC3 = T3_MPFC_Neg_Neut) |>
  tidyr::pivot_longer(
    cols = !id_orig,
    names_to = "wave",
    names_prefix = "MPFC", 
    values_to = "MPFC",
    names_transform = list(wave = as.integer)
  ) 

df_brain <- read_spss("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\data\\202204_Research_EmotionRegulation_MRI.sav") |>
  select(SubjectID, ends_with("DLPFC_Neg_Neut")) |>
  dplyr::rename(id_orig = SubjectID,
                DLPFC1 = T1_DLPFC_Neg_Neut,
                DLPFC2 = T2_DLPFC_Neg_Neut,
                DLPFC3 = T3_DLPFC_Neg_Neut) |>
  tidyr::pivot_longer(
    cols = !id_orig,
    names_to = "wave",
    names_prefix = "DLPFC", 
    values_to = "DLPFC",
    names_transform = list(wave = as.integer)
  ) |> 
  dplyr::right_join(df_INS, by = c("id_orig", "wave")) |>
  dplyr::right_join(df_MPFC, by = c("id_orig", "wave"))|>
  dplyr::right_join(df_age, by = c("id_orig", "wave")) |>
  tidyr::pivot_longer(
    cols = c("AI", "MPFC", "DLPFC"),
    names_to = "RoI",
    values_to = "response"
  )
  
ggplot(df_brain, aes(x = age, y = response)) +
  geom_line(aes(color = factor(id_orig)), size = .2) + 
    geom_smooth(aes(group = 1), 
    size = 2, 
    method = "lm", 
    se = FALSE,
    colour = "black") + 
  labs(x = "Age (years)", y = "Observed DLPFC response") + 
  ggplot2::facet_wrap(. ~ RoI, ncol = 3) + 
  theme_light() + 
  theme(text = element_text(size = 12), legend.position = "none")
```

## Social well-being
At wave 6, a social well-being questionnaire was filled in, consisting of 5 subscales: Ten items from the Adolescent Wellbeing Paradigm (AWP), ten items from the World Health Organization Quality of Life Scale (WHOQoL), and three subscales (each five items) from the Harter’s Self-Perception Profile for Adolescents (SPPA), specifically the subscales Social Competence (SC), Close Friendships (CF) and Global Self-worth (GS). All items were answered on a four-point Likert scale, with low scores indicating low social well-being and high scores indicating high social well-being. Instructions in each of the subscale manuals were followed for the handling of missing data and scoring of subscale scores. As such simple mean scores were computed per subscale and used for further statistical analyses. 

# Statistical models {#sec-models}
To explore Aim 1, the behavioral and neural responses were analyzed using Bayesian multilevel growth curve models, and specified using the R package `brms` [@burkner_brms_2017] in R version 4.2.2 [@r_core_team_r_2019]. The Bayesian estimation framework was chosen because it can easily accommodate censoring by integrating out censored values [@gelman_bayesian_2013]. Moreover, missing data can be imputed as part of the model fitting process under the assumption of missing at random (MAR) [@burkner_brms_2017]. The multilevel framework was chosen because it treats the time of observation as data (rather than as a parameter, as in a latent growth curve model), and can therefore more easily handle individually-varying times of observations [@mehta_putting_2000]. For Aims 2 and 3, a single structural equation model was specified in which social well-being was predicted from the estimated individual growth components extracted from the Bayesian multilevel growth curve model. 

## Intercept-only models 
To assess the proportion of observed variance that can be explained by each of the levels (expressed in terms of an intraclass correlation), intercept-only models were run first. The tabs below contain R code for estimating each of these intercept-only models using `brms`. The results are displayed in Table 2 of the main paper. 

::: {.panel-tabset}
### `INS`

```{.r eval=FALSE}
m_INS_ICC <- brm(
  formula = INS ~ 1 + (1 | fam/id),
  data = df_brain,
  family = gaussian(),
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random", 
  cores = 2,
  thin = 4,
  seed = 20221005
)

ICC_id <- "sd_fam:id__Intercept^2 / (sd_fam__Intercept^2 + sd_fam:id__Intercept^2 + sigma^2) > 0"
ICC_fam <- "sd_fam__Intercept^2 / (sd_fam__Intercept^2 + sd_fam:id__Intercept^2 + sigma^2) > 0"

hypothesis(m_INS_ICC, hypothesis = ICC_id, class = NULL, robust = FALSE)
hypothesis(m_INS_ICC, hypothesis = ICC_fam, class = NULL, robust = FALSE)
```

### `MPFC`

```{.r eval=FALSE}
m_MPFC_ICC <- brm(
  formula = MPFC ~ 1 + (1 | fam/id),
  data = df_brain,
  family = gaussian(),
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random", 
  cores = 2,
  thin = 4,
  seed = 20221005,
  control = list(adapt_delta = 0.9)
)

hypothesis(m_MPFC_ICC, hypothesis = ICC_id, class = NULL, robust = FALSE)
hypothesis(m_MPFC_ICC, hypothesis = ICC_fam, class = NULL, robust = FALSE)
```

### `DLPFC`

```{.r eval=FALSE}
m_DLPFC_ICC <- brm(
  formula = DLPFC ~ 1 + (1 | fam/id),
  data = df_brain,
  family = gaussian(),
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random", 
  cores = 2,
  thin = 4,
  seed = 20221005,
  control = list(adapt_delta = 0.9)
)

hypothesis(m_DLPFC_ICC, hypothesis = ICC_id, class = NULL, robust = FALSE)
hypothesis(m_DLPFC_ICC, hypothesis = ICC_fam, class = NULL, robust = FALSE)
```

### `noise`

```{.r eval=FALSE}
m_noise_gaus_ICC <- brm(
  formula = noise | cens(censored) ~ 1 + (1 | fam/id/wave),
  data = df_noise,
  family = gaussian(),
  warmup = 2000,
  iter = 6000,
  chains = 2,
  init = "random",
  cores = 2,
  thin = 6,
  seed = 20221005,
  control = list(adapt_delta = 0.9)
)

ICC_noise_fam <- "sd_fam__Intercept^2 / (sd_fam__Intercept^2 + sd_fam:id__Intercept^2 + sd_fam:id:wave__Intercept^2 + sigma^2) > 0"
ICC_noise_id <- "sd_fam:id__Intercept^2 / (sd_fam__Intercept^2 + sd_fam:id__Intercept^2 + sd_fam:id:wave__Intercept^2 + sigma^2) > 0"
ICC_noise_wave <- "sd_fam:id:wave__Intercept^2 / (sd_fam__Intercept^2 + sd_fam:id__Intercept^2 + sd_fam:id:wave__Intercept^2 + sigma^2) > 0"

hypothesis(m_noise_gaus_ICC, hypothesis = ICC_noise_fam, class = NULL, robust = TRUE)
hypothesis(m_noise_gaus_ICC, hypothesis = ICC_noise_id, class = NULL, robust = TRUE)
hypothesis(m_noise_gaus_ICC, hypothesis = ICC_noise_wave, class = NULL, robust = TRUE)
```

:::

## Bayesian multilevel growth curve model (Aim 1)
The growth curve models were build up as follows. Let $noise_{fiwt}$ be the noise blast duration (in seconds) for an individual $i$ from family $f$, at measurement wave $w$ and trial $t$. At level 1, the model is given by
$$
noise_{fiwt} = \beta_{0fiw} + \beta_{1fiw} \; x_{fiwt} + e_{fiwt}. 
$$
$x_{fiwt}$ is the condition (type of feedback, with 0 = neutral, negative = 1), $\beta_{0fiw}$ is a random mean representing noise blast duration for neutral feedback for individual $i$ from family $f$ at wave $w$, $\beta_{1fiw}$ is a random slope representing the difference between the neutral and negative condition for individual $i$ from family $f$ at wave $w$, and $e_{fiwt}$ represents the residual. 

At level 2, the model is specified as
$$\beta_{0fiw} = \gamma_{00fi} + u_{0fiw},$$
$$\beta_{1fiw} = \gamma_{10fi} + \gamma_{11fi} \; a_{fiw} + \gamma_{12fi} \; a_{fiw}^{2} + u_{1fiw},$$
where $\gamma_{00fi}$ is a random intercept representing noise blast for the neutral condition, and $u_{0fiw}$ is the random effect with variance $Var[u_{0fiw}]$, representing a family-, individual-, and wave-specific deviation from $\gamma_{00fi}$. To describe development in response over time, $\beta_{1fiw}$ is predicted from age $a_{fiw}$, and age-squared $a^{2}_{fiw}$. The predictor is grand-mean centered at 0 to prevent issues with multicollinearity, such that $a_{fiw} = 0$ now represents an age of approximately 9 years and 9 months. The coefficients $\gamma_{10fi}$, $\gamma_{11fi}$, and $\gamma_{12fi}$ are the growth components of interest, describing the mean response, linear development in response, and quadratic development in response of individual $i$ from family $f$ at 9 years and 9 months of age, respectively. $u_{1fiw}$ is a random component with variance $Var[u_{1fiw}]$, representing a family-, individual-, and wave-specific deviations from the expected $\beta_{1fiw}$.

At level 3, the model is specified as 
$$\gamma_{00fi} = \delta_{0} + v_{0fi},$$
$$\gamma_{10fi} = \delta_{1f} + v_{1fi},$$
$$\gamma_{11fi} = \delta_{2f} + v_{2fi},$$
$$\gamma_{12fi} = \delta_{3f} + v_{3fi},$$
where $\delta_{0f}$ represents the mean noise blast of family $f$ for the neutral condition, and $v_{0fi}$ represents the individual-specific deviation from $\delta_{00}$ with variance $Var[v_{0fi}]$. $\delta_{1f}$, $\delta_{2f}$, and $\delta_{3f}$ represent the family means of the growth components described on level 2, and $v_{1fi}$, $v_{2fi}$, and $v_{3fi}$ represent individual-specific (i.e., within-family) deviations from these means with variances $Var[v_{1fi}]$, $Var[v_{2fi}]$, and $Var[v_{3fi}]$. 

Finally, at level 4, the model is given by
$$\delta_{0f} = \theta_{0} + z_{0f},$$
$$\delta_{1f} = \theta_{1} + z_{1f},$$
$$\delta_{2f} = \theta_{2} + z_{2f},$$
$$\delta_{3f} = \theta_{3} + z_{3f}.$$
where $\theta_{0}$ is the grand mean (across families, individuals, waves, and trials) of noise blast for the neutral condition. $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ are the grand means of the growth components. $z_{1f}$, $z_{2f}$, and $z_{3f}$ are family-specific differences herein, with variances $Var[z_{1f}]$, $Var[z_{2f}]$, and $Var[z_{3f}]$. In total, 15 parameters are estimated in this model, namely 4 fixed effects and 4 random effect variances at level 4, 4 random effect variances at level 3, 2 random effect variances at level 2, and 1 residual variance at level 1. 

Missing data were imputed as part of model fitting. Due to the setup of the study (with a planned two-year interval between subsequent measurement waves), missing values for age $a^{mis}_{fiw}$ could be predicted from the individual's age at the first measurement wave `ageW1c` and the wave number. Missing values for age-squared $a^{2, mis}_{fiw}$ could then be created simply by squaring $a^{mis}_{fiw}$. This resulted in the exclusion of a single family (two individuals) from the sample, who had missing values at their first measurement wave. For the outcome, missing values were predicted from feedback condition $X$, age $A$, and age-squared $A^{2}$. For this procedure to work well, missing at random (MAR) has to be assumed, implying that the degree of missingness only depends on the observed data. While, the plausibility of this assumption can be questioned, it is a weaker assumption than missing completely at random (MCAR), which would be assumed when the missing values were simply deleted from the data [@gelman_bayesian_2013]. 

For the brain data, mean response is already available from the preprocessing of the MRI data. Therefore, in the model above, the trial level is removed, as well as the equation for $\beta_{0fiw}$ at level 2. The interpretation of the remaining parameters, and the procedure used to impute for missing values, are equivalent. 

R code for model fitting is printed below, for each outcome in a separate tab. `brms` allows for models to be specified using `lme4`-like syntax. The `mi()` wrapper around variables denotes that missing values are imputed. For each outcome, `iter = 4000` values were sampled from the posterior distribution, of which the first `warmup = 2000` were discarded. `chains = 2` chains were used, and convergence was assessed using the split-$\hat{R}$ value [@gelman_bayesian_2013]. Only every `thin = 4`$^{th}$ simulation draw was kept to keep model size manageable, resulting in a total of 1000 posterior simulation draws. `Stan`'s default uninformative priors were used because no independent prior information was available on the development of behavioral and brain emotion regulation as measured using the SNAT. 

::: {.panel-tabset}

### `noise` 
Noise blast measurements have a lower limit of 0, and are restricted from above at 3500 ms. In theory, this information could be included in the model in the form of a bounded prior on the range of values that the fixed effects $\theta$ could take on. However, the `brms` manual states that such a procedure is rarely useful, and hence default uninformative priors were used. The `cens()` function was used to denote which values were censored. Two models where fitted. The first assumed that the outcome was Gaussian distributed. Second, given the complexity of this 4-level model, and the restricted range of `noise`, a model assuming a *t*-distributed outcome was fitted as well to allow for more robust inference [@gelman_bayesian_2013]. 

```{r}
#| label: model-noise
#| eval: FALSE

# Gaussian response distribution
f_noise_gaus <- 
  bf(noise | cens(censored) + mi() ~ 1 + con + con:mi(agec) + con:mi(agec2) + 
       (1 + con | fam/id/wave) + (-1 + con:mi(agec) + con:mi(agec2) | fam/id)) + 
  gaussian() + 
  bf(agec | mi() ~ -1 + ageW1c + wave0) +
  bf(agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02)

m_noise_gaus <- brm(f_noise_gaus,
  data = df_noise,
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random",
  cores = 2,
  thin = 4,
  seed = 20221005,
  control = list(adapt_delta = 0.85)
)

# Student-t response distribution
f_noise_stud <- 
  bf(noise | cens(censored) + mi() ~ 1 + con + con:mi(agec) + con:mi(agec2) + 
       (1 + con | fam/id/wave) + (-1 + con:mi(agec) + con:mi(agec2) | fam/id)) + 
  student() + 
  bf(agec | mi() ~ -1 + ageW1c + wave0) +
  bf(agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02)

m_noise_stud <- brm(f_noise_stud,
  data = df_noise,
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random",
  cores = 2,
  thin = 4,
  seed = 20221005,
  control = list(adapt_delta = 0.88)
)
```

### `DLPFC` 

```{r DLPFC, eval=FALSE}
f_DLPFC_3L <- 
  bf(DLPFC | mi() ~ 1 + mi(agec) + mi(agec2) + (1 + mi(agec) + mi(agec2) | fam/id)) + 
  gaussian() +
  bf(agec | mi() ~ -1 + ageW1c + wave0) +
  bf(agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02)

m_DLPFC_3L <- brm(f_DLPFC_3L,
  data = df_brain,
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random", 
  cores = 2,
  thin = 4,
  seed = 20221005
)
```

### `AI`

The term `AI` in the R code below actually refers to the anterior insula `AI`. 

```{r INS, eval=FALSE}
f_INS_3L <- 
  bf(INS | mi() ~ 1 + mi(agec) + mi(agec2) + (1 + mi(agec) + mi(agec2) | fam/id)) + 
  gaussian() +
  bf(agec | mi() ~ -1 + ageW1c + wave0) +
  bf(agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02)

m_INS_3L <- brm(f_INS_3L,
  data = df_brain,
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random", 
  cores = 2,
  thin = 4,
  seed = 20221005
)
```

### `MPFC` 

```{r MPFC, eval=FALSE}
f_MPFC_3L <- 
  bf(MPFC | mi() ~ 1 + mi(agec) + mi(agec2) + (1 + mi(agec) + mi(agec2) | fam/id)) + 
  gaussian() + 
  bf(agec | mi() ~ -1 + ageW1c + wave0) +
  bf(agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02)

m_MPFC_3L <- brm(f_MPFC_3L,
  data = df_brain,
  warmup = 2000,
  iter = 4000,
  chains = 2,
  init = "random", 
  cores = 2,
  thin = 4,
  seed = 20221005
)
```

:::

### Deviations from preregistration
The preregistration by @achterberg_individual_2022 stated that correlations between the random coefficients at each level would be estimated. However, due to low level-specific sample sizes (i.e., three measurements per individual at level 2, two individuals per family at level 3), this was not statistically feasible. Therefore, these coefficients were excluded from the model. 

## Structural equation model (Aims 2 and 3)
The analyses described here are related to both Aims 2 and 3. For Aim 2, associations between the growth components of the four behavioral and brain outcomes were investigated. For Aim 3, estimated individual-level growth components were used to predict social well-being at wave 6. In this section, I first discus how growth components for each individual are extracted from the growth curve models. Second, a one-factor confirmatory factor analysis model is introduced to investigate the dimensionality of the five social well-being subscales. Third, a structural equation model is introduced to estimate correlations between the growth components, and to predict social well-being during adolescence.

### Creating individual-level growth components
The results from the fitted Bayesian multilevel growth curve models are based on 1000 draws of the poster distribution for each parameter, which includes both fixed effects and random effects. Therefore, at each iteration, the sampled fixed effects and individually drawn random effects can be extracted from the model and used compute personalized growth components. Specifically, for an individual $i$ from family $f$, and given a specific iteration, the growth components are computed as
$$\gamma_{10fi} = \theta_{1} + z_{1f} + v_{1if},$$
$$\gamma_{11fi} = \theta_{2} + z_{2f} + v_{2if},$$
$$\gamma_{12fi} = \theta_{3} + z_{3f} + v_{3if}.$$
Computing these three parameters, for each of the four outcomes, and for each iteration, results in 1000 data sets of $4 \times 3 = 12$ personalized growth components. These serve as input for predicting social well-being. R code for these computations is too long to display here, but is available at [GitHub](https://github.com/JeroenDMulder/social-emotion-regulation). 

### Dimensionality of social well-being
Social well-being data consists of 35 items, belonging to one of five subscales. The data were preprocessed by computing factor scores for each subscale by taking the mean of the respective items (as instructed in the subscale manuals). Using the R package *lavaan* (version 0.6.12), a one-factor confirmatory factor analysis model was fitted with the subscales scores as indicators [@rosseel_lavaan_2012].

```{r}
#| label: cfa-wellbeing
#| eval: FALSE
# Confirmatory factor analysis social well-being using lavaan.
m_cfa_well <- '
  # 1-factor model for social well-being
  well =~ 1*co + who + sc + cf + gs
'
```

### Multivariate regression/prediction model
Structural equation modeling was used to estimate the relations between the development behavioral and brain responses, as well as later social well-being. A multiple regression model was specified with the social well-being subscales as the outcome, and the growth components for each outcome as predictors. The predictors are allowed to covary freely with each other to investigate the associations between development in response of behavioral and brain outcomes. The estimated regression coefficients then represent the relationship between development in social emotion regulation and social well-being in adolescence. `lavaan` model syntax can be found below.

```{r}
#| label: lavaan
#| eval: FALSE
m_regression <- '
  # Regression
  co ~ gamma10_INS + gamma11_INS + gamma12_INS + gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  who ~ gamma10_INS + gamma11_INS + gamma12_INS + gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  sc ~ gamma10_INS + gamma11_INS + gamma12_INS + gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  cf ~ gamma10_INS + gamma11_INS + gamma12_INS + gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  gs ~ gamma10_INS + gamma11_INS + gamma12_INS + gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  
  # Variances 
  gamma10_INS ~~ gamma10_INS
  gamma11_INS ~~ gamma11_INS
  gamma12_INS ~~ gamma12_INS
  
  gamma10_DLPFC ~~ gamma10_DLPFC
  gamma11_DLPFC ~~ gamma11_DLPFC
  gamma12_DLPFC ~~ gamma12_DLPFC
  
  gamma10_MPFC ~~ gamma10_MPFC
  gamma11_MPFC ~~ gamma11_MPFC
  gamma12_MPFC ~~ gamma12_MPFC
  
  gamma10_noise ~~ gamma10_noise
  gamma11_noise ~~ gamma11_noise
  gamma12_noise ~~ gamma12_noise
  
  # Residual variances
  co ~~ co 
  who ~~ who 
  sc ~~ sc
  cf ~~ cf
  gs ~~ gs
  
  # Covariances
  gamma10_INS ~~ gamma11_INS + gamma12_INS + gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise 
  gamma11_INS ~~ gamma12_INS + gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise 
  gamma12_INS ~~ gamma10_DLPFC + gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise 
  
  gamma10_DLPFC ~~ gamma11_DLPFC + gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise 
  gamma11_DLPFC ~~ gamma12_DLPFC + gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  gamma12_DLPFC ~~ gamma10_MPFC + gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  
  gamma10_MPFC ~~ gamma11_MPFC + gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise 
  gamma11_MPFC ~~ gamma12_MPFC + gamma10_noise + gamma11_noise + gamma12_noise
  gamma12_MPFC ~~ gamma10_noise + gamma11_noise + gamma12_noise
  
  gamma10_noise ~~ gamma11_noise + gamma12_noise
  gamma11_noise ~~ gamma12_noise
  
  # Residual covariances
  co ~~ who + sc + cf + gs
  who ~~ sc + cf + gs
  sc ~~ cf + gs
  cf ~~ gs
' 
```

To take the uncertainty of the individual growth component estimates from the Bayesian multilevel growth curve models into account, the SEM model was fitted to each of the 1000 data sets using the R package `semTools` [@jorgensen_semtools_2022]. Point and standard error estimates were pooled according to rules by @rubin_multiple_1987. 

# Results

## Development of behavioral and brain responses (Aim 1)
In this section, model fit is evaluated using a posterior predictive check, and the results are presented. The results are organized below per outcome separately. 

::: {.panel-tabset}

### `noise` 

The posterior predictive checks (ppc) of two models for `noise` are shown in @fig-noise-ppc: The left panel is a ppc-plot from a model assuming a Gaussian-distributed outcome; the right panel is a ppc-plot from a model assuming a *Student t*-distributed outcome. The comparison shows that using a *t*-distribution, the model is better able to capture the kurtosis in the observed density, and that the modes of the posterior draws better align with the peak of the observed data. Therefore, the model assuming a *t*-distributed outcome is selected for further interpretation. 

```{r}
#| label: fig-noise-ppc
#| echo: FALSE
#| fig-cap: "Posterior predictive checks (ppc) for the behavioural (aggresive) response (noise). The model underlying the left plot assumes a normally-distributed outcome; the model underlying the right plot assumes a Student t-distributed outcome."
#| fig-width: 10
#| fig-height: 5
#| warning: false
knitr::include_graphics("figures\\ppc-noise.png")
```
The right panel of @fig-noise-ppc shows that the model predicts numerous values higher than 3.5. To assess how well the proportion of *predicted* values higher than 3.5 represents the proportion of *observed* censored values, a third posterior predictive check was performed. @fig-noise-ppc-rcens shows that the model underestimates the proportion of noise blasts longer 3.5 seconds. However, in absolute sense, this misfit is small enough to move forward with this model. 

```{r}
#| label: fig-noise-ppc-rcens
#| echo: FALSE
#| fig-cap: "Posterior predictive checks (ppc) for the proportion right censored outcomes for behavioural (aggresive) response. A Student t-distributed outcome was assumed."
#| fig-width: 10
#| fig-height: 5
#| warning: false
knitr::include_graphics(".\\figures\\ppc-noise-rcens.png")
```

Results for the `noise` model (with Student *t*-distributed response) are shown below. The fixed effects estimated $\hat{\theta}_{1}$, $\hat{\theta}_{2}$, and $\hat{\theta}_{3}$ can be found under `Population-Level Effects` as `noise_con`, `noise_con:miagec`, and `noise_con:miagec2`, respectively. Standard deviation estimates of random effects at the family-level (i.e., $z$ parameters) are reported under `~fam`, standard deviation estimations of random effects at the individual-level (i.e., $v$ parameters) are reported under `~fam:id`, and standard deviation estimates of random effecst at the wave-level (i.e., $u$ parameters) are reported under `fam:id:wave`.

```{r}
#| label: results-noise
#| eval: TRUE
#| echo: FALSE
#| fig-column: page-right
load("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\brms\\m_noise_stud.Rds")
summary(m_noise_stud)
```

The $\hat{R}$-value for the random within-family effect of the mean noise response $SD(v_{1})$ is not below the recommended value 1.05, indicating that the posterior sampling chains have not converged sufficiently yet to a stable solution. Moreover, the bulk effective sample size `Bulk_ESS` for this parameter, as well as for the random between-family effect of the mean noise response $SD(z_{1})$, is (far) below the recommended number of effective samples (100 times the number of chains). This speaks to the complexity of the model, especially in relation to the relatively low sample sizes at the family- and individual-level. 

The fixed effects show that there is 95% certainty that expected response in `noise` at 9 years and 9 months of age $\theta_{1}$ lies between $1.31$ and $1.49$ seconds. The random effects imply that there are significant differences *between* individuals (within families) herein, $SD(v_{1}) = [0.03 - 0.42]$, as well as significant between-family differences, $SD(z_{1}) = [0.15 - 0.43]$. Linear development of `noise` at mean age $\theta_{2}$ lies between  $-0.05$ and $0.02$. The random effect estimates suggests slight between family differences herein, $SD(z_{2}) = [0.01 - 0.16]$. The quadratic development of `noise` at mean age $\theta_{3}$ is estimated to lie between $-0.06$ and $-0.03$, and no significant differences between individuals, $SD(v_{3}) = [0.00 - 0.03]$, or between families, $SD(z_{3}) = [0.00 - 0.05]$, have been found. 

```{r}
#| label: fig-noise-predicted
#| echo: FALSE
#| fig-cap: "Model predicted development in behavioral response."
#| fig-width: 10
#| fig-height: 5
#| warning: false
knitr::include_graphics(".\\figures\\plot_noise_predicted.png")
```

### `DLPFC`

Results for the `DLPFC` model are shown below. 

```{r}
#| label: results-DLPFC
#| eval: TRUE
#| echo: FALSE
#| fig-column: page-right
load("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\brms\\m_DLPFC_3L.Rds")
summary(m_DLPFC_3L)
```

The $\hat{R}$-values are below the recommended 1.05, indicating adequate mixing of the posterior sampling chains. The posterior predictive checks in @fig-DLPFC-ppc below imply that the model shows adequate fit for both the imputed values of age and age$^{2}$, as well as for the response variable `DLPFC`.  

```{r}
#| label: fig-DLPFC-ppc
#| echo: FALSE
#| fig-cap: "Posterior predictive check for a) the DLPFC response, b) imputed values of age, and c) imputed values of age-squared. Dark blue lines represents observed data. Light blue lines represent model predicted (replicated) values."
#| warning: false
knitr::include_graphics("figures\\ppc-DLPFC.png")
```

The fixed effects show that there is 95% certainty that expected response in `DLPFC` at 9 years and 9 months of age $\theta_{1}$ lies between $-0.65$ and $-0.09$. The random effects imply that there are significant *between* families differences herein, $SD(z_{1}) = [0.01 - 0.84]$, as well as significant *within* family differences, $SD(v_{1}) = [0.02 - 0.97]$. Linear development of `DLPFC` at mean age $\theta_{2}$ lies between  $-0.14$ and $0.11$. Again, the random effect estimates suggests differences between families herein, $SD(z_{2}) = [0.01 - 0.43]$, as well as within families, $SD(v_{2}) = [0.02 - 0.53]$. The quadratic development of `DLPFC` at mean age $\theta_{3}$ is estimated to lie between $-0.02$ and $0.11$, and no significant differences between families, $SD(z_{3}) = [0.00 - 0.16]$, or within families, $SD(v_{3}) = [0.00 - 0.17]$, have been found. 

These results are visualized in @fig-DLPFC-predicted. The thick black line represents predicted development in `DLPFC` response over time based on the fixed effects estimates. The lighter black lines denote the uncertainty around this estimated development. The red dotted line denotes the mean age of the sample. 

```{r}
#| label: fig-DLPFC-predicted
#| echo: FALSE
#| fig-cap: "Model predicted development in DLPFC response."
#| warning: false
knitr::include_graphics(".\\figures\\plot_DLPFC_predicted.png")
```

### `AI`

Results for the `AI` model are shown below. Note that the term `INS` here refers to the anterior insula `AI`. 

```{r}
#| label: results-AI
#| eval: TRUE
#| echo: FALSE
#| fig-column: page-right
load("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\brms\\m_INS_3L.Rds")
summary(m_INS_3L)
```

The estimation procedure resulted in 1 warning about a single divergent transition after the warm-up phase. Divergences can occur when the sampling procedure has difficulties exploring the target posterior distribution, which is an issue when the curvature of the posterior is highly variable. However, because there is only a single divergence, and there are good $\hat{R}$ values, and effective sample sizes for the bulk and tail of the distribution, there is enough confidence in the resulting posterior to move forward. 

The posterior predictive checks in @fig-AI-ppc below imply that the model shows adequate fit for both the imputed values of age and age$^{2}$, as well as for the response variable `AI`.  

```{r}
#| label: fig-AI-ppc
#| echo: FALSE
#| fig-cap: "Posterior predictive check for a) the AI response, b) imputed values of age, and c) imputed values of age-squared."
#| warning: false
knitr::include_graphics("figures\\ppc-INS.png")
```

The fixed effects show that there is 95% certainty that expected response in `AI` at 9 years and 9 months of age $\theta_{1}$ lies between $0.42$ and $1.03$. The random effects imply that there are significant differences *between* individuals within families herein, $SD(z_{1}) = [0.01 - 0.81]$, as well as significant *within* family differences, $SD(v_{1}) = [0.02 - 1.01]$. Linear development of `AI` at mean age $\theta_{2}$ lies between  $-0.29$ and $-0.02$. Again, the random effect estimates suggests differences between families herein, $SD(z_{2}) = [0.01 - 0.35]$, as well as within families, $SD(v_{2}) = [0.03 - 0.66]$. The quadratic development of `AI` at mean age $\theta_{3}$ is estimated to lie between $-0.01$ and $0.13$, and no significant differences between families, $SD(z_{3}) = [0.00 - 0.14]$, or within families, $SD(v_{3}) = [0.00 - 0.17]$, have been found. 

These results are visualized in @fig-AI-predicted. The thick black line represents predicted development in `AI` response over time based on the fixed effects estimates. The lighter black lines denote the uncertainty around this estimated development. The red dotted line denotes the mean age of the sample. 

```{r}
#| label: fig-AI-predicted
#| echo: FALSE
#| fig-cap: "Model predicted development in AI response."
#| warning: false
knitr::include_graphics(".\\figures\\plot_INS_predicted.png")
```

### `MPFC` 

Results for the `MPFC` model are shown below:

```{r}
#| label: results-MPFC
#| eval: TRUE
#| echo: FALSE
#| fig-column: page-right
load("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\brms\\m_MPFC_3L.Rds")
summary(m_MPFC_3L)
```

The estimation procedure resulted in 1 warning about a single divergent transition after the warm-up phase. Divergences can occur when the sampling procedure has difficulties exploring the target posterior distribution, which is an issue when the curvature of the posterior is highly variable. However, because there is only a single divergence, and there are good $\hat{R}$ values, and effective sample sizes for the bulk and tail of the distribution, there is enough confidence in the resulting posterior to move forward. 

The posterior predictive checks in @fig-MPFC-ppc below imply that the model shows adequate fit for both the imputed values of age and age$^{2}$, as well as for the response variable `MPFC`.  

```{r}
#| label: fig-MPFC-ppc
#| echo: FALSE
#| fig-cap: "Posterior predictive check for a) the MPFC response, b) imputed values of age, and c) imputed values of age-squared."
#| warning: false
knitr::include_graphics("figures\\ppc-MPFC.png")
```

The fixed effects show that there is 95% certainty that expected response in `MPFC` at 9 years and 9 months of age $\theta_{1}$ lies between $0.52$ and $1.20$. The random effects imply that there are significant differences *between* individuals within families herein, $SD(z_{1}) = [0.13 - 1.62]$, as well as significant *within* family differences, $SD(v_{1}) = [0.03 - 1.37]$. Linear development of `MPFC` at mean age $\theta_{2}$ lies between  $-0.17$ and $0.10$. In contrast to the response variables `DLPFC` and `MPFC`, the random effect estimates suggests no differences between families in linear development, $SD(z_{2}) = [0.00 - 0.31]$, but only individual-level differences within families herein, $SD(v_{2}) = [0.01 - 0.56]$. The quadratic development of `MPFC` at mean age $\theta_{3}$ is estimated to lie between $-0.01$ and $0.13$, and no significant differences between families, $SD(z_{3}) = [0.00 - 0.27]$, or within families, $SD(v_{3}) = [0.00 - 0.23]$, have been found. 

These results are visualized in @fig-MPFC-predicted. The thick black line represents predicted development in `MPFC` response over time based on the fixed effects estimates. The lighter black lines denote the uncertainty around this estimated development. The red dotted line denotes the mean age of the sample. 

```{r}
#| label: fig-MPFC-predicted
#| echo: FALSE
#| fig-cap: "Model predicted development in MPFC response."
#| warning: false
knitr::include_graphics(".\\figures\\plot_MPFC_predicted.png")
```

:::

## Associations between social emotion regulation development and social well-being (Aims 2 and 3)
The fitted multivariate regression model, predicting social well-being subscale from the individual growth components extracted from the multilevel growth curve model, is saturated (there are 0 degrees of freedom). Hence, model fit is perfect. Analyses on 47 datasets of plausible values for the individual-level growth components did not result in convergence, such that the pooled results below are based on 953 data sets. In the output, `gamma10_`, `gamma11_` and `gamma12_` represent $\hat{\gamma}_{10}$, $\hat{\gamma}_{11}$, and $\hat{\gamma}_{12}$, respectively, with the text after the underscore `_` denoting the outcome that the growth component belongs to. The term `INS` refers to the anterior insula `AI`. 

```{r phaseII, eval=TRUE, echo=FALSE}
load("C:\\Users\\5879167\\surfdrive\\Research\\Applied - Social emotion regulation\\brms\\out_phaseII.Rds")
summary(out_phaseII)
```

These results are interpreted in the main paper.