Online supplementary materials
``Behavioral and neural responses to social rejection: Individual differences in developmental trajectories across childhood and adolescence’’ by Mulder, Dobbelaar, and Achterberg (2024)
Introduction
This document contains the online supplementary materials to Mulder, Dobbelaar, and Achterberg (2024). 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.
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.
Research aims
This study has three main aims:
- To describe developmental trajectories of behavioral (aggressive) and brain responses of negative social feedback throughout childhood and emerging adolescence,
- to estimate the associations between development in behavioral aggression regulation and development in neural responses to social rejection, and
- 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 et al. (2022), and in Dobbelaar et al. (2022). For details on the Leiden Consortium Individual Development project, the reader is referred to Crone et al. (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, Figure 1 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.
Fourth, there is a fair amount of missing data. For the behavioral response measurements, this exists almost entirely of dropout. The noise columns in Figure 2 display the different missing data patterns across time for these measurements.
noise. DLPFC, INS, and MPFC are the neural responses in the dorsolateral prefrontal cortex, anterior insula, and medial prefrontal cortext, respectively.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 Figure 2. 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.
Statistical 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 (Bürkner 2017) in R version 4.2.2 (R Core Team 2019). The Bayesian estimation framework was chosen because it can easily accommodate censoring by integrating out censored values (Gelman et al. 2013). Moreover, missing data can be imputed as part of the model fitting process under the assumption of missing at random (MAR) (Bürkner 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 and West 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.
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)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)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)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 et al. 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 et al. 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.
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 et al. 2013).
# 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)
)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
)The term AI in the R code below actually refers to the anterior insula AI.
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
)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 et al. (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.
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.
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, Pornprasertmanit, and Schoemann 2022). Point and standard error estimates were pooled according to rules by Rubin (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.
The posterior predictive checks (ppc) of two models for noise are shown in Figure 4: 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.
The right panel of Figure 4 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. Figure 5 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.
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.
Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors. Family: MV(student, gaussian, gaussian)
Links: mu = identity; sigma = identity; nu = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: 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)
agec | mi() ~ -1 + ageW1c + wave0
agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02
Data: df_noise (Number of observations: 48976)
Draws: 2 chains, each with iter = 4000; warmup = 2000; thin = 4;
total post-warmup draws = 1000
Group-Level Effects:
~fam (Number of levels: 256)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(noise_Intercept) 0.13 0.06 0.01 0.24 1.02 159
sd(noise_con) 0.31 0.07 0.15 0.43 1.04 89
sd(noise_con:miagec) 0.09 0.04 0.01 0.16 1.01 168
sd(noise_con:miagec2) 0.02 0.01 0.00 0.05 1.00 255
Tail_ESS
sd(noise_Intercept) 502
sd(noise_con) 171
sd(noise_con:miagec) 260
sd(noise_con:miagec2) 417
~fam:id (Number of levels: 512)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(noise_Intercept) 0.31 0.05 0.21 0.39 1.04 180
sd(noise_con) 0.23 0.10 0.03 0.42 1.06 55
sd(noise_con:miagec) 0.04 0.03 0.00 0.11 1.02 208
sd(noise_con:miagec2) 0.01 0.01 0.00 0.03 1.00 447
Tail_ESS
sd(noise_Intercept) 228
sd(noise_con) 231
sd(noise_con:miagec) 334
sd(noise_con:miagec2) 506
~fam:id:wave (Number of levels: 1302)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(noise_Intercept) 0.73 0.02 0.69 0.77 1.00 423 838
sd(noise_con) 0.95 0.03 0.89 1.01 1.03 181 614
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
noise_Intercept 1.68 0.03 1.63 1.73 1.01 716 732
noise_con 1.40 0.05 1.31 1.49 1.00 692 759
agec_ageW1c 1.02 0.00 1.02 1.03 1.00 906 686
agec_wave0 2.22 0.00 2.22 2.22 1.00 915 999
agec2_ageW1c2 1.00 0.00 1.00 1.01 1.00 1014 920
agec2_ageW1c_wave0 4.50 0.01 4.49 4.51 1.00 919 828
agec2_wave02 5.05 0.01 5.04 5.06 1.00 873 733
noise_con:miagec -0.01 0.02 -0.05 0.02 1.00 869 856
noise_con:miagec2 -0.04 0.01 -0.06 -0.03 1.00 757 850
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_noise 0.25 0.00 0.25 0.26 1.00 908 977
sigma_agec 0.18 0.00 0.18 0.18 1.00 993 947
sigma_agec2 0.89 0.00 0.88 0.89 1.00 950 921
nu_noise 1.01 0.01 1.00 1.03 1.00 910 735
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).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.
Results for the DLPFC model are shown below.
Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: DLPFC | mi() ~ 1 + mi(agec) + mi(agec2) + (1 + mi(agec) + mi(agec2) | fam/id)
agec | mi() ~ -1 + ageW1c + wave0
agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02
Data: df_brain (Number of observations: 1536)
Draws: 2 chains, each with iter = 4000; warmup = 2000; thin = 4;
total post-warmup draws = 1000
Group-Level Effects:
~fam (Number of levels: 256)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(DLPFC_Intercept) 0.31 0.23 0.01 0.84 1.00
sd(DLPFC_miagec) 0.18 0.12 0.01 0.43 1.00
sd(DLPFC_miagec2) 0.06 0.05 0.00 0.16 1.00
cor(DLPFC_Intercept,DLPFC_miagec) -0.08 0.50 -0.89 0.86 1.00
cor(DLPFC_Intercept,DLPFC_miagec2) -0.23 0.50 -0.94 0.79 1.00
cor(DLPFC_miagec,DLPFC_miagec2) 0.01 0.50 -0.86 0.90 1.00
Bulk_ESS Tail_ESS
sd(DLPFC_Intercept) 519 749
sd(DLPFC_miagec) 475 746
sd(DLPFC_miagec2) 424 683
cor(DLPFC_Intercept,DLPFC_miagec) 687 775
cor(DLPFC_Intercept,DLPFC_miagec2) 550 768
cor(DLPFC_miagec,DLPFC_miagec2) 741 900
~fam:id (Number of levels: 512)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(DLPFC_Intercept) 0.38 0.26 0.02 0.97 1.00
sd(DLPFC_miagec) 0.25 0.14 0.02 0.53 1.02
sd(DLPFC_miagec2) 0.06 0.04 0.00 0.17 1.00
cor(DLPFC_Intercept,DLPFC_miagec) -0.15 0.46 -0.90 0.78 1.00
cor(DLPFC_Intercept,DLPFC_miagec2) -0.13 0.49 -0.89 0.83 1.00
cor(DLPFC_miagec,DLPFC_miagec2) 0.02 0.48 -0.83 0.86 1.01
Bulk_ESS Tail_ESS
sd(DLPFC_Intercept) 395 605
sd(DLPFC_miagec) 196 540
sd(DLPFC_miagec2) 404 562
cor(DLPFC_Intercept,DLPFC_miagec) 477 799
cor(DLPFC_Intercept,DLPFC_miagec2) 717 749
cor(DLPFC_miagec,DLPFC_miagec2) 507 841
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
DLPFC_Intercept -0.38 0.15 -0.65 -0.09 1.00 770 908
agec_ageW1c 1.02 0.00 1.02 1.03 1.00 978 960
agec_wave0 2.22 0.01 2.21 2.23 1.01 1010 991
agec2_ageW1c2 1.00 0.01 0.99 1.02 1.00 1030 960
agec2_ageW1c_wave0 4.49 0.03 4.43 4.55 1.00 1053 922
agec2_wave02 5.04 0.03 4.98 5.10 1.00 1079 1028
DLPFC_miagec -0.01 0.06 -0.14 0.11 1.00 889 876
DLPFC_miagec2 0.05 0.03 -0.02 0.11 1.00 894 1025
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_DLPFC 3.23 0.09 3.04 3.41 1.00 482 723
sigma_agec 0.18 0.00 0.17 0.18 1.01 955 939
sigma_agec2 0.87 0.02 0.84 0.91 1.00 932 822
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The \(\hat{R}\)-values are below the recommended 1.05, indicating adequate mixing of the posterior sampling chains. The posterior predictive checks in Figure 7 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.
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 Figure 8. 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.
Results for the AI model are shown below. Note that the term INS here refers to the anterior insula AI.
Warning: There were 1 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: INS | mi() ~ 1 + mi(agec) + mi(agec2) + (1 + mi(agec) + mi(agec2) | fam/id)
agec | mi() ~ -1 + ageW1c + wave0
agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02
Data: df_brain (Number of observations: 1536)
Draws: 2 chains, each with iter = 4000; warmup = 2000; thin = 4;
total post-warmup draws = 1000
Group-Level Effects:
~fam (Number of levels: 256)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(INS_Intercept) 0.29 0.21 0.01 0.81 1.00
sd(INS_miagec) 0.14 0.10 0.01 0.35 1.01
sd(INS_miagec2) 0.05 0.04 0.00 0.14 1.00
cor(INS_Intercept,INS_miagec) -0.12 0.51 -0.91 0.85 1.00
cor(INS_Intercept,INS_miagec2) -0.18 0.52 -0.93 0.83 1.00
cor(INS_miagec,INS_miagec2) -0.11 0.51 -0.92 0.89 1.00
Bulk_ESS Tail_ESS
sd(INS_Intercept) 613 610
sd(INS_miagec) 417 848
sd(INS_miagec2) 513 695
cor(INS_Intercept,INS_miagec) 820 910
cor(INS_Intercept,INS_miagec2) 794 861
cor(INS_miagec,INS_miagec2) 959 943
~fam:id (Number of levels: 512)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(INS_Intercept) 0.46 0.28 0.02 1.01 1.01
sd(INS_miagec) 0.37 0.17 0.03 0.66 1.00
sd(INS_miagec2) 0.06 0.05 0.00 0.17 1.00
cor(INS_Intercept,INS_miagec) -0.44 0.46 -0.96 0.76 1.00
cor(INS_Intercept,INS_miagec2) -0.07 0.50 -0.90 0.88 1.00
cor(INS_miagec,INS_miagec2) -0.36 0.48 -0.97 0.75 1.00
Bulk_ESS Tail_ESS
sd(INS_Intercept) 273 544
sd(INS_miagec) 251 534
sd(INS_miagec2) 344 651
cor(INS_Intercept,INS_miagec) 184 304
cor(INS_Intercept,INS_miagec2) 549 658
cor(INS_miagec,INS_miagec2) 544 798
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
INS_Intercept 0.73 0.15 0.42 1.03 1.00 768 876
agec_ageW1c 1.02 0.00 1.02 1.03 1.00 1065 979
agec_wave0 2.22 0.01 2.21 2.23 1.00 880 905
agec2_ageW1c2 1.00 0.01 0.99 1.02 1.00 940 987
agec2_ageW1c_wave0 4.49 0.03 4.43 4.55 1.00 909 907
agec2_wave02 5.04 0.03 4.98 5.10 1.00 928 838
INS_miagec -0.16 0.07 -0.29 -0.02 1.00 1156 1030
INS_miagec2 0.06 0.03 -0.01 0.13 1.00 987 941
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_INS 3.24 0.11 3.02 3.44 1.00 300 667
sigma_agec 0.18 0.00 0.17 0.18 1.00 975 860
sigma_agec2 0.87 0.02 0.84 0.91 1.00 971 875
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).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 Figure 9 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.
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 Figure 10. 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.
Results for the MPFC model are shown below:
Warning: There were 1 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: MPFC | mi() ~ 1 + mi(agec) + mi(agec2) + (1 + mi(agec) + mi(agec2) | fam/id)
agec | mi() ~ -1 + ageW1c + wave0
agec2 | mi() ~ -1 + ageW1c2 + ageW1c_wave0 + wave02
Data: df_brain (Number of observations: 1536)
Draws: 2 chains, each with iter = 4000; warmup = 2000; thin = 4;
total post-warmup draws = 1000
Group-Level Effects:
~fam (Number of levels: 256)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(MPFC_Intercept) 0.93 0.39 0.13 1.62 1.01
sd(MPFC_miagec) 0.11 0.08 0.00 0.31 1.00
sd(MPFC_miagec2) 0.10 0.07 0.00 0.27 1.02
cor(MPFC_Intercept,MPFC_miagec) 0.06 0.49 -0.85 0.89 1.00
cor(MPFC_Intercept,MPFC_miagec2) -0.49 0.46 -0.97 0.69 1.01
cor(MPFC_miagec,MPFC_miagec2) -0.01 0.49 -0.87 0.89 1.00
Bulk_ESS Tail_ESS
sd(MPFC_Intercept) 258 622
sd(MPFC_miagec) 661 946
sd(MPFC_miagec2) 189 195
cor(MPFC_Intercept,MPFC_miagec) 750 635
cor(MPFC_Intercept,MPFC_miagec2) 423 582
cor(MPFC_miagec,MPFC_miagec2) 598 808
~fam:id (Number of levels: 512)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(MPFC_Intercept) 0.55 0.38 0.03 1.37 1.01
sd(MPFC_miagec) 0.24 0.14 0.01 0.56 1.01
sd(MPFC_miagec2) 0.09 0.06 0.00 0.23 1.00
cor(MPFC_Intercept,MPFC_miagec) 0.22 0.47 -0.74 0.91 1.01
cor(MPFC_Intercept,MPFC_miagec2) -0.20 0.51 -0.94 0.84 1.00
cor(MPFC_miagec,MPFC_miagec2) 0.22 0.47 -0.77 0.91 1.01
Bulk_ESS Tail_ESS
sd(MPFC_Intercept) 217 345
sd(MPFC_miagec) 251 360
sd(MPFC_miagec2) 234 305
cor(MPFC_Intercept,MPFC_miagec) 274 599
cor(MPFC_Intercept,MPFC_miagec2) 326 697
cor(MPFC_miagec,MPFC_miagec2) 396 758
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
MPFC_Intercept 0.87 0.17 0.52 1.20 1.00 882 964
agec_ageW1c 1.02 0.00 1.02 1.03 1.00 909 834
agec_wave0 2.22 0.01 2.21 2.23 1.00 972 887
agec2_ageW1c2 1.00 0.01 0.99 1.02 1.01 831 594
agec2_ageW1c_wave0 4.49 0.03 4.43 4.55 1.00 825 986
agec2_wave02 5.04 0.03 4.98 5.10 1.00 773 682
MPFC_miagec -0.04 0.07 -0.17 0.10 1.00 842 912
MPFC_miagec2 0.06 0.04 -0.01 0.13 1.00 796 593
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_MPFC 3.42 0.11 3.19 3.64 1.00 328 694
sigma_agec 0.18 0.00 0.17 0.18 1.00 1024 1029
sigma_agec2 0.87 0.02 0.84 0.91 1.00 865 887
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).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 Figure 11 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.
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 Figure 12. 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.
References
Achterberg, Michelle, Jeroen Mulder, Simone Dobbelaar, Stephan Heunis, and Eveline Crone. 2022. “Individual Differences in Developmental Trajectories of Social Emotion Regulation from Childhood to Emerging Adolescence.” Preregistration. https://doi.org/10.17605/OSF.IO/HDRZC.
Bürkner, Paul-Christian. 2017. “Brms : An r Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1). https://doi.org/10.18637/jss.v080.i01.
Crone, Eveline A., Michelle Achterberg, Simone Dobbelaar, Saskia Euser, Bianca Van Den Bulk, Mara Van Der Meulen, Lina Van Drunen, Lara M. Wierenga, Marian J. Bakermans-Kranenburg, and Marinus H. Van IJzendoorn. 2020. “Neural and Behavioral Signatures of Social Evaluation and Adaptation in Childhood and Adolescence: The Leiden Consortium on Individual Development (L-CID).” Developmental Cognitive Neuroscience 45 (October): 100805. https://doi.org/10.1016/j.dcn.2020.100805.
Dobbelaar, Simone, Michelle Achterberg, Lina Van Drunen, Anna C. K Van Duijvenvoorde, Marinus H Van IJzendoorn, and Eveline A Crone. 2022. “Development of Social Feedback Processing and Responses in Childhood: An fMRI Test-Replication Design in Two Age Cohorts.” Social Cognitive and Affective Neuroscience, June, nsac039. https://doi.org/10.1093/scan/nsac039.
Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. Chapman; Hall/CRC. https://doi.org/10.1201/b16018.
Jorgensen, Terrence D., Sunthud Pornprasertmanit, and Alexander M. Schoemann. 2022. “semTools: Useful Tools for Structural Equation Modeling.” https://CRAN.R-project.org/package=semTools.
Mehta, Paras D., and Stephen G. West. 2000. “Putting the Individual Back into Individual Growth Curves.” Psychological Methods 5 (1): 23–43. https://doi.org/10.1037/1082-989X.5.1.23.
R Core Team. 2019. “R: A Language and Environment for Statistical Computing.” Vienna, Austria: R Foundation for Statistical Computing.
Rosseel, Yves. 2012. “Lavaan : An r Package for Structural Equation Modeling.” Journal of Statistical Software 48 (2). https://doi.org/10.18637/jss.v048.i02.
Rubin, Donald B., ed. 1987. Multiple Imputation for Nonresponse in Surveys. Wiley Series in Probability and Statistics. Hoboken, NJ, USA: John Wiley & Sons, Inc. https://doi.org/10.1002/9780470316696.
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.