Resilient phenotypes among bereaved youth: a comparison of trajectory, relative, and cross-domain approaches

Abstract

Background

Bereavement is a common traumatic event associated with adverse health outcomes across the life course. Despite these risks, not all bereaved individuals experience these negative effects. Limited scientific consensus exists on how to define resilience in individuals who have experienced the death of a loved one.

Methods

Using a sample of N = 3766 youth from the Avon Longitudinal Study of Parents and Children birth cohort, we identified bereavement of a family member between ages 7 and 8.5. We derived and compared three different approaches to assess resilience among bereaved youth. Trajectory-based psychological resilience identified sub-groups with similar psychological symptom profiles between ages 6 and 16 using latent growth mixture models. Relative psychological resilience at age 16 leveraged standardized residuals from a model regressing psychological symptoms on bereavement to determine better-than-expected psychological functioning relative to bereavement status. Relative cross-domain resilience around age 16 was a sum score of the residuals approach applied to eight unique domains of health. Predictive validity of each approach was assessed using depressive symptoms at age 17.5

Results

Overall, N = 877 (23%) youth were bereaved of a family member between ages 7 and 8.5. Using latent growth mixture models, a three-class solution described 84% of bereaved youth with low and stable psychological symptoms over time, 8% with worsening symptoms, and 8% with improving yet elevated symptoms. Each relative resilience score was largely concordant with the trajectory-based approach in identifying individuals as resilient or not, though relative psychological resilience demonstrated a stronger degree of concordance than the cross-domain score. Relative psychological and cross-domain resilience exhibited moderate to low correlation, depending on the domains included (r = 0.14–0.43). For each approach, resilience significantly predicted lower depressive symptoms at age 17.5, highlighting predictive validity of these measures.

Conclusions

Psychological symptom trajectories among bereaved youth aligned with those previously identified among bereaved adults. The residual-based approach to defining resilience exhibited limited utility in the context of bereavement. When identifying risk and resilience after bereavement, researchers and clinicians must address the interplay across psychosocial and physical health domains, as bereaved youth considered resilient from a mental health perspective may benefit from intervention in other domains.

Background

The death of a loved one is among the most common and impactful traumatic events reported globally [1, 2]. In line with extensive literature documenting the widespread impacts of adverse childhood experiences [3, 4], bereavement among children and adolescents associates with adverse psychosocial outcomes across the life course, including disruptions to social and familial support structures [5,6,7] and increased risk of substance use [8], depressive symptoms and behaviors [9,10,11,12], and sleep and school problems [12]. Bereavement during childhood is also associated with adverse physical health outcomes, including cardiometabolic concerns [13] and cortisol abnormalities [14]. The high risks for these negative outcomes are not constrained to the loss of a parent or primary caregiver but extend to the loss of first- and second-degree relatives and close friends as well [8, 9]. Despite the potential risks following bereavement, not all youth who face the loss of a loved one experience these adverse effects. Evidence suggests that most individuals who experience bereavement integrate grief without lasting adverse health outcomes [15, 16]. Resilience, or the ability to adapt or maintain healthy levels of functioning in the face of trauma exposure [16,17,18], is especially salient in the context of bereavement because of the high prevalence of this exposure and the multi-dimensional adverse impacts associated. Understanding resilience and risk is critical to better informing prevention and intervention efforts among bereaved youth.

A significant barrier to progress is that approaches to identifying and quantifying resilience proliferate with little scientific consensus on the relative advantages and disadvantages of each, especially among bereaved individuals [19,20,21]. Furthermore, understanding resilience in the context of bereavement presents a unique set of considerations that do not necessarily intersect with other traumatic exposures. First, the loss of a loved one can have direct impacts on structural factors (e.g., loss of familial income, changes in education setting induced by relocation) and the social ecosystem on which children depend (e.g., reduced social network). Second, post-bereavement pathogenesis may be distinguishable from other post-trauma conditions as evidenced by extant literature demonstrating unique syndromic profiles associated with grief symptoms specifically [22, 23]. Given the evidence that bereavement during childhood is likely to affect multiple health pathways, including biological regulatory systems, multidimensional evaluations of the impact of this common exposure are critical [24]. Finally, resilience is often defined by the absence of psychopathology, assuming homogeneity in adjustment to trauma and overlooking a more comprehensive state of well-being that includes both positive psychosocial functioning and physical health [25, 26]. Approaches commonly used to address these limitations employ methods such as latent growth mixture models (LGMMs), which suggest relative consistency in psychological outcome trajectories, with most individuals sustaining low symptoms over time following a common stressor [16, 19, 27]. Other approaches assume that resilience is continuously distributed and model resilience in terms of deviation from what is expected relative to trauma exposure status [28, 29]. However, congruence across approaches is not commonly investigated, resulting in limited comparability across studies and settings [21, 30, 31]. Additionally, given the narrow focus on psychological functioning and internalizing disorders in the extant literature, individuals who may be asymptomatic in some domains of well-being, but struggling to adjust to loss in other unmeasured domains, are not typically identified [32, 33].

In response, the present analysis develops and compares three resilience constructs on the same data to determine who is resilient under multiple definitions vs. one, with each construct informed by the strengths and limitations of existing approaches. We opt to compare different constructs given the lack of research about resilience among bereaved youth and the need for additional evaluation of different resilience constructs and related issues of misclassification. To address this, we vary (1) the domains of functioning or well-being through which we quantify resilience, (2) the timing of our outcome assessments, and (3) the assumptions underlying the distribution of resilience. Specifically, we construct trajectory-based psychological resilience, relative psychological resilience, and relative cross-domain resilience. The trajectory-based approach uses LGMMs to identify classes of bereaved youth with distinct longitudinal psychological symptom profiles. Relative psychological resilience draws on standardized residuals from a linear model regressing psychological symptoms on bereavement to determine who is doing better than expected given their bereavement status; relative cross-domain resilience is based on a sum score of the residuals approach applied to eight unique domains of health. We examine concordance and correlations between constructs as well as the predictive validity of each one. The aims of these analyses are to evaluate the performance of different resilience constructs using the same data among a bereaved youth sample, and to characterize how multiple domains of functioning and different assumed distributions of resilience result in potential misclassification of post-bereavement well-being.

Methods

Study population

The sample consisted of participants enrolled in the Avon Longitudinal Study of Parents and Children (ALSPAC), a birth cohort epidemiological study of parents and their children. All pregnant women resident in Avon, UK with expected delivery dates between April 1, 1991 and December 31, 1992 were invited to participate [34,35,36]. Of the 15,447 pregnancies enrolled, there were 15,658 foetuses, 14,901 of which were alive at one year of age. Detailed health and socio-demographic data were collected via self-, maternal-, and paternal-report questionnaires as well as in-person assessment clinics. Please note that the study website contains details of all the data that is available through a fully searchable data dictionary and variable search tool: http://www.bristol.ac.uk/alspac/researchers/our-data/.

At 8 years 7 months, 8,304 mothers of youth filled out a questionnaire that inquired about, among other information, the loss of a youth’s family member since the age of 7. Of those respondents, 8,195 mothers (99%) answered questions on the youth’s loss and thus were included in this analysis. See Additional file 1 for a flow chart of participant retention across analytic samples. Selection bias could occur at this stage of bereavement ascertainment if mothers of those who were bereaved were less likely to respond to the questionnaire at 8 years 7 months than mothers of those who were not bereaved. However, given our focus on child outcomes rather than parent outcomes and the fact that participation rate by mothers remained high until mid-late adolescence [37], we did not use inverse probability weighting in our analyses.

Among 8195 youth with bereavement status ascertained, 4191 participated in an ALSPAC clinic at age 17 (51% participation rate in this clinic did not differ by bereavement status). We further restricted analyses to these youth, as many of our outcomes were measured at this clinic. Finally, we restricted to 3766 youth with at least three Strengths and Difficulties Questionnaire and three Moods and Feelings Questionnaire measures across five timepoints. The latter step was used to ensure sufficient observed outcome data for sound imputation of outcomes and covariates (more details on imputation under Statistical Analysis); given that the proportions of bereaved and non-bereaved individuals who met this last criterion were similar (90% non-bereaved vs. 88% bereaved), we did not use inverse probability weighting in our analyses.

Measures

Bereavement

At 8 years and 7 months, mothers reported on the death of a youth’s family member: “Since his 7th birthday… somebody in the family died.” Response options included: ‘Yes and he was very upset’; ‘Yes and he was quite upset’; ‘Yes and he was a bit upset’; ‘Yes but he wasn’t upset’; and ‘No did not happen.’ Relationship to person who died was not asked. To define bereavement for these analyses, responses to this question were recoded into binary loss versus no loss categories, avoiding reliance on qualitative judgments from someone other than the bereaved child to quantify the impact of this exposure [38].

Mental health

Symptoms of psychopathology were assessed using the Strengths and Difficulties Questionnaire (SDQ) [39]. The SDQ is a 25-item behavioral screening tool composed of four sub-scales capturing hyperactivity, emotional symptoms, conduct problems, and peer problems. It was administered to the mother/caregiver and completed on behalf of the study youth at the ages of 6, 9, 11, 13, and 16. Responses were recorded on a 3-point scale ranging from ‘not true’ to ‘certainly true.’ Item responses were reverse coded, and total scores were derived by summing item responses. Prorated scores–scores weighted according to the non-missing items for comparability to those based on fully observed data–were used if no more than 8 items were missing. If more than 8 items were missing, the score was set to missing. Reliability of the SDQ total difficulties score was high across administration waves at ages 6, 9, 11, and 13 (Cronbach’s alpha by time point: 0.76, 0.79, 0.79, 0.73). Reliability was not calculated for the SDQ at age 16 because item-level data were not available.

Symptoms of depression in late adolescence were assessed using the Short Moods and Feelings Questionnaire (MFQ) [40]. At 17.5 years, study adolescents were administered the MFQ through an online survey. Items responses in the ALSPAC questionnaires were on a three-point scale ranging from ‘true’ to ‘not true.’ Items were reverse coded such that higher scores indicated greater depression, and a total depression score ranging from 0 to 26 was derived by summing all items. As per the ALSPAC team’s derivation, scores were not prorated, instead requiring responses to all 13 questions. Cronbach’s alpha was 0.90, demonstrating excellent internal consistency in the current sample.

Alcohol use

Alcohol use was assessed using the Alcohol Use Disorders Identification Test (AUDIT) [41]. At 17, youth were administered the AUDIT, a 10-item questionnaire about alcohol related disorders. Questions assess frequency and quantity of alcohol consumption, dependence behaviors, and harm resulting from alcohol use. Total scores range from 0 to 40, with higher scores indicating unhealthy alcohol use. Reliability was not calculated for the AUDIT because item-level data were not available.

Functional status

School attendance/absenteeism was assessed as the frequency of school attendance. School absenteeism has been shown to correlate strongly with academic performance and is often studied as an outcome in and of itself [42, 43]. At 16, youth were asked to describe the percentage that they attended school, with response options including: none, 10%, 20%, 40%, 60%, 80%, 100%, or not registered. School attendance was transformed as a continuous measure, grouping ‘not registered’ and ‘none’ responses at 0%.

Social functioning was assessed as the number of self-reported close friends. At 17, youth were asked to report their number of close friends, with response options including: 0, 1, 2–4, 5–9, 10–13, 15–19, or 20 + . Number of friends was transformed as a continuous measure. For responses 0, 1, and 20 + , an individual was assigned values 0, 1, and 20, respectively. For response 2–4, an individual was randomly assigned a number between 2 and 4, with equal probability of each number being selected. This within-range random assignment was repeated for the remaining response options as well (5–9, 10–13, and 15–19).

Sleep duration was assessed using average amount of sleep on a school night. At 15.5 years, youth were asked about the length of time they sleep on a normal school night. Reports in hours and minutes were converted to hours with decimal values accepted.

Physical health

Cardiometabolic health was defined by both systolic blood pressure and body mass index (BMI). At 17, resting blood pressure was measured using a DINAMAP 9301 machine and was taken twice on each arm. In these analyses, we used the average of two measures for right arm systolic blood pressure, though if only one measure was taken, that was used. Body mass index (BMI) was derived from height and weight measurements as kg/m². Height and weight were measured using a Harpenden stadiometer (to the last complete mm) and the Tanita Body Fat Analyzer (Model TBF 401A; to the nearest 50 g), respectively.

Finally, body inflammation was assessed by C-reactive protein (CRP) levels. At 17, CRP was measured in mg/l from a blood assay.

Resilience constructs

We derived three different methods to define resilience: a trajectory-based approach, a single-domain relative resilience approach, and a cross-domain relative resilience approach. See Fig. 1 for a schematic of these constructs.

Fig. 1

Schematic of three resilience constructs

For the trajectory-based approach, we conducted latent growth mixture modeling (LGMM) to identify heterogeneous sub-populations with distinct psychological response patterns to bereavement. This approach has been widely utilized among bereaved adults, with results typically identifying two to five core response trajectories, including chronically elevated symptoms, acute symptoms followed by recovery, and sustained low symptoms [19, 20, 44]. In the present investigation, we modeled symptoms of psychopathology across time using the SDQ total difficulties scores from pre-loss (age 6) to post-loss (ages 9, 11, 13, and 16).

To calculate relative psychological resilience, we followed methods used previously among trauma-exposed adults [21, 28, 29] as well as in the ALSPAC sample specifically with regards to count of adverse childhood experiences [45]. We regressed SDQ total difficulties score at age 16 on bereavement and utilized the normed standardized residuals from that regression as an indicator of relative resilience. With this approach, a positive residual will reflect a healthier outcome than would be expected given bereavement exposure. A negative residual, on the other hand, would reflect a poorer outcome relative to this exposure. An advantage to this approach is that the resulting continuous score is specific to confederates in a specific sample population. We selected SDQ at age 16 (versus earlier timepoints) to harmonize with the cross-domain measures (see below) assessed around age 16. We selected SDQ over MFQ to optimize comparability with trajectory-based resilience, exploring two different ways to define resilience within one domain of health.

The third construct, what we call a relative cross-domain resilience score, defines resilience continuously and additively across eight domains of health assessed around age 16: psychological functioning, social functioning, school attendance, alcohol use, sleep duration, cardiometabolic health (BMI and blood pressure), and body inflammation. After obtaining normed standardized residuals from regressing each domain separately on bereavement, we summed these into a single cross-domain score, such that a higher score suggests better cross-domain resilience. Specifically, lower SDQ score, AUDIT score, systolic blood pressure, CRP, and BMI, as well as higher number of friends, school attendance, and sleep, reflected a positive health outcome. We considered more sleep to be positive given the higher amount of sleep deemed appropriate for children and adolescents (8–11 h) versus for adults (7–9 h) and the extremely low prevalence of sleep over 11 h reported in our analytic sample [46]. Additionally, we considered lower BMI to be positive given the relationship between BMI and cardiometabolic health that was of interest and the extremely low prevalence of underweight in our analytic sample [47].

Covariates

Sociodemographic covariates were measured during pregnancy, at birth, or before bereavement exposure. Child sex (male or female) was determined from the birth notification. Mother’s prenatal financial difficulties score (range 0–15) was assessed at 32 weeks’ gestation based on self-reported difficulties in affording the following items: food, clothing, heating, rent/mortgage, and baby items. Each item was ranked on a scale of 0 (not difficult to afford) to 3 (very difficult to afford). Mother’s highest educational attainment was self-reported at 32 weeks’ gestation. Father’s highest educational attainment was reported by the mother at this same questionnaire. Maternal smoking during pregnancy was considered present if she reported smoking one or more cigarettes in response to a query on current smoking status in a questionnaire sent between 24- and 41-weeks’ gestation. Quintiles of the Townsend deprivation score, a population measure of material deprivation, were derived based on a mother’s reported address during pregnancy. Finally, youth’s prior life events exposure between the ages of 5 and 7 was measured via maternal response to questions about 18 potentially upsetting events in the youth’s life (e.g., child admitted to hospital, child changed caretaker), assessing both event occurrence and child reaction on a 5-point scale from ‘did not happen’ to ‘yes, very upset.’ Items were coded such that higher scores indicated greater impact on child, and the total life events score ranging from 0 to 72 was derived by summing all items. This score was prorated, allowing at most half of components to be missing.

Statistical analysis

To identify latent trajectories of psychopathology symptoms, captured by the SDQ total difficulties score between ages 6 and 16, using LGMM, we specified models following Jung and Wickrama [48]. First, we fit a single class growth model without covariates to ensure convergence. Second, we progressively fit models with two to six classes, assessing model fit according to a variety of indices: Akaike information criterion (AIC); Bayesian information criterion (BIC); sample-size adjusted Bayesian information criterion (SSBIC); and entropy, which identifies how clearly the model divides the classes. Two likelihood ratio tests were also used to adjudicate model fit: Lo-Mendell-Rubin likelihood test (LMR-LRT) [49] and Bootstrap Likelihood Ratio Test (BLRT). Third, covariates were added as predictors of latent class membership in a multinomial logistic regression. Step two was then repeated using conditional models since unconditional models could result in distorted results if covariates have strong effects on intercepts and class membership. Covariates, all numerically entered, include child’s sex (female vs. male), child’s life events score between 5 and 7, maternal financial difficulties at 32 weeks of gestation, paternal and maternal educational attainment, maternal smoking status around birth (yes vs. no), and Townsend deprivation index quintiles. Final model specification was determined by fit indices, interpretability of class assignment, and model parsimony. To facilitate model stability and proper class assignment that could be generalized to different datasets from targeted populations, we fixed the variance of the slope to zero as recommended in prior literature [48, 50]. Following model development and class assignment, each class was subjectively characterized based on mean class trajectories and categorized as resilient or non-resilient accordingly. We drew on the prototypical trajectories of adjustment following a potentially traumatic event outlined by Bonanno, whereby low symptoms, or positive adaptation (allowing for temporary disruptions to these patterns following trauma), are considered resilient [16, 27].

To define relative psychological resilience, we regressed SDQ total difficulties score at age 16 (square root transformed due to high skew) on bereavement, obtaining normed standardized residuals from that regression for each bereaved individual. As mentioned, we negated all residual values for better interpretability such that a positive residual would indicate a better mental health outcome than expected relative to bereavement exposure within this cohort. A negative residual for a bereaved individual, on the other hand, would suggest poorer mental health than expected given their exposure status. With a binary exposure that induces only moderate risk for adverse psychological outcomes, the expected variance explained in our outcome is minimal. Thus, in sensitivity analyses, we added youth’s life events score between 5 and 7 as a continuous covariate in the regression to capture more variance in the stressor exposure. Specifically, the addition of this covariate results in a residual that patterns SDQ response relative to both bereavement and a count of life stressors prior to bereavement.

To define relative cross-domain resilience around age 16, we first obtained normed standardized residuals from a separate regression of each outcome on bereavement for the remaining seven outcomes. We transformed variables with high skew using square root (AUDIT, school attendance, friends) or log (BMI, CRP) transformations. As was done with the SDQ and in order to maintain both strong interpretability (i.e., higher score = stronger resilience) and consistent valence across domains, we negated residuals for AUDIT, blood pressure, BMI, and CRP domains. Finally, we summed an individual’s eight residuals (no weighting scheme used), obtaining a single summary score that captures relative resilience across eight domains. This score assumes that resilience across domains is additive and equally weighted, a convention commonly employed elsewhere in measures such as the Townsend Index [51]. In sensitivity analyses, we explored a seven-domain score that excluded SDQ to identify the extent to which that psychological domain was responsible for driving observed associations with other definitions of resilience.

To compare our continuous approaches to classifying resilience, we estimated correlation coefficients between relative psychological resilience and cross-domain score methods. To compare these two continuous resilience measures with trajectory-based classes, we explored concordance between continuous constructs and class membership using simple percentiles and summary tabulations. Finally, given that these resilience constructs have not yet been applied to a sample of bereaved children, we evaluated the predictive ability of each resilience construct by regressing depressive symptoms at 17.5 years (as measured by MFQ, square root transformed and standardized) on resilience and baseline covariates, one resilience measure at a time, within the bereaved cohort. We expect positive resilience following bereavement to significantly predict lower depressive symptoms later in adolescence, reflecting predictive validity of the construct in question. This validation technique was used previously in the ALSPAC sample with relative resilience applied to adverse childhood experiences exposure [45].

We used multiple imputation with chained equations (m = 10 imputations) to impute all outcome and covariate information among those with bereavement status ascertained. All analyses, unless otherwise specified, drew on these m = 10 imputed datasets, combining results according to Rubin’s Rules [52]. LGMMs were specified using Mplus 8.0 [53] with one randomly selected imputed dataset, as we wanted a consistent imputation model across analyses and Mplus does not support pooling of results from multiply imputed datasets, to our knowledge. All other analyses were run using R 4.1.3 [54].

Results

Among our analytic sample of N = 3766 youth, 877 (23.3%) were bereaved of a family member between ages 7 and 8.5. Overall, 44.2% of youth were male sex, including 44.7% of non-bereaved youth and 42.6% of bereaved youth (Table 1). Compared to their non-bereaved counterparts, bereaved children exhibited a higher number of stressful life events prior to bereavement (3.2 vs. 2.6), had parents who were less educated (25% vs. 29% of fathers with university degree), and had mothers who were less likely to smoke around birth (9.7% vs. 11.1%); however, these differences were small and largely non-significant.

Table 1 Cohort characteristics and missingness

Trajectory-based psychological resilience

Unconditional model

To describe longitudinal change in symptoms of psychopathology among the bereaved cohort, we first estimated a simple growth model without any covariates and compared successive solutions from one to six classes. The three-class model demonstrated superior fit to the two-class solution according to all model fit metrics except for the BLRT (Additional file 2). Additionally, while the four-class model demonstrated marginally lower AIC, BIC, and SSBIC, it exhibited lower entropy and a non-significant BLRT test compared to the three-class solution. Thus, based on fit indices as well as theoretical interpretability, we adopted the three-class solution for subsequent examination in conditional models.

Conditional model

In the next step, we introduced relevant covariates to the unconditional models to adjudicate improvement in model fit. Entropy increased marginally from 0.886 in the unconditional model to 0.889 in the conditional model (Additional file 2). Inclusion of covariates produced class membership probabilities that were nearly identical to the unconditional model. According to this three-class solution, most bereaved youth (N = 741; 84%) demonstrated low and stable SDQ total difficulties scores between ages 6 and 16 (Fig. 2). Only 8% (N = 69) demonstrated worsening SDQ scores over time, and 8% (N = 67) demonstrated improving yet elevated SDQ scores over time.

Fig. 2

LGMM class mean SDQ trajectories (trajectory-based psychological resilience). A higher SDQ score reflects worse emotional and behavioral difficulties. Grey line at 13 separates normal scores (≤13) from borderline/abnormal scores (>13) according to the parent completed SDQ. Grey shaded region reflects timeline of bereavement exposure

Finally, we examined covariates that predicted trajectory class assignment. Bereaved children whose mothers had higher prenatal financial difficulties were more likely to belong to the elevated class (vs. the low/stable class; p-value < 0.001; Additional file 3). Males and those with a higher number of stressful life events at 5 were also more likely to belong to the elevated class (p-value = 0.05 for sex; p-value = 0.02 for life events). Bereaved children whose fathers were more educated were more likely to belong to the worsening class (p-value = 0.04), whereas those whose mothers were more educated were less likely to belong to the worsening class (p-value = 0.04). Other covariates included, namely maternal smoking around birth and Townsend Index quintiles, were not significantly associated with class membership among this bereaved cohort. Overall, the trajectory-based approach to defining resilience identified three sub-groups of bereaved youth with distinct psychopathology symptoms over time and with unique sociodemographic characteristics differentiating those three classes.

Relative psychological resilience

To assess psychological symptoms relative to bereavement exposure, we regressed SDQ total difficulties score at 16 on bereavement exposure between ages 7 and 8.5 and recorded the normed standardized residuals from this regression. On average, bereavement was associated with a 0.08 SD increase in SDQ score at 16 (p-value < 0.05; Additional file 4). Among bereaved individuals, residuals ranged from -3.3 (suggestive of poor resilience relative to this bereaved cohort) to 2.3 (strong resilience). Critically, bereavement exposure explained less than 1% of the variance in SDQ score at 16, resulting in residuals that were nearly perfectly correlated with observed SDQ score (r = − 1.00, p-value < 0.001). Consequently, all interpretations of the relative psychological resilience construct must be made with caution. In sensitivity analyses, regressing SDQ score on bereavement and life events score resulted in similar conclusions but with increased variance explained (Additional file 4). These residuals were also strongly correlated with observed SDQ score (r = − 1.00, p-value < 0.001). Overall, we identified a weak association between bereavement and psychological symptoms, resulting in a resilience score that aligned closely with the psychological symptom score but that did not depend on bereavement status much at all.

Relative cross-domain resilience

Our relative cross-domain resilience score examined functioning relative to bereavement status across eight unique domains of health. This cross-domain approach revealed wide variability in multi-domain outcomes around age 16 within the bereaved cohort. First, we obtained standardized residuals from independent regressions with the seven additional outcomes, following the same process that we did for psychological functioning. For each domain, bereavement explained minimal variance in the observed outcome, resulting in residuals that were highly correlated with the observed outcome. Overall, relative resilience across these eight domains was at most moderately correlated. Our two markers of resilience in cardiometabolic health (residuals for BMI and blood pressure) were positively corelated (r = 0.32; p-value < 0.05). Body inflammation positively correlated with BMI (r = 0.27; p-value < 0.05), whereas alcohol use and social functioning were negatively correlated (r = − 0.14; p-value < 0.05; Table 2). Weaker correlations existed between other resilience domains, including psychological functioning with school attendance and social functioning, school attendance with alcohol use, and BMI with sleep duration. Standardized residuals across these eight domains were then summed to obtain a cross-domain relative resilience score. Cross-domain resilience was normally distributed, with mean 0 and standard deviation 3.4 among bereaved individuals (Additional file 5). This score covaried most strongly with the BMI domain (r = 0.55; p-value < 0.05) and least strongly with the social functioning domain (r = 0.30; p-value < 0.05). A 7-domain score that excluded psychological functioning was similarly distributed around mean 0, but with smaller standard deviation 3.1 among bereaved individuals (Additional file 6). Correlation between the 7- and 8-domain scores was high at 0.95, as expected. Overall, resilience in one domain of health did not correlate strongly with resilience in other domains of health.

Table 2 Correlation coefficients between continuous resilience constructs (relative psychological and cross-domain resilience)

Comparing resilience constructs

We found moderate to low correlation between our two cross-sectional, continuous resilience constructs. Specifically, we found a correlation of 0.43 (p-value < 0.05) between relative psychological resilience at age 16 and the 8-domain relative resilience score around age 16 that includes psychological functioning (Table 2). However, when we excluded psychological functioning from the cross-domain score, the correlation with relative psychological resilience was only 0.14 (p-value < 0.05), suggesting a small but meaningful association between relative psychological resilience and resilience in other domains. The correlation of cross-domain scores with relative psychological resilience adjusted for life events prior to bereavement followed a similar pattern.

Next, we examined the concordance of these constructs in identifying individuals as resilient by comparing the proportion of each latent class in the 50th and 90th percentiles of each continuous resilience construct. Assuming meaningful latent class assignment and continuous resilience constructs, we expect the low/stable psychological symptom class to be over-represented (e.g., > 50% in the 50th percentile) and the other two classes under-represented (e.g., < 10% in the 90th percentile) among the most resilient individuals identified by each continuous resilience construct. In other words, we define concordance here as the extent to which these constructs identify the same set of individuals as resilient or not (i.e., concordance of “resilient” group membership, rather than correlation of outcomes); this assessment will be useful in conjunction with the methodological similarities and differences of each construct. Within the low/stable symptom class, 62% (16%) were in the 50th (90th) percentile of relative psychological resilience scores, and 54% (10%) were in the 50th (90th) percentile of cross-domain scores (Table 3). Within the worsening symptom class, 3% (0%) were in the 50th (90th) percentile of relative psychological resilience scores, and 20% (2%) were in the 50th (90th) percentile of cross-domain scores. Finally, within the elevated symptom class, 11% (0%) were in the 50th (90th) percentile of relative psychological resilience scores, and 22% (3%) were in the 50th (90th) percentile of cross-domain scores. The low/stable symptom class was over-represented and the worsening and elevated symptom classes underrepresented in the top percentiles of both continuous scores, indicating meaningful concordance of trajectory-based resilience with both relative psychological and cross-domain resilience. While both scores exhibited a degree of concordance with class assignment, relative psychological resilience was more concordant with class membership than was relative cross-domain resilience.

Table 3 Concordance of continuous resilience constructs and class assignment

Finally, we assessed the predictive validity of each resilience measure by regressing depressive symptoms at 17.5 (as measured by the MFQ) on each resilience construct and other covariates. All resilience measures were significantly predictive of depressive symptoms at age 17.5 (Table 4). Class assignment and relative psychological resilience more strongly predicted depressive symptoms compared to the cross-domain score, though differences in the coefficients of determination were minor. Those in the worsening and elevated classes exhibited more depressive symptoms, on average, compared to those in the low/stable class (worsening class: β = 0.43; p-value < 0.01). Higher relative psychological resilience strongly predicted fewer depressive symptoms as well (β = − 0.19; p-value < 0.001). On the other hand, a one-SD increase in relative cross-domain resilience (8 domains) was only associated with a 0.06 lower square root MFQ score (p-value < 0.001); the association between 7-domain resilience and depressive symptoms at age was very similar.

Table 4 Predictive validity of resilience constructs with MFQ at 17.5

Discussion

To our knowledge, this is the first study to compare different measures of resilience—trajectory-based psychological resilience, relative psychological resilience, and a relative cross-domain resilience score—among bereaved youth in the same longitudinal dataset. Growth models identified three bereavement response trajectories capturing low/stable (84%), worsening (8%), and elevated (8%) psychological symptom profiles. We found at most moderate correlation (r <  = 0.32) between eight individual resilience domains, low correlation (r = 0.14) between relative psychological resilience and relative cross-domain resilience related to physical and social health (but excluding psychological health), relatively strong concordance between relative psychological resilience and corresponding trajectory classification, and relatively poor concordance between relative cross-domain resilience and trajectory classification. Each resilience measure displayed significant predictive validity with depressive symptoms at age 17.5, though the performance of the relative resilience constructs must be interpreted with caution given the high correlation between residuals and observed outcomes.

Our investigation of trajectory-based psychological resilience in a bereaved youth cohort yielded two key insights that align with extant literature on bereavement and psychopathology. First, trajectory-based results in the present study were consistent with prior work using similar methods among other bereaved adult samples [55,56,57]. Most bereaved youth in our sample maintained low and stable psychological symptom profiles over time, aligning with prior evidence demonstrating that a significant proportion of bereaved individuals will sustain low symptom trajectories over time [20]. In addition, we note that psychiatric symptoms prior to bereavement–not just at the last assessment point–differentiated psychological functioning trajectories among bereaved individuals rather significantly. This finding is supported by prior studies linking pre-trauma psychopathology to worse post-trauma sequalae and underscores the value of incorporating pre-bereavement risk into assessments of resilience following bereavement [58]. The application of the growth model approach to other domains of health and functioning (e.g., academic performance) may similarly facilitate identification of individuals at-risk for adverse outcomes after bereavement, even before the stressor occurs.

Next, the application of the residual-based approach to defining both psychological and cross-domain resilience highlighted two important limitations of existing approaches to defining resilience after trauma. First, we identify limited use of the residual-based approach to defining resilience in studies where the relationship between exposure and outcome is not meaningfully correlated. In this sample, bereavement was associated with minimal psychological risk, explaining little variance in observed psychological symptoms at age 16 and producing residuals that were nearly perfectly correlated with observed SDQ outcome values. Consequently, relative psychological resilience was just a proxy for the observed SDQ score, and interpretations of concordance and predictive validity shifted away from our resilience target and towards the SDQ scale itself. The low magnitude of association observed between bereavement and psychological symptoms may be explained by our broad definition of bereavement, which included the death of any family member; this limitation is explained in further detail below. Building on prior work that applied this residual-based method, we recommend consistent reporting of the amount of variance explained by the stressor to aid interpretability of findings [21, 29]. We further caution against future application of this resilience construct for a binary trauma exposure that exerts a small effect, as the assumption of linearity between exposure and outcome will likely not be satisfied.

Second, the low correlation observed between relative psychological resilience and the relative cross-domain resilience score illustrates that health in any one resilience domain does not reliably associate with health in a different resilience domain. It is evident that integration of psychological, physical, and social domains of functioning in future studies of resilience is critical to better capturing heterogeneity in post-bereavement well-being. We discourage classification of individuals as resilient based on psychological functioning alone because post-bereavement decrements in health may manifest across social and/or physical domains [33]. To date, recommendations for a more comprehensive assessment of resilience, including examination of a critical set of outcomes or of composite scores, remain empirically understudied [18, 33, 59]. This may be explained, in part, by the methodological complexities of modeling longitudinal, multi-level outcomes [17]. In our analysis, we weighted each resilience domain equally and assumed that domains contributed additively and independently towards the overall cross-domain resilience score. However, a burgeoning literature suggests that adolescents who face adversity yet demonstrate positive psychological adjustment are more likely to exhibit higher physiological stress over time [60, 61]. Accordingly, the influence of these domains may ultimately be modeled differently, with observed correlations and interactions, timing of outcomes assessment, and prior evidence of relative impact all taken into account. From a clinical perspective, this cross-domain heterogeneity suggests that bereaved children may benefit from assessment of social and physical health outcomes even in the absence of mental health symptom endorsement. Youth who appear resilient from a mental health perspective may require services or interventions specific to other key domains of health and functioning.

Numerous study limitations should be considered. First, we were not able to determine the relationship of the study child to the deceased family member, and we expect this relationship to considerably shape subsequent risk of adverse outcomes. For example, we expect the loss of a parent to be more impactful to a child than the loss of a distant relative–and the variance explained in psychiatric symptoms to be larger, as a result. In addition, we could not disentangle time-varying estimates of acute post-bereavement symptoms vs. longer term effects without a more unified timing for our exposure or more outcome measurements. The latter two limitations likely help to explain the modest statistical association observed between bereavement and psychiatric symptoms. Third, the residual approach required outcomes to exhibit positive versus negative functioning valence along just one axis or direction. We selected low (vs. high) BMI and long (vs. short) sleep duration to represent the “healthy” valence or direction. In reality, of course, very low BMI values or very long sleep duration are not necessarily representative of positive functioning. Fourth, prior to our imputation, SDQ scores were prorated, whereas MFQ scores were not. However, given the high internal consistency for both SDQ and MFQ, coupled with our use of imputation for missing values, we do not suspect this to be an important limitation. Finally, we imputed outcomes and covariates using a robust set of auxiliary variables but required bereavement status to be ascertained for analytic sample inclusion, without the use of inverse probability weighting. We cannot rule out selection bias, especially as it relates to bereavement status ascertainment.

In conclusion, we compared three resilience constructs on the same longitudinal data to better understand the distribution of resilience among bereaved youth. Using LGMMs, we found evidence that psychological symptom trajectories among bereaved youth are similar to those identified among bereaved adults and that these trajectories are predictive of depressive symptoms in late adolescence. The residual-based approach, on the other hand, exhibited limited utility and interpretability in the context of bereavement. Despite this limitation, we showed that individual-level resilience varied greatly across domains of health. Our explorations of cross-domain resilience call for better address of the interplay across the psychosocial and physical health domains that shape resilience after trauma. Overall, this study contributes to the growing literature documenting important heterogeneity in the impact of and response to bereavement among youth.