Introduction
The construct of resilience involves two key components: experiencing adversity and demonstrating positive or better-than-expected adaptation following exposure to adversity (Luthar et al., Reference Luthar, Cicchetti and Becker2000; Luthar, Reference Luthar2015; Rutter, Reference Rutter2012). Adversity refers to negative life events that are empirically found to increase risk for later adaptation problems. Positive adaptation has been assessed in resilience research using indicators of competent functioning (e.g., success in age-appropriate role performance) and low levels of mental health, physical health, or behavioral problems. Consistent with the principle of multifinality (Cicchetti & Rogosh, Reference Cicchetti and Rogosch1996), children who are exposed to major adversities such as bereavement (Dowdney, Reference Dowdney2000), divorce (Amato, Reference Amato2014), and abuse (Xiao et al., Reference Xiao, Murat Baldwin, Wong, Obsuth, Meinck and Murray2023) are at increased risk for multiple adjustment problems. How to account for the multiple adjustment problems following adversity poses an important challenge for resilience research. In this paper, we present an approach to assessing positive adaptation over time among parentally bereaved youth using indicators of competence and adjustment in five domains of functioning, which prior research has found to be associated with the adversity of the death of a parent during childhood.
An important methodological advance in resilience research over the past two decades has been the use of a person-focused approach, such as growth mixture modeling (GMM; Muthén et al., Reference Muthén, Brown, Masyn, Jo, Khoo, Yang, Wang, Kellam, Carlin and Liao2002), to group individuals with similar trajectories of adapting to adversities over time (Infurna & Grimm, Reference Infurna and Grimm2018). GMM or an alternative, latent class growth analysis (LCGA; Nagin, Reference Nagin2005), is a statistical technique used to identify unobserved, distinct trajectories of individuals’ outcomes based on similar patterns of responses over multiple time points. This approach differs from a variable-focused approach of growth curve modeling which examines the trajectory of change on a variable over time, assuming that all individuals are drawn from a single population with common population parameters (i.e., similar baseline scores and growth trajectories). GMM allows for differences of growth parameters across unobserved subpopulations (Muthén, Reference Muthén and Kaplan2004). By using GMM, we can assess trajectories of competence and adjustment in multiple domains of functioning for bereaved youth following the death of their parent.
A meta-analysis of 54 GMM studies of children and adults who had experienced a range of traumatic events (Galatzer-Levy et al., Reference Galatzer-Levy, Huang and Bonanno2018) identified four distinct trajectories of common mental health problems over time. The most common trajectory, resilience, was characterized by stable low levels of mental health problems, with a pooled prevalence rate of 65.7% across studies. Studies of children within this same review (n = 6 studies) also found that resilience was the most common trajectory across multiple single outcomes, with a 52% pooled prevalence rate. Based on these findings, the authors concluded that resilience is a common and modal trajectory following a wide range of adversities.
One caveat to this conclusion is that the review only reported on the single “primary outcomes” and did not include reports on resilient trajectories on multiple domains of problems or developmental competence. Prior research has argued that children’s resilience following adversity needs to account for variability in adjustment across different domains of functioning, including positive developmental outcomes (Luthar & Eisenberg, Reference Luthar and Eisenberg2017). Luthar et al. (Reference Luthar, Cicchetti and Becker2000) noted that heterogeneity across domains of functioning is to be expected and found that children who exhibit successful functioning in one domain often exhibit problems in other domains. Illustratively, Luthar et al. (Reference Luthar, Doernberger and Zigler1993) in a sample of inner-city adolescents found that resilience in academic functioning was not matched by resilience in emotional well-being. As children are at risk for problems in multiple domains of functioning following adversity, Luthar and colleagues (Reference Luthar, Cicchetti and Becker2000) proposed that it makes sense to assess resilience in each domain separately. We follow this approach in the current study of resilience across domains of functioning for parentally bereaved youth. Key research questions concern the prevalence of resilience for each at-risk domain as well as across domains, and the extent to which domain-specific and multidomain resilience predict later adjustment.
Prior research with adult populations has demonstrated that while resilience is often the modal trajectory for single outcomes following adversity, its prevalence varies by outcome, and resilience across multiple outcomes is rare (referred to as multidimensional resilience, Infurna & Luthar, Reference Infurna and Luthar2017a; Reference Infurna and Luthar2017b). In a series of studies of resilience for bereaved adults (Infurna & Luthar, Reference Infurna and Luthar2017a, Reference Infurna and Luthar2017b), the authors found differential rates of resilience across multiple domains of physical and mental health functioning. Specifically, the percentage of people who showed a resilient trajectory to spousal loss varied across five domains: life satisfaction (66%), negative affect (19%), positive affect (26%), general health (37%), and physical functioning (28%). Very few people (8%) had resilient trajectories across all five domains, and 20% did not show a single resilient trajectory in any of the domains (Infurna & Luthar, Reference Infurna and Luthar2017a). Similar findings were observed for parents who had experienced child loss (Infurna & Luthar, Reference Infurna and Luthar2017b). Although they report high prevalence of resilience in some domains, they concluded that resilience across multiple outcome domains is rare.
The current study assesses trajectories of resilience across multiple domains of problems and competencies on which they are at risk of problem outcomes. We use both self- and caregiver-reports of outcomes where available in a sample of parentally bereaved youth. This is only the third study to assess trajectory patterns of outcomes for parentally bereaved youth. In a sample of youth who experienced the sudden death of a parent, Melhem et al. (Reference Melhem, Porta, Shamseddeen, Walker Payne and Brent2011) identified three distinct trajectories of children’s grief over 33 months: a modal trajectory of low levels of grief that decreased over time (58.8%), moderate levels of grief that decreased over time (30.8%), and high levels of grief that remained high over time (10.4%). Using a sample of 244 parentally bereaved youth, we partially replicated these findings on trajectories of grief and extended them over six years (Sandler et al., Reference Sandler, Tein, Hoppe, Uhlman and Wolchik2025).
In the current study, we use this sample to extend research on resilience of parentally bereaved youth by assessing trajectories of 10 outcomes in five domains for which parentally bereaved youth had previously been found to be at elevated risk, including higher levels of prolonged intrusive grief (Sandler et al. Reference Sandler, Tein, Hoppe, Uhlman and Wolchik2025; Kaplow et al., Reference Kaplow, Layne, Pynoos, Cohen and Lieberman2012; Melhem et al., Reference Melhem, Porta, Shamseddeen, Walker Payne and Brent2011), internalizing problems (Berg et al., Reference Berg, Rostila and Hjern2016; Böckerman et al., Reference Böckerman, Haapanen and Jepsen2023; McKay et al., Reference McKay, Kilmartin, Meagher, Cannon, Healy and Clarke2022; Simbi et al., Reference Simbi, Zhang and Wang2020), and externalizing problems (Hamdan et al., Reference Hamdan, Mazariegos, Melhem, Porta, Payne and Brent2012; Kaplow et al., Reference Kaplow, Saunders, Angold and Costello2010; Melhem et al., Reference Melhem, Porta, Shamseddeen, Walker Payne and Brent2011; Luecken & Rubinov, Reference Luecken and Roubinov2012), as well as lower levels of academic competence (Brent et al., Reference Brent, Melhem, Masten, Porta and Payne2012; Liu et al., Reference Liu, Grotta, Hiyoshi, Berg and Rostila2022; Oosterhoff et al., Reference Oosterhoff, Kaplow and Layne2018), and peer competence (Brent et al., Reference Brent, Melhem, Masten, Porta and Payne2012). We study the prevalence of domain-specific resilience for each of these outcomes and multidomain resilience as the cumulative number of domains on which youth had resilient trajectories. Similar to research with adults (Infurna & Luthar, Reference Infurna and Luthar2017a, Reference Infurna and Luthar2017b), we expected that resilience may be the modal trajectory for individual outcomes, but very few people would be considered resilient across all domains or have no domain being considered as resilient. Finally, we test the long-term predictive value of both domain-specific and multidomain resilience for outcomes assessed 15 years after baseline.
Method
Participants
Participants included 244 children and adolescents (hereafter referred to as youth) and their caregivers from 156 families who participated in a randomized controlled trial of the Family Bereavement Program (FBP). Families had at least one youth between 8 and 16 years of age (M = 11.29, SD = 2.43, 54% boys) who experienced the death of a parent between 3 and 30 months prior to entering the trial (M = 10.81, SD = 6.35) and were fluent in English. Youth identified their race/ethnicity as follows: non-Hispanic Caucasian (67%), Hispanic (16%), African American (7%), Native American (3%), Asian or Pacific Island (1%), and Others (6%). The reported causes of death included illness (67%), accident (20%), and violence (suicide or homicide; 13%).
Procedure
All study procedures were approved by the university’s Institutional Review Board. Full details on eligibility, recruitment, and study procedure were described elsewhere (Sandler et al., Reference Sandler, Ayers, Wolchik, Tein, Kwok, Haine, Twohey-Jacobs, Suter, Lin, Padgett-Jones, Weyer, Cole, Kriege and Griffin2003; Sandler et al., Reference Sandler, Ayers, Tein, Wolchik, Millsap, Khoo, Kaplan, Ma, Luecken, Schoenfelder and Coxe2010; Sandler et al., Reference Sandler, Gunn, Mazza, Tein, Wolchik, Kim, Ayers and Porter2018). Briefly, families were recruited between 1996 and 1999 through local community agencies that had contact with bereaved youth and families (e.g., schools, churches, hospices), newspapers, and media presentations. After the baseline assessment, eligible families were randomized into either the FBP intervention condition (n = 135 youth from 90 families) or the self-study condition (n = 109 youth from 66 families). Families were interviewed and assessed on all measures used in the current study at baseline (W1) and at three months (W2; 98% retention), 14 months (W3; 90% retention), and 15-years (77% retention) post baseline. The 15-year outcomes were assessed between 2011 to 2014. Families and youth received compensation for participating in each assessment. All participants provided informed consent or assent at each assessment.
Measures
Ten measures were administered from W1 – W3 for our assessment of GMM trajectories. Measures were classified within five domains of functioning, including grief, internalizing problems, externalizing problems, academic competence, and peer competence. For all domains except grief, measures were administered to both youth and caregivers, allowing us to classify youth as resilient on each domain using the trajectories showing agreement across both reporters and by either reporter. We reported internal consistency (Cronbach’s α) of each scale at baseline. The internal consistency scores of the later waves were similar (i.e., different by the second decimal place) and thus are not reported here.
Domain 1: Grief
Intrusive grief, youth-report
Youth intrusive grief was measured by the Intrusive Grief Thoughts Scale developed for the FBP study (Sandler et al., Reference Sandler, Ma, Tein, Ayers, Wolchik, Kennedy and Millsap2010; α = .89) at W1 – W3. Youth rated 10 statements reflecting the frequency of experiencing intrusive and unwanted negative or impairing thoughts regarding the death of their parent on a 5-point Likert scale from 1 (not at all) to 5 (several times a day). Higher mean scores indicated more intrusive grief symptoms.
Domain 2: Internalizing problems
Depression, youth-report
Youth depression was measured by the Children’s Depression Inventory (CDI; Kovacs, Reference Kovacs1981; α = .87). Youth rated 27 statements from 0 (absence of symptoms) to 2 (definite symptom) in the past two weeks. Higher sum scores indicated more severe depressive symptoms.
Anxiety, youth-report
Youth anxiety was measured by the Revised Children’s Manifest Anxiety Scale (RCMAS; Reynolds & Richmond, Reference Reynolds and Richmond1978; α = .90). Youth responded with 1 (no) or 2 (yes) to 28 statements reflecting anxiety symptoms. Higher sum scores indicated more severe anxiety symptoms.
Internalizing problems, caregiver-report
Caregivers reported children’s internalizing problems using the Child Behavioral Checklist (CBCL; Achenbach, Reference Achenbach1991a) for youth under 18 years old. Age- and sex- adjusted T-scores were calculated (M = 50, SD = 10).
Domain 3: Externalizing problems
Externalizing problems, youth-report
Youth reported their own externalizing problems using the Youth Self Report (YSR; Achenbach, Reference Achenbach1991b). Age- and sex- adjusted T-scores were calculated (M = 50, SD = 10).
Externalizing problems, caregiver-report
Caregivers reported children’s externalizing problems using the Child Behavioral Checklist (CBCL; Achenbach, Reference Achenbach1991a) for youth under 18 years. Age- and sex- adjusted T-scores were calculated (M = 50, SD = 10).
Domain 4: Academic competence
Academic competence, youth- and caregiver-report
Youth self-report and caregiver-report of academic competence were measured by the subscale from the Coatsworth Competency Scale (Coatsworth & Sandler, Reference Coatsworth and Sandler1993). Youth and caregivers rated six parallel statements about the youth’s academic competence on a 4-point Likert Scale from 1 (not at all) to 4 (very much). This subscale showed good internal consistency in the current sample (α = .79 for youth-report, α = .89 for caregiver-report).
Domain 5: Peer competence
Peer competence, youth- and caregiver-report
Similar to academic competence, youth self-report and caregiver-report of peer competencies were measured by the 7-item subscale from the Coatsworth Competency Scale (Coatsworth & Sandler, Reference Coatsworth and Sandler1993). The subscale showed acceptable internal consistency in the current sample (α = .64 for youth-report, α = .78 for caregiver-report). While the youth-report showed lower-than-ideal internal consistency, it was moderately associated with caregiver report (r = .27) and youth reports of anxiety and depression (both r = .30).
Outcomes at 15 Years for analysis of long-term prediction by domain-specific and multidomain resilience
Intrusive grief thoughts, youth-report
Intrusive grief was measured by the same 10-item Intrusive Grief Thoughts Scale (α = .90) developed for the FBP study (Sandler et al., Reference Sandler, Ma, Tein, Ayers, Wolchik, Kennedy and Millsap2010).
Social detachment/Insecurity, youth-report
Social detachment/insecurity was measured by the 7-item subscale derived as a specific grief dimension that was not correlated with general grief based on a bi-factor analysis of the Inventory of Traumatic Grief (Prigerson & Jacobs, Reference Prigerson, Jacobs, Stroebe, Hansson, Stroebe and Schut2001; Sandler et al., Reference Sandler, Ma, Tein, Ayers, Wolchik, Kennedy and Millsap2010; α = .85). Youth rated each item on a 5-point Likert scale, and higher scores indicated more grief-related social detachment and insecurity.
Personal growth through grief, youth-report
Youth personal growth through grief was measured using the 12-item personal growth subscale from the Hogan Grief Reaction Checklist (HGRC; Hogan et al., Reference Hogan, Greenfield and Schmidt2001; α = .94). The scale assessed youth’s sense of having become more compassionate, tolerant, forgiving, and hopeful post-bereavement on a 5-point Likert scale from 1 (does not describe me at all) to 5 (describes me very well). Higher scores indicate more personal growth.
Internalizing problems, youth-report
Internalizing problems during the past month were assessed using the Adult Self-Report (ASR; Achenbach & Rescorla, Reference Achenbach and Rescorla2003; α = .93).
Major depressive disorder (MDD)
The binary MDD diagnosis variable was assessed by the World Health Organization World Mental Health Composite International Diagnostic Interview (CIDI; Robins et al., Reference Robins, Wing, Wittchen, Helzer, Babor, Burke, Farmer, Jablenski, Pickens, Regier, Sartorius and Towle1988). Using computerized algorithms, the CIDI defines whether the individual meets all criteria for major depressive disorder as defined by the DSM-IV and ICD-10 (0 = no MDD diagnosis, 1 = MDD diagnosis).
Externalizing problems, youth-report
Externalizing problems during the past month were assessed using the Adult Self-Report (ASR; Achenbach & Rescorla, Reference Achenbach and Rescorla2003; α = .90).
Polysubstance use, youth-report
Youth polysubstance use was measured by the sum scores of Tobacco, Alcohol, and Drug use subscales from the Adult Self-Report (ASR; Achenbach & Rescorla, Reference Achenbach and Rescorla2003). Internal consistency was not applicable for this measure.
Suicidality, youth and caregiver report
Suicidal risk was assessed using a binary variable (0 = no suicidal risk, 1 = experiencing suicidal ideation or attempts) by two items from YASR (Achenbach, Reference Achenbach1997) reported by youth and key informants (nominated by the youth). Youth reported on the items that happened in the past month (i.e., “Deliberately harms self or attempts suicide” and “Talks about killing self”) while the key informants reported on the parallel items that happened in the last six months (see Sandler et al., Reference Sandler, Tein, Wolchik and Ayers2016). Suicidal ideation or attempts were considered present if endorsed by either the youth or key informants.
Mastery, youth-report
Mastery was assessed by a modified 10-item Mastery scale (Pearlin & Schooler, Reference Pearlin and Schooler1978; α = .86). Youth rated items on a 4-point Likert scale from 1 (strongly disagree) to 4 (strongly agree). Higher scores indicate higher mastery.
General health, youth-report
Youth rated one item from SF-12 (Ware et al., Reference Ware, Kosinski and Keller1996), “how would you describe your general health” from 1 (Poor) to 5 (Excellent).
Demographic variables
Youth age at parental death (in years), sex (0 = boy, 1 = girl), cause of parental death (two dummy-coded variables for 1 = accident vs. 0 = illness; 1 = violent vs. 0 = illness), and intervention status (0 = control, 1 = FBP) were assessed at baseline.
Analytic plan
We conducted analyses in four phases: (1) examined trajectory patterns across W1 – W3 for each of the grief, mental health, and competency variables using GMM, (2) defined the resilience group for each of the 10 variables as better than expected across each of the three assessments in the trajectory using national norms or scores for bereaved children (Luthar, Reference Luthar2015; Rutter, Reference Rutter2012); identified resilience for each of the five domains (i.e., domain-specific resilience) based on agreement across reporters and reports by either reporter; and calculated multidomain resilience scores as the number of domains on which youth were considered resilient based on agreement across trajectories of both reporters or resilient trajectory assessed by either reporter, (3) examined the relationship between multidomain resilience scores and baseline demographics as well as intervention status using Poisson regressions, and (4) examined the associations between grief, mental health, mastery and general health outcomes at 15-years post-baseline and the prevalence of domain-specific and multidomain resilience using linear and logistic regressions as appropriate. Phases 1, 3, and 4 were conducted in Mplus (Muthén & Muthén, Reference Muthén and Muthén1998-2017).
Phase 1: Growth trajectory patterns
For each variable, we fit a series of GMMs to identify trajectory patterns over time following the date of parental death and estimate the posterior probability of each individual being a member of each profile. Instead of using the conventional fixed time for the assessments (i.e., W1 – W3), the time since death was the time indicator. In the current study, youth experienced parental death between 3 and 30 months prior to baseline. As a result, the timeframe for all participants ranged from 3 to 44 months. For example, youth who lost a parent three months before baseline were assessed at approximately 3, 6, and 17 months since death whereas those whose loss occurred 30 months before the baseline were assessed at 30, 33, and 44 months since death. For ease of modeling, we compressed the time unit to every 4 months. Mplus settings for conducting GMMs leverage all available data. Thus, data from youth having at least one assessment on the variable were included in the growth models. However, the exact time of death for four youth was missing, and those cases were excluded from the GMM. We accounted for family clustering by computing robust standard errors using a sandwich estimator (Muthén, 1998).
GMM is an exploratory data analysis; we started with a 1-profile model and successively increased the number of profiles by one until model fit indices leveled off or we encountered convergence issues (Ram & Grimm, Reference Ram and Grimm2009). We tested both linear and quadratic growth models. When the two growth models were comparable in fit, we selected the quadratic model only if one or more quadratic terms within the latent classes were statistically significant. To avoid getting local maximum solutions, we repeated models with multiple sets of start values and ensured that the best log-likelihood value was replicated (Muthén, Reference Muthén and Kaplan2004). There is no single fit index that is best for optimal class solution for the mixture modeling (Nylund-Gibson & Choi, Reference Nylund-Gibson and Choi2018; Tein et al., Reference Tein, Coxe and Cham2013). We thus determined the optimal number of profiles based on several fit indices and likelihood ratio tests: Bayesian information criterion (BIC; Schwarz, Reference Schwarz1978), sample-size adjusted Bayesian information criterion (saBIC; Sclove, Reference Sclove1987), and Vuong-Lo-Mendell-Rubin likelihood ratio test (LMR; Lo et al., Reference Lo, Mendell and Rubin2001). We could not perform a parameter bootstrapped likelihood ratio test (Peel & McLachlan, Reference Peel and McLachlan2000) due to the clustering effect of siblings. We also relied on entropy to gauge whether the latent profiles were highly discriminating (Nylund et al., Reference Nylund, Asparouhov and Muthén2007; Ram & Grimm, Reference Ram and Grimm2009) and substantive interpretations for model selection criteria. Good entropy is generally defined as > .80, with lower entropy being acceptable (e.g., > .70) if the other fit indices favor the solution; but there is no agreed-upon cutoff point for poor entropy (Muthén, Reference Muthén2017). Similarly, for the LMR test, a p-value ≤ .05 is preferred to indicate that the K0-class solution is significantly better than the K-1- class solution. However, in case LMR did not agree with other fit indices, we rely on information criterion (BIC, saBIC) and interpretability (Muthén, Reference Muthén2019) to identify the optimal solution. Each youth was assigned to the most likely profile based on the estimated posterior probabilities for each profile for the computation of multidomain resilience scores in the next phase.
Phase 2: Prevalence of domain-specific and multidomain resilience
After identifying the best-fitting model for each variable, we defined the resilience groups for each of the 10 outcome variables across the five domains. We defined resilience based on the widely used conceptualization of resilience as stable low problems or stable high competence (Reference Infurna and Luthar2016b; Bonanno & Diminich, Reference Bonanno and Diminich2013; Galatzer-Levy et al., Reference Galatzer-Levy, Huang and Bonanno2018; Galatzer-Levy & Bonanno, Reference Galatzer-Levy and Bonanno2016, Infurna & Luthar, Reference Infurna and Luthar2016a). While the classification of “stable low” or “stable high” trajectories is somewhat subjective, we adopt Rutter’s (Reference Rutter2012) notion of “better-than-expected” outcomes for our operationalized definition. Consistent with Melhem et al. (Reference Melhem, Porta, Shamseddeen, Walker Payne and Brent2011), “expected” outcomes were operationalized as the baseline sample mean, with additional reference to subclinical (when available) or clinical cutoffs and the substantive meaning of the rating scale when applicable. Thus, we categorized a group as resilient if the estimated trajectory remained entirely below the baseline mean for symptom measures or above the baseline mean for competence measures, with additional reference when applicable. Specifically, for the variables with clinical cutoff scores (i.e., CDI and RCMAS) or standardized norms (i.e., CBCL and YSR T-scores), the entire estimated trajectory for the resilience group was required to remain below the (marginal) clinical cutoff scores. For those without clinical cutoff scores or national norms (grief and competence scales), substantive meaning of the rating scales and a score above the mean of academic and social competence and below the mean on grief in this sample was used to satisfy the resilience definition. For each variable, youth received a score of 1 for being in the resilience group or 0 for non-resilience. This definition emphasized that resilience is not a static construct assessed at a single time point but reflects the stability of positive or better-than-expected outcomes over time.
Four of the five domains (excluding grief) included both youth- and caregiver-report measurement. As there had been long-term debate on who served as better reporter of child behaviors (e.g., Achenbach et al., Reference Achenbach, McConaughy and Howell1987), we calculated two types of domain-specific resilience: (a) resilience based on agreement across reporters (1 = agreement across youth and caregiver reports of resilience; 0 = at least one report of non-resilience), and (b) resilience by either reporter (1 = at least one report from either youth or caregiver indicated resilience; 0 = non-resilience by the other reporter). For example, youth were classified as resilient in the externalizing problems domain across reporters when both youth and caregiver reports indicated resilience, and as resilient by either reporter when resilience was indicated by at least one of the reporters. Youth were considered resilient on internalizing problems across reporters when they received a score of 1 on all three measures: youth-report anxiety, youth-report depression, and caregiver-report of internalizing problems. Youth were considered resilient on internalizing problems by either reporter if they received a score of 1 on either of the three measures. Similarly, two types of multidomain resilience scores were calculated by summing the resilience across the five domains for (a) agreement across reporters and (b) either-reporter definitions (range = 0 – 5). (Table 1–4)
Table 1. Summary of growth mixture models analyses
Note. Youth = youth self-report, CG = caregiver report.
Table 2. Frequency distribution of multidomain resilience indexes
Table 3. Poisson regression analyses of demographic variables on multidimensional resilience indexes
Note. P FDR = False Discovery Rate (FDR) adjusted p-value.
*p < .05, ** p < .01, *** p < .001.
Table 4. Regression analyses of mental health, grief, health and mastery 15-years post baseline on domain-specific and multidomain resilience
Note. All outcome measures are based on youth self-reports, except for the measure of suicidality, which combines both youth self-report and key informant report.
All regression analyses controlled for youth age at parental death, youth sex, intervention condition, and cause of parental death.
Adjusted p < .05 are bolded.
Phase 3: Baseline prediction of multidomain resilience
We tested the relations between multidomain resilience scores and demographics, including youth age at parental death, youth sex, and cause of parental death. We included the intervention condition as a control variable. The intervention condition is not expected to affect latent classifications in a randomized trial; its effect would be to lower or raise the slope within each class trajectory (Brincks et al., Reference Brincks, Perrino, Howe, Pantin, Prado, Huang, Cruden and Brown2018; Muthén et al., Reference Muthén, Brown, Masyn, Jo, Khoo, Yang, Wang, Kellam, Carlin and Liao2002). The count scores of the five-domain multidomain resilience scores for agreement across reporters and by either reporter were the dependent variables, respectively, using Poisson regression. We applied false discovery rate correction (FDR; Benjamini et al., Reference Benjamini, Drai, Elmer, Kafkafi and Golani2001) across all the p-values to adjust multiple tests separately for predicting resilience scores across reporters and by either reporter. The missing data in Phases 3 and 4 were handled with the full information maximum likelihood method (Muthén & Muthén, Reference Muthén and Muthén1998–2017).
Phase 4: Domain-specific resilience and multidomain resilience as predictors of outcomes 15-years post-baseline
We examined associations of grief, mental health, mastery and general health outcomes 15-years post-baseline with domain-specific and multidomain resilience scores. Linear regressions were used for continuous outcome variables and logistic regressions for dichotomous outcome variables. Covariates included youth age, youth sex, cause of parental death, and intervention condition. We applied FDR across all the p-values to adjust for multiple tests separately for predictions from domain-specific and multidomain resilience scores based on agreement across reporters (see Table 4, upper section) and either reporter (see Table 4, lower section).
Out of the 240 youth, 164 (68.3%) had complete data on all 10 outcomes. Missing data on single 15-year outcome measures ranged from 26.3 to 30%. Little’s MCAR test showed that 15-year outcomes data were not missing completely at random, X2(35) = 73.42, p = .0002. Independent sample t-tests and chi-square tests showed youth with (n = 76) or without missing data (n = 164) on 15-year outcomes did not significantly differ based on youth baseline age, time since parental death, youth sex, cause of death, or group assignment (ps > .05). Youth with missing data on 15-year MDD diagnosis were less likely to experience a parental death due to accident (X2(1) = 5.33, p = .02, Odds Ratio [OR] = 0.42). This indicated that the missingness could be associated with observed variables (i.e., Missing at Random). No other missingness on any single outcome variables was related to youth baseline age, time since parental death at baseline, youth sex, cause of death, or group assignment (ps > .05).
Results
Growth trajectory patterns
GMM for all 10 outcomes identified at least one resilient group. As shown in Table 1, significant variability in the proportion of youth classified in resilient groups on 10 single outcomes was found, ranging from 36.3% (n = 87) for youth-report intrusive grief to 83% (n = 198) for youth-report depression. With the exception of intrusive grief, resilience was the modal trajectory for the other nine variables as assessed by single outcomes, with a pooled prevalence of 62.3% (median = 62.1%) across 10 variables. The fit statistics, proportions, observed individual trajectories, and estimated mean trajectory with the growth parameters from the Mplus output in each class are presented in the supplemental materials.
Significant variability was found across the five domains in the proportion of youth classified as resilient, with the lowest rate being grief (36.3%). On average, 43.2% of youth (median = 44.6%) were classified as domain-specific resilient (i.e., resilient on a single domain) on a given domain based on across-reporter agreement, while 74.1% (median = 79.6%) were classified as domain-specific resilient on a given domain based on either-reporter coding. The internalizing problems domain showed the largest discrepancy, with 45% classified as resilient across reporters and 95.8% classified as resilient by either reporter. For externalizing problems, academic competence, and peer competence about one-third to half of youth were classified as resilient based on agreement across reporters, whereas more than three-quarters were classified as resilient by either reporter.
Domain 1: Grief
Intrusive grief, youth-report
Fit statistics for GMM and sample posterior probabilities in each latent class are shown in Table S1; the observed individual trajectories and mean trajectory for each class can be found in Figure S1; and the observed individual trajectories and estimated mean trajectory for each class are shown in Figure 1(a). The 3-class linear model was selected since it has lowest or near lowest BIC and saBIC, a significant LMR test, and highest entropy. A resilient (n = 87, 36%) group was identified, demonstrating lower than baseline mean (M W1 = 3.49, SD W1 = 0.86) of self-reported intrusive grief thoughts (i.e., occurrence ranged between not at all and less than once a week) across W1 – W3.
Figure 1. Growth mixture models illustrating the estimated means for each profile across all 10 variables.
Domain 2: Internalizing problems
Depression, youth-report
The 3-class quadratic growth model (Figure 1(b); Table S2; Figure S2) provided the optimal solutions based on near lowest BIC, acceptable saBIC, and highest entropy. No LMR test was significant. A resilient group (n = 198, 82.5%) was identified, demonstrating lower than baseline mean (M W1 = 9.84, SD W1 = 7.55) across W1 – W3, which is lower than the clinical cutoff score of 16 (Timbremont et al., Reference Timbremont, Braet and Dreessen2004) across W1 – W3.
Anxiety, youth-report
The 3-class quadratic growth model (Figure 1(c); Table S3; Figure S3) provided the optimal solutions based on acceptable BIC and saBIC and highest entropy. No LMR test was significant. A resilience group (n = 154, 64.2%) was identified, demonstrating lower than baseline mean (M W1 = 37.62, SD W1 = 6.88) across W1 – W3. It is noteworthy that it is also lower than the clinical cutoff of 38 (Stallard et al., Reference Stallard, Velleman, Langsford and Baldwin2001) across W1 – W3.
Internalizing problems, caregiver-report
The 4-class quadratic growth model (Figure 1(d); Table S4; Figure S4) provided the optimal solution based on near lowest BIC and saBIC, significant LMR test, and high entropy. Two resilience groups were identified (n = 149, 62.1%), demonstrating lower than baseline mean (M W1 = 54.91, SD W1 = 10.87) and below the marginal clinical level across W1 – W3.
Domain 3: Externalizing problems
Externalizing problems, youth-report
The 2-class quadratic growth model (Figure 1(e); Table S5; Figure S5) provided the optimal solution based on near lowest BIC and saBIC, highest entropy, and significant LMR test. One resilient group (n = 149, 62.1%) was identified, demonstrating lower than baseline mean (M W1 = 51.06, SD W1 = 11) across W1 – W3.
Externalizing problems, caregiver-report
The 4-class linear growth model provided the optimal solution based on lowest or near lowest BIC and saBIC, significant LMR test, and highest entropy (Figure 1(f); Table S6; Figure S6). Two resilient groups were identified (n = 137 combined across the two groups, 57.1%), demonstrating lower than baseline mean (M W1 = 53.38, SD W1 = 10.68) and below the marginal clinical level from W1 – W3.
Domain 4: Academic competence
Academic competence, youth-report
The 2-class linear growth model (Figure 1(g); Table S7; Figure S7) provided the optimal solutions based on acceptable BIC and saBIC, significant LMR test, and highest entropy. Additionally, a comparison of the 2-class linear versus 2-class quadratic models indicated that the quadratic term was not significant in any of the two classes. A resilient (n = 138, 57.5%) group was identified, demonstrating higher than baseline mean (M W1 = 3.05, SD W1 = 0.64) across W1 – W3.
Academic competence, caregiver-report
The 2-class linear growth model (Figure 1(h); Table S8; Figure S8) provided the optimal solutions based on acceptable BIC and saBIC, significant LMR test, and near highest entropy. Additionally, a comparison of the 2-class linear versus 2-class quadratic models indicated that the quadratic term was not significant in any of the two classes. A resilient group (n = 164, 68.3%) was identified, demonstrating higher than baseline mean (M W1 = 3.09, SD W1 = 0.73) across W1 – W3.
Domain 5: Peer competence
Peer competence, youth-report
The 2-class linear growth model (Figure 1(i); Table S9; Figure S9) provided the optimal solutions based on acceptable BIC and saBIC, significant LMR test, and near highest entropy. Additionally, a comparison of the 2-class linear versus 2-class quadratic models indicated that the quadratic term was not significant in any of the two classes. A resilient group (n = 182, 75.8%) was identified, demonstrating higher than baseline mean (M W1 = 3.36, SD W1 = 0.49) across W1 – W3.
Peer competence, caregiver-report
The 3-class quadratic growth model (Figure 1(j); Table S10; Figure S10) provided the optimal solutions based on acceptable BIC and saBIC, significant LMR test, and acceptable entropy. A resilient group (n = 137, 57.1%) was identified, demonstrating higher than baseline mean (M W1 = 3.30, SD W1 = 0.5) across W1 – W3.
Multidomain resilience
We assessed multidomain resilience as the count of the number of the five domains in which the youth were classified as resilient. Table 2 presents the frequency of the multidomain resilience scores based on agreement across reporters and by either reporter. Using the criteria of resilient by both reporters, on average, youth exhibited across-reporter classification as resilient on 2.16 (SD = 1.45) domains (median = 2, skewness = 0.29, kurtosis = −0.86; Figure 2(a)). Only 17 youth (7.1%) were classified as resilient across all five domains, and 31 youth (12.9%) were classified as non-resilient across all five domains. For resilience by either reporter, youth were classified as resilient by either reporter on an average of 3.70 (SD = 1.12) domains (median = 4, skewness = −0.92, kurtosis = 0.59; Figure 2(b)). Sixty-one youth (25.4%) were classified as resilient across all five domains, and two youth (0.8%) were classified as non-resilient across all five domains by either reporter.
Figure 2. Distribution of multidomain resilience scores across reporters (a, top) and by either reporter (b, bottom). Note. Although the distribution does not meet the statistical criteria for normality, we have included normal distribution curves (shown in black) in the figure for reference. This is to illustrate that the observed distribution does not strongly deviate from a normal distribution.
Baseline predictors of multidomain resilience
As shown in Table 3, none of the demographic or intervention condition variables significantly predicted multidomain resilience scores based on agreement across reporters and by either reporter after FDR adjustment.
Domain-specific and multidomain resilience as predictors of 15-year outcomes
The upper part of Table 4 presents results of predicting 15-year outcomes from the domain-specific resilient and multidomain resilient scores, where the resilient classification within each domain was based on agreement across reporters. Resilience in grief was associated with two grief-related 15-year outcomes, one indicating intrusive grief thoughts and the other indicating positive growth from grief. Resilience in internalizing problems and academic competence each predicted multiple measures of problem and positive functioning at the 15-year outcomes. Resilience in externalizing problems was predictive of lower grief-related social detachment/insecurity, externalizing problems, and polysubstance use. Resilience in peer competence did not predict any problem or positive functioning outcomes. Multidomain resilience, as reported by both reporters, predicted a broad range of grief, mental health, and mastery outcomes. Although no longer statistically significant after FDR correction (p = .055), the association between multidomain resilience and suicidality may be of interest for further investigation.
The lower part of Table 4 presents results of predicting outcomes using the criterion of being classified as resilient by either reporter. The predictions of 15-year outcomes by resilience in grief were the same since both methods of assessing the grief resilience domain only used youth report. Resilience in peer competence only predicted one of the 10 outcomes. Domain-specific resilience on internalizing problems, externalizing problems, academic competence, and multidomain resilience as classified by either reporter predicted somewhat fewer 15-year problem and positive functioning variables as compared to resilience as classified across reporters.
Discussion
The current study adds to prior literature on multidimensional resilience (Luthar et al., Reference Luthar, Cicchetti and Becker2000) and the interrelations of domains of functioning across development (Masten & Cicchetti, Reference Masten and Cicchetti2010) by describing the prevalence of problem and competence resilient trajectories across multiple domains, the distribution of cumulative resilience across domains, and the predictive relations between these resilience domains and outcomes 15-years post-baseline in a sample of parentally bereaved youth. The three main findings from this research are that bereaved youth experience differential prevalence of resilient trajectories across different domains in which they are at risk, that few bereaved youth are resilient across all such domains or have no domain on which they are resilient, and that resilient trajectories in specific domains differentially predict outcomes 15 years following baseline assessment. These findings are discussed in terms of their implications for estimating the prevalence of resilience, how understanding theoretical processes underlying resilience can inform the design of interventions to promote it across domains of functioning, and the study’s limitations and directions for future research on resilience among youth exposed to major life adversities.
An important advance in resilience assessment has been the estimation of trajectories of outcomes over repeated observations over time using growth mixture modeling (GMM; e.g., Bonanno & Diminich, Reference Bonanno and Diminich2013). Prior estimates of resilience prevalence using GMM concluded that resilience is very common based on multiple studies of resilient trajectories of functioning on a single outcome (Galatzer-Levy et al., Reference Galatzer-Levy, Huang and Bonanno2018). Others have argued that resilience is rare when assessed across multiple outcome following adversity (Infurna & Luthar, Reference Infurna and Luthar2017a; Reference Infurna and Luthar2017b). Our findings are consistent with the conclusion, that resilience is indeed uncommon, if assessed as a dichotomous variable across all five domains in which bereaved youth are at risk (i.e. multidomain resilience). Only 7.1% of youth met the criteria for resilience across all five domains based on agreement across reporters, and only 25% met the criteria for resilience based on either reporter. Yet very few youth were non-resilient across domains. Only 12.9% of youth did not meet criteria for resilience based on agreement across reporters and 0.8% of youth did not meet criteria for resilience by either reporter. Overall, cumulative multidomain resilience following the major adversity of parental bereavement followed a near-normal distribution, suggesting that most youth adapt well in some domains but not others, and that few youth have no area of resilience nor are resilient across all domains. Accordingly, multidomain resilience is better understood as a count of the number of discrete domains on which youth are resilient rather than a dichotomous phenomenon requiring resilience on all domains for which youth are at risk. To further unpack nuances and complexity of cross-domain resilience, future research could use person-centered approaches, such as latent profile analysis, to identify common patterns of resilient and non-resilient domains, which may reveal qualitatively distinct resilience profiles with different implications for intervention and long-term outcomes.
An important methodological strength of the current study is to assess resilience within domains based on agreement across multiple reporters’ resilient trajectories. Drawing from the multi-informant assessment literature (e.g., Achenbach et al., Reference Achenbach, McConaughy and Howell1987; De Los Reyes et al., Reference De Los Reyes, Augenstein, Wang, Thomas, Drabick, Burgers and Rabinowitz2015), we view the cross-reporter agreement resilience as a more stringent criterion than single reporter resilience, as it reflects consensus of observation across distinct contexts. Cross-reporters resilience may indicate generalized or cross-context resilience, while resilience reported by only one informant may reflect resilience as observed within a specific context. Prior scholarship has also emphasized that ratings from others can be especially valuable in resilience research because they capture competencies in salient developmental domains more reliably than self-reports alone (Infurna & Luthar, Reference Infurna and Luthar2018; Luthar & Zelazo, Reference Luthar, Zelazo and Luthar2003; Luthar, Reference Luthar2015). However, we present prevalence estimates for both methods to allow comparisons with prior studies that have used only a single reporter to assess trajectories of functioning over time. Not surprisingly resilience prevalence was much higher when based on either single report compared to the more stringent criterion of agreement across reporters. We also acknowledge that differences between reporters are shaped by the combination of distinct characteristics and contexts in which they observe children’s behavior (e.g., Bauer et al., Reference Bauer, Howard, Baldasaro, Curran, Hussong, Chassin and Zucker2013; Cai et al., Reference Cai, Georgeson, Basha, Lee, Xu and Gewirtz2025). Future research should aim to disentangle rater bias from meaningful contextual variation, in order to better understand the sources and implications of informant discrepancies in resilience research.
Consistent with prior research (see Luthar, Reference Luthar2015 for a review of studies on diversity of resilience across domains), we also find that resilience prevalence differs across different domains of problem and competent functioning. There was a broad range of resilient trajectories across five domains when assessed by agreement across reporters (36.3% to 48.3%; M = 2.16; median = 2) and by either reporter (36.3% to 91.3%; M = 3.70; median = 4). Similar to Galatzer-Levy et al., Reference Galatzer-Levy, Huang and Bonanno2018, these findings indicate that on some specific domains and reporters, resilience is common. Domains on which bereaved youth are resilient may provide them resources to adapt following the death. Those who were not resilient on any domain of functioning (12.9% based on agreement across reporters), are likely at increased risk for long-term difficulties of adaptation. Our conclusion is that there is no single estimate of resilience prevalence that does justice to the variety of outcomes across domains for bereaved youth. Rather, the variety of outcomes reflects a dynamic process in which adaptive processes (Masten, Reference Masten2014) are more or less successful in impacting different areas of functioning. The lack of associations between demographic and death-related factors and resilience may reflect that these predictors do not capture the underlying processes that lead to resilience.
We agree with Infurna and Luthar (Reference Infurna and Luthar2016b) that estimates of resilience prevalence have consequences for the allocation of resources to treat or prevent problem outcomes following adversity. If resilience is seen as common, then interventions to promote resilience may not be needed for the great majority of people (e.g., Bonanno et al., Reference Bonanno, Westphal and Mancini2011), and intervention resources should focus only on the few who are at risk of suffering. Our finding that most youth show resilience in some domains but struggle in others suggests that interventions should target the specific domains in which bereaved youth face difficulties, particularly those linked to long-term adverse outcomes.
Our study found that some domains of resilience, including academic competence and internalizing problems, as well as multidomain resilience predicted multiple outcomes 15-years post-baseline, whereas resilience on the grief domain had a specific effect to predict only long-term grief outcomes. We will focus our discussion on those findings that have the strongest implications for interventions and the underlying processes that may be the most productive targets for interventions.
Luthar et al. (Reference Luthar, Cicchetti and Becker2000) propose that resilience-promoting interventions should target the processes that drive positive or problematic outcomes. Research has identified empirically supported protective individual-, family-, and community-level processes that are associated with resilience of bereaved youth (Alvis et al., Reference Alvis, Zhang, Sandler and Kaplow2022; Hoppe et al., Reference Hoppe, Sandler, Williams and Winter2024; Luthar et al., Reference Luthar, Crossman, Small, Lerner and Lamb2015; Sandler et al., Reference Sandler, Wolchik and Ayers2007; Sandler et al., 2024; Tein et al., Reference Tein, Sandler, Ayers and Wolchik2006). Overall, the interactions of individual-, family-, and community-level processes shape youths’ perceptions of their world as supportive or threatening, influencing both their sense of self (e.g., internalizing problems) and their engagement with their environments (e.g., family and school). Interventions that target malleable processes within and across these domains offer proximal leverage points for strengthening resilience among bereaved youth. For example, the finding that resilience in the grief domain predicts lower intrusive grief and higher growth through grief 15-years post-baseline supports the long-term benefit of early interventions to promote adaptive grief. At the individual level, these interventions may focus on psychoeducation to normalize grief, receiving the support from others who are grieving, and promoting effective grief focused coping and emotion regulation skills (e.g. Hill et al., Reference Hill, Oosterhoff, Layne, Rooney, Yudovich, Pynoos and Kaplow2019; Sandler et al., Reference Sandler, Ayers, Tein, Wolchik, Millsap, Khoo, Kaplan, Ma, Luecken, Schoenfelder and Coxe2010). In addition, interventions could also be designed to promote family level processes that are related to intrusive grief such as caregiver facilitation of children’s grief discussion (Alvis et al., Reference Alvis, Zhang, Sandler and Kaplow2022; Hoppe et al., Reference Hoppe, Alvis, Oosterhoff and Kaplow2025).
Our finding that resilience on internalizing problems and academic competence domains predict multiple mental health, grief, and mastery outcomes indicates that multi-level interventions to promote processes that account for resilience on these domains may have long-term impacts on functioning across multiple domains of problem and positive functioning. For example, a multi-component program to teach youth effective coping and emotion regulation skills, and to teach caregivers effective parenting and grief facilitation skills had positive impact to reduce internalizing problems including major depression (Tein et al., Reference Tein, Sandler, Ayers and Wolchik2006; Sandler et al., Reference Sandler, Tein, Zhang and Wolchik2024) and improve academic outcomes (Schoenfelder et al., Reference Schoenfelder, Tein, Wolchik and Sandler2015). In a study using a nationally representative sample, Oosterhoff et al. (Reference Oosterhoff, Kaplow and Layne2018) found that bereavement due to sudden death was associated with lower academic functioning, weaker beliefs that teachers are fair, reduced school belongingness and enjoyment, and greater concentration difficulties compared to non-bereaved peers. The authors as well as others recommended school-based programs, including teacher training and academic tutoring, to support the academic success of bereaved youth.
The current study had multiple limitations. The study adopted a conceptualization of resilient trajectories as consistently lower than expected levels of problems or higher levels of competent functioning (Luthar et al., Reference Luthar, Crossman, Small, Lerner and Lamb2015; Rutter, Reference Rutter2012). This conceptualization aligns with prior research on resilience using a GMM approach (Bonanno, Reference Bonanno2004), including a meta-analysis of 54 studies with children and adults (Galatzer-Levy et al., Reference Galatzer-Levy, Huang and Bonanno2018) and was the focus of an ongoing debate in the literature regarding the prevalence of resilience (e.g., Reference Infurna and Luthar2016b, Galatzer-Levy & Bonanno, Reference Galatzer-Levy and Bonanno2016; Infurna & Luthar, Reference Infurna and Luthar2016a). However, we acknowledge that other trajectories, including those reflecting recovery of positive functioning over time, should also be considered indicators of resilience. The prevalence and the correlates of recovery trajectories in different domains is an important area of future research.
The current study investigates resilience of youth who experienced a specific major adversity, the death of a parent. Further research is needed to assess whether the findings concerning domains of resilience generalize across other acute adversities such as parental divorce, accidents, or exposure to violence, as well as to more chronic adversities, such as growing up in poverty. Moreover, and although this is not a limitation per se, it should be noted that the study investigated resilience in terms of domains of outcomes, rather than the processes or mechanisms underlying these outcomes. Future research should study these underlying mechanisms and explore how they lead to domain-specific resilience or resilience across multiple domains of functioning.
Although assessing resilience using a multidomain approach based on the GMM derived trajectories provides a useful perspective on resilience, it is important to not overinterpret the meaning of the individual trajectories. GMM is an exploratory analysis tool and using it to approximate trajectories is subject to a certain degree of uncertainty and may differ depending on the assumptions specified in the model (Infurna & Grimm, Reference Infurna and Grimm2018). To address common convergence issues with limited sample sizes, we constrained most variances of the growth factors (i.e., intercept, linear slope, quadratic slope) to zero (Hox, Reference Hox2002; Nagin, Reference Nagin1999) which might have affected the class solutions or parameter estimates. This underscores the importance of replicating the study with large samples. Furthermore, we decided the optimal solutions based on a combination of criteria since there are no definitive rules for selecting the best model (Hair & Black, Reference Hair, Black, Grimm and Yarnold2000; Ram & Grimm, Reference Ram and Grimm2009) and assigned each youth to the most likely trajectory including the resilience group for each measure. More studies of multidomain resilience in youth are needed with other samples of bereaved youth and those experiencing other adversities to replicate and extend the findings from this study. In addition, future research should move beyond the study of individual domains to explore a configural approach to identify patterns of resilient domains for which youth are at risk. Methods such as latent class analysis may illuminate these cross-domain resilience patterns and clarify which patterns of resilience are associated with what long-term problems and competencies over time.
Despite these limitations, the study highlights the value of a multidomain approach for understanding the prevalence of resilience and the associations between domain-specific and multidomain resilience with long-term outcomes among youth who have experienced one of the most stressful adversities, the death of a parent (Yamamoto et al., Reference Yamamoto, Davis, Dylak, Whittaker, Marsh and van der Westhuizen1996; McKay et al., Reference McKay, Kilmartin, Meagher, Cannon, Healy and Clarke2022). The findings significantly contribute to our understanding of the challenges youth encounter as they adapt following a major adversity and directions for the development of future research to inform the development of interventions to promote their healthy, satisfying and productive lives.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0954579425101107.
Acknowledgments
The first, second, and third authors are equal contributors to this work and are designated as co-first authors.
Funding statement
This work was supported by the National Institute of Mental Health (R21MH127288-01; 2R01 MH049155-11A) and the National Institute on Drug Abuse Grant (5T32DA039772).
Competing interests
The authors declare no competing interests.
Availability of data, code, and other materials
The analysis code and materials necessary to reproduce the analysis and replicate the findings are available from the first authors upon reasonable request.
AI
AI was not used in the development or writing of this paper.
Pre-registration
We did not qualify for pre-registration because the growth trajectories of one variable had been previously analyzed in a separate published study described in this paper.
Open access
