In Silico Psycho-oncology: Understanding Resilience Pathways in Breast Cancer - Determinants of Longitudinal Depression and Quality of Life Trajectories

1. Introduction

Patient resilience is a critical determinant of cancer treatment outcomes, influencing both psychological adaptation and overall quality of life [1]. Higher levels of resilience have been associated with younger age, female sex, higher socioeconomic status and being married, as well as with internal psychological resources such as optimism, self-efficacy and adaptive emotion regulation strategies and external resources such as social support [2,3,4]. Positive cognitive strategies, including acceptance and positive thinking, are linked to enhanced psychological well-being, whereas maladaptive strategies such as rumination and catastrophizing are consistently associated with adverse emotional outcomes [5,6,7,8,9,10,11]. Among women with breast cancer, optimism and self-efficacy predict better psychological well-being, improved quality of life and lower distress over time [12,13], while perceived social and family support facilitates adaptive coping and resilience [14,15,16,17].

Within this context, the EU-funded BOUNCE project (https://www.bounce-project.eu/; accessed 20 December 2025) conducted a prospective multicenter clinical study across Finland, Italy, Portugal and Israel to investigate predictors of resilience trajectories among women with breast cancer. The study collected a thorough range of sociodemographic, clinical, psychological and functional factors with the aim of improving the ability to predict resilience in response to breast cancer and, ultimately, to inform personalized interventions that promote effective psychological recovery.

The aim of the current study is twofold: (1) to identify distinct longitudinal trajectories of depression and quality of life (QoL) over the 18 months following diagnosis and/or surgery and (2) to identify factors associated with these trajectories, with particular focus on patterns reflecting resilience, recovery, deterioration and persistent impairment. To address the first aim, we apply latent class growth analysis (LCGA) and growth mixture modelling (GMM), methods widely used in psycho-oncology to capture heterogeneity in longitudinal symptom and QoL patterns. To address the second aim, we use logistic regression models to examine associations between candidate predictors and trajectory membership, allowing results to remain clinically interpretable. In addition to univariate analyses, we perform multivariable feature selection using elastic net regularization. Feature selection is used to identify parsimonious sets of predictors that jointly differentiate trajectory groups, while accounting for correlations between variables and limiting overfitting. This approach moves beyond isolated univariate associations and supports clearer interpretation of factors relevant for personalized risk assessment and intervention.

To capture clinically meaningful heterogeneity in longitudinal mental health and quality-of-life (QoL) outcomes, we focus on selected comparisons between trajectory groups that address complementary research questions. Specifically, we examine:

• Who remains resilient over time, compared with individuals who experience persistently elevated depressive symptoms or persistently lower QoL, thereby capturing differences in long-term outcome levels.
• Who is at risk for persistently poor depression or QoL outcomes, providing insight into profiles associated with sustained vulnerability.
• Who deteriorates despite early resilience, a contrast that is less confounded by baseline outcome levels and enables the identification of early warning markers relevant to preventive strategies, clinical monitoring and early intervention.
• Who recovers among individuals with comparable baseline levels, a comparison that is likewise less influenced by baseline outcome levels and highlights factors associated with improvement rather than symptom burden, with potential implications for therapeutic intervention.

2. Materials and Methods

2.1. Participants

The BOUNCE project (“Predicting Effective Adaptation to Breast Cancer to Help Women to BOUNCE Back”) is a multicenter prospective study designed to identify factors associated with psychological adaptation following breast cancer diagnosis and treatment (https://www.bounce-project.eu/; accessed 20 December 2025). Women diagnosed with breast cancer were recruited from four European countries: Finland (Helsinki University Hospital), Israel (Shaare Zedek Medical Center and Rabin Medical Center, coordinated by the Hebrew University of Jerusalem), Italy (European Institute of Oncology) and Portugal (Champalimaud Clinical Centre). Participants were enrolled approximately three to four weeks from diagnosis.

Eligible participants were female patients aged 40–70 years at diagnosis with histologically confirmed, operable invasive breast cancer (tumor stage I–III), who were receiving surgery as part of local treatment and any form of systemic therapy for breast cancer and who provided written informed consent signed by both the patient and the treating physician. Exclusion criteria included refusal to consent; presence of metastatic disease or a history of another malignancy within the previous five years; a history of severe psychiatric, neurological, (before 40 years of age) or other serious medical conditions; major surgery within the four weeks preceding enrollment; ongoing treatment for another invasive cancer or major illness; and pregnancy or breastfeeding at the time of recruitment. Ethical approval was obtained from the Ethics Committee of the European Institute of Oncology (Approval No. R868/18-IEO916) and from the corresponding ethics committees of all participating centers, while all participants provided written informed consent prior to inclusion.

Data were collected at seven assessment points at three-month intervals, from baseline (M0) to 18 months of follow-up (M18). Psychological symptoms and subjective health status were assessed at all time points, whereas sociodemographic, lifestyle and medical or disease-related variables considered in this analysis were collected at baseline [study protocol at [18]. After applying eligibility criteria and data completeness requirements, the final analytic sample consisted of 538 patients with baseline psychological assessments, at least one long-term follow-up assessment (month 12, 15, or 18) and no more than three missing assessments across the 18-month follow-up period.

2.2. Measures

2.2.1. Outcome Variables

Two outcome variables are considered reflecting psychological distress and health-related quality of life.

Psychological distress was measured using the Depression subscale of the Hospital Anxiety and Depression Scale [19] (HADS). It derives from seven items of the HADS questionnaire, with higher scores (range 0-3 – mean scale) indicating greater severity of depressive symptoms. HADS Depression score was interpreted using established thresholds, with scores of 0–1 indicating low depressive symptoms (corresponding to 0–7 on the 0–21 scale), 1.14–1.42 moderate symptoms (8-10 on the 0-21 scale) and ≥1.57 high symptomatology (≥11 on the 0-21 scale) [20]. Longitudinal changes were interpreted based on published minimally important difference criteria, whereby changes of approximately 0.2–0.43 points were considered clinically meaningful and changes of ≥0.43 points indicative of moderate to large clinical change [21,22,23].

Overall health-related quality of life was assessed using the Global Health Status/Quality of Life [24] (GHS/QoL) scale of the European Organisation for Research and Treatment of Cancer Quality of Life Questionnaire C30 (EORTC QLQ-C30). This scale reflects patients’ self-reported overall health and quality of life and is derived from two items of the EORTC QLQ-C30, with higher scores (range 0–100) indicating better well-being. As no official clinical cut-offs are defined, interpretation of absolute GHS/QoL levels was guided by published empirical thresholds, with scores ≥70 considered indicative of good to excellent quality of life, scores between approximately 50 and 69 reflecting moderate quality of life and scores <50 indicating poor quality of life [25,26]. Longitudinal changes in GHS/QoL were interpreted using established minimally important difference criteria, with changes of 5–10 points considered small but clinically meaningful and changes ≥10 points considered moderate to large [27].

2.2.2. Sociodemographic, Lifestyle and Clinical Data

The analysis incorporated socioeconomic and lifestyle factors, health-related background variables and tumor and treatment characteristics collected at baseline (Table 1).

Socioeconomic and lifestyle data collected at baseline included age, body mass index (BMI), ethnicity (Portugal; Italy; Finland; Israel), educational attainment (non-University; University), marital status (single/engaged; married/common-law; divorced/widowed), employment status (employed full-/part-time or self-employed; unemployed/housewife; retired) and monthly household income (low; middle; high). Low monthly income was defined as ≤1,000 EUR in moderate-income countries (Portugal, Italy) and ≤1,500 EUR in higher-income countries (Finland, Israel), whereas high income was defined as >3,000 EUR and >3,500 EUR, respectively.

Lifestyle factors included physical activity level (none; low/moderate; heavy), dietary pattern (no specific diet; Mediterranean/vegetarian-type; special diet), alcohol consumption (no consumption; moderate consumption; heavy consumption) and smoking status (current smoker; never smoker; former smoker). Heavy alcohol consumption was defined as intake of more than three drinks on any day or more than seven drinks per week. Heavy physical activity was defined as ≥200 minutes per week of moderate aerobic activity or ≥100 minutes per week of vigorous aerobic activity; or ≥5 weekly strength-training sessions; or combined aerobic and strength activity meeting ≥100–180 minutes per week of moderate (or ≥50–90 minutes per week of vigorous) aerobic exercise with ≥1–4 weekly strength sessions.

Health history variables included presence of chronic diseases (no/yes), metabolic diseases (no/yes), mental illness (no/yes), exposure to negative life events (none; one event; two or more events) and family history of breast cancer (no/yes). These variables were considered background health factors potentially, treatment tolerance and psychological outcomes.

Clinical and cancer-related data included menopausal status prior to cancer diagnosis (pre-/peri-menopausal; postmenopausal), use of hormone replacement therapy before diagnosis (no/yes), tumor stage (I, II, III), tumor grade (I, II, III) and histological subtype (ductal; lobular; other). Tumor biomarker characteristics comprised estrogen receptor (ER) status (negative/positive), progesterone receptor (PR) status (negative/positive), human epidermal growth factor receptor 2 (HER2) status (negative/positive) and Ki67 proliferation index (<20%; ≥20%). Molecular subtypes were defined as Luminal A-like (ER+, PR+, HER2−, low Ki67 <20%); Luminal B-like (HER2−: ER+, PR+/− with high Ki67 ≥20% or ER+, PR− with any Ki67; or HER2+: ER+, PR+/− with any Ki67); HER2-positive non-luminal (ER−, PR−, HER2+); and triple-negative (ER−, PR−, HER2−). Treatment-related variables included type of surgery (lumpectomy; mastectomy), receipt of radiotherapy (no/yes), systemic treatment modality (chemotherapy only [± anti-HER2 therapy], endocrine therapy only, or combined chemotherapy and endocrine therapy [± anti-HER2 therapy]), use of anti-HER2 therapy (no/yes) and administration of neoadjuvant chemotherapy (no/yes).

2.2.3. Psychological Scales

Psychological variables were assessed using validated self-report questionnaires administered at predefined time points throughout follow-up. Relatively stable personality and psychosocial characteristics were assessed at baseline, including optimism using the Life Orientation Test–Revised [28] (LOT-R; 10 items), sense of coherence using the Sense of Coherence Scale [29] (SOC; 13 items assessing meaningfulness, comprehensibility and manageability), trait resilience using the Connor–Davidson Resilience Scale [30] (CD-RISC; 10 items), dispositional mindfulness using the Mindful Attention Awareness Scale [31] (MAAS; 15 items) and general coping capacity using the Perceived Ability to Cope with Trauma Scale [32] (PACT; 20 items). Cancer coping self-efficacy was measured using the brief Cancer Behavior Inventory [33] (CBI-B; 12 items) and Fear of cancer recurrence was also assessed using the Fear of Cancer Recurrence Scale–Short Form [34] (FCR-SF; 9 items), both at baseline. Coping responses to cancer were assessed at three months post-diagnosis using the Mini–Mental Adjustment to Cancer Scale [35] (Mini-MAC; 29 items assessing helplessness–hopelessness, anxious preoccupation, cognitive avoidance and fighting spirit; the fatalism dimension was excluded from analyses). Post-traumatic growth was measured at three months using the Post-Traumatic Growth Inventory–Short Form [36] (PTGI-SF; 10 items).

Anxiety symptoms were assessed every three months using the the Hospital Anxiety and Depression Scale [19] (HADS; 7 items). Emotional functioning was assessed longitudinally every three months using the Positive and Negative Affect Schedule [37] (PANAS; 20 items), while health-related quality of life, including patients’ functioning status and symptom burden (e.g., arm symptoms and treatment-related side effects), was assessed every three months using the European Organisation for Research and Treatment of Cancer Quality of Life Questionnaire-Core 30 [24] (EORTC QLQ-C30) and its breast cancer–specific module [24] (EORTC QLQ-BR23).

Several single-item measures were additionally collected to capture specific psychosocial and behavioral aspects. These included a single item assessing what patients reported doing to cope with cancer, a general self-efficacy item and a perceived social support item, all assessed every three months. Adherence to medical advice was assessed using item 5 from the Medical Outcomes Study [38] (MOS) adherence questionnaire.

2.3. Statistical Analysis

2.3.1. Missing Data

Missing data were addressed using multiple imputation by chained equations (mice R package [39]), generating 30 or 50 imputed datasets, depending on the proportion of missingness. All candidate predictors, outcomes and the clinical site variable were included in the imputation models to preserve associations among variables. The random forest–based method was selected.

2.3.2. Derivation of Mental Health and GHS/QoL Trajectories

Latent class growth analysis (LCGA) and growth mixture modelling (GMM) were conducted using the lcmm package in R [40] to identify discrete longitudinal trajectories among breast cancer survivors based on EORTC QLQ-C30 GHS/QoL and HADS Depression scores over an 18-month follow-up period. Separate analyses were performed for each outcome. LCGA is a restricted form of GMM in which within-class variances of the random intercept and slope are fixed to zero, allowing heterogeneity only between classes. As random effects are not permitted within classes, LCGA typically requires a larger number of latent classes to adequately capture variability in the data compared with GMM.

To account for non-normal outcome distributions, a six-knot I-spline link function was applied. The choice of link function was informed by preliminary analyses of the null (single-class) model (Appendix A).

Quadratic models of the change across time were considered. For GMM, models with one to six latent classes were estimated, while for LCGA models with one to eight latent classes were considered. Each model was estimated multiple times using a grid of 250 sets of initial values to reduce the risk of convergence to local maxima. For each specified number of latent classes, the model with the highest log-likelihood was retained.

Model selection and determination of the optimal number of latent classes were based on a combination of statistical fit indices, class separation, trajectory shapes, minimum class size, interpretability and clinical relevance [41]. Fit was assessed using the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC) and Integrated Complete Likelihood (ICL), with lower values indicating better fit [41]. Both AIC and BIC balance model fit and complexity, but BIC imposes a stronger, sample-size–dependent penalty and therefore typically selects more parsimonious solutions with fewer latent classes. The ICL further penalizes solutions with poor class separation, jointly accounting for model fit and classification quality. Class separation was further evaluated using entropy and average posterior class membership probabilities, with values greater than 0.60 and 0.70, respectively, indicating adequate separation [42,43]. A minimum class size of at least 5% of the total sample was required.

The lcmm framework accommodates missing outcome data, therefore, no imputation of outcome variables was performed.

2.3.3. Determinants of Mental Health and GHS/QoL Trajectories

Once the trajectory groups were identified, binomial logistic regression was conducted to determine which factors at baseline, three-month and six-month follow-up are associated with trajectory group membership.

Initially, associations between individual predictors and the outcome were examined using logistic regression models adjusted for clinical site across imputations. Clinical site was included as a fixed-effect adjustment variable to account for between-site heterogeneity and results were pooled using Rubin’s rules to obtain odds ratios and 95% confidence intervals that reflect uncertainty due to missing data.

Variable selection was conducted using penalized logistic regression with an elastic net penalty, implemented in R using the glmnet package [44] (glmnet and cv.glmnet functions). The elastic net approach was chosen to balance variable selection and coefficient shrinkage while accommodating correlated predictors, which are common in clinical, quality-of-life, and psychosocial data. This strategy supported the identification of robust predictors while maintaining model interpretability and reducing the risk of overfitting. The elastic net mixing parameter was set to α = 0.5, providing an equal balance between L1 (LASSO) and L2 (ridge) penalties to achieve stable variable selection in the presence of correlated predictors [45]. Candidate predictors were encoded using model matrices (model.matrix), such that categorical variables were represented by indicator variables. The regularization parameter (λ) was selected via 10-fold cross-validation with stratified folds to preserve outcome class proportions within each fold and mitigate the impact of class imbalance. To promote stable and parsimonious models and reduce the risk of overfitting, the one-standard-error rule (λ₁se) was used for λ selection [46]. Cross-validation was performed using deviance as the optimization criterion, corresponding to the negative log-likelihood of the logistic model, thereby favoring predictors that improve probabilistic calibration rather than threshold-dependent classification accuracy. The selection procedure was repeated 3 times across all multiple imputed datasets generated using the mice package (see 3.3.1) and predictor importance was quantified by selection frequency across imputations and 3 repeats. To evaluate the robustness of our predictors, we conducted the stability analysis under two conditions. Clinical sites were treated as standard predictors subject to the L1/L2 penalty and clinical sites were forced into the model by setting their penalty factor to zero, ensuring their inclusion in every iteration regardless of the regularization threshold. This allowed us to distinguish site-independent predictors from variables reflecting differences in healthcare systems, socioeconomic contexts, and cultural factors. Predictors with high selection frequency (>60%) under both penalized and unpenalized site conditions were considered robust. Descriptive penalized odds ratios were reported to indicate the direction and relative magnitude of associations at the selected regularization parameter under the initial penalized site condition.

Model performance was evaluated using cross-validated predictions from elastic net logistic regression models, with metrics computed on held-out data and averaged across folds, thus reflecting out-of-sample performance. Performance was assessed using ROC-AUC, log-loss, and the Brier score, capturing complementary aspects of predictive quality. ROC-AUC quantifies discrimination, whereas log-loss and the Brier score assess probabilistic accuracy and calibration, with lower values indicating better performance. Although no strict cut-offs exist, ROC-AUC values of 0.60-0.70 indicate poor-to-fair discrimination, 0.70-0.80 acceptable, 0.80-0.90 very good, and >0.90 excellent discrimination [47]; Brier score and log loss were interpreted against chance-level (null-model) values, corresponding to predictions equal to the observed outcome prevalence.

Clinical site was included among candidate predictors during the elastic net selection procedure to assess its stability and potential predictive contribution.

3. Results

3.1. Baseline Demographic and Clinical Characteristics

The study cohort comprised women with a mean age of 55.4 years (40-70) recruited across four countries, most commonly Finland (38.1%), followed by Portugal (24.9%), Israel (19.3%) and Italy (17.7%). The majority had a University education (60.7%), were married or living with a partner (74.9%) and were employed (72.9%), with most reporting middle income levels (61.6%). Regarding lifestyle characteristics, two thirds reported engaging in exercise (66.3%), almost half followed a type of diet (45.4%), most consumed alcohol in moderation (68.2%) and the majority were never smokers (67.4%).

Clinically, most participants were postmenopausal (61.5%), had stage I–II disease (91.0%), grade II tumors (52.2%) and ductal histology (77.9%). Tumors were predominantly hormone receptor–positive (ER-positive 89.6%; PR-positive 79.8%) and HER2-negative (81.8%), with luminal A-like (36.9%) and luminal B-like HER2-negative (39.0%) subtypes being most frequent. Breast-conserving surgery was the most common surgical approach (74.6%), radiotherapy was administered to 80.6% of patients and systemic treatment most often consisted of endocrine therapy alone (47.3%) or combined chemotherapy and endocrine therapy (37.7%).

Table 2. Baseline clinical and cancer characteristics of the study participants. Total number of patients n=538.

Variable	n (%)	Variable	n (%)
Negative Life Events		Estrogen receptor Positivity	467 (89.6%)
None	58 (12%)	Progesterone receptor Positivity	410 (79.8%)
One event	239 (49.6%)	HER2 Positivity	89 (18.2%)
Two or more events	185 (38.4%)	Ki67 levels ≥20%	293 (56.7%)
Chronic diseases	191 (35.7%)	Subtypes1
Metabolic diseases		Luminal A-like	175 (36.9%)
Mental illness		Luminal B-like (HER2 -)	185 (39%)
Family history of beast cancer	330 (64.3%)	Luminal B-like (HER2 +)	68 (14.3%)
Menopausal status pre		Her2-positive (non luminal)	20 (4.2%)
Pre/Peri-menopausal	202 (38.5%)	Triple-negative	26 (5.5%)
Postmenopausal	322 (61.5%)	Lumpectomy	391 (74.6%)
HRT before diagnosis	105 (21.6%)	Mastectomy	133 (25.4%)
Cancer stage		Radiotherapy	424 (80.6%)
I	251 (48.2%)	Systemic Therapy
II	223 (42.8%)	Chemotherapy only (± anti-HER2)	78 (14.9%)
III	47 (9%)	Endocrine therapy only	247 (47.3%)
Cancer grade		Chemo + Endocrine therapy (± anti-HER2)	197 (37.7%)
I	91 (17.5%)	Anti-HER2 therapy	82 (15.4%)
II	271 (52.2%)	Neoadjuvant Chemotherapy	84 (16%)
III	157 (30.3%)
Cancer histological type
Ductal	408 (77.9%)
Lobular	80 (15.3%)
Other	36 (6.9%)

¹Luminal A-like: ER+, PR+, HER2-, low Ki67 (<20%); Luminal B-like (HER2 negative): ER+, PR+/-, HER2-, high Ki67 (≥20%) or ER+, PR-, HER2-, Ki67 any; Luminal B-like (HER2 positive): ER+, PR+/-, HER2+, any Ki67; Her2-positive (non luminal) ER-, PR-, HER2+, any Ki67; Triple-negative (ER-, PR-, HER2-, any Ki67) Abbreviations: HRT: Hormone Replacement Therapy; HER2: human epidermal growth factor receptor 2.

Table 2. Baseline clinical and cancer characteristics of the study participants. Total number of patients n=538.

Variable	n (%)	Variable	n (%)
Negative Life Events		Estrogen receptor Positivity	467 (89.6%)
None	58 (12%)	Progesterone receptor Positivity	410 (79.8%)
One event	239 (49.6%)	HER2 Positivity	89 (18.2%)
Two or more events	185 (38.4%)	Ki67 levels ≥20%	293 (56.7%)
Chronic diseases	191 (35.7%)	Subtypes1
Metabolic diseases		Luminal A-like	175 (36.9%)
Mental illness		Luminal B-like (HER2 -)	185 (39%)
Family history of beast cancer	330 (64.3%)	Luminal B-like (HER2 +)	68 (14.3%)
Menopausal status pre		Her2-positive (non luminal)	20 (4.2%)
Pre/Peri-menopausal	202 (38.5%)	Triple-negative	26 (5.5%)
Postmenopausal	322 (61.5%)	Lumpectomy	391 (74.6%)
HRT before diagnosis	105 (21.6%)	Mastectomy	133 (25.4%)
Cancer stage		Radiotherapy	424 (80.6%)
I	251 (48.2%)	Systemic Therapy
II	223 (42.8%)	Chemotherapy only (± anti-HER2)	78 (14.9%)
III	47 (9%)	Endocrine therapy only	247 (47.3%)
Cancer grade		Chemo + Endocrine therapy (± anti-HER2)	197 (37.7%)
I	91 (17.5%)	Anti-HER2 therapy	82 (15.4%)
II	271 (52.2%)	Neoadjuvant Chemotherapy	84 (16%)
III	157 (30.3%)
Cancer histological type
Ductal	408 (77.9%)
Lobular	80 (15.3%)
Other	36 (6.9%)

¹Luminal A-like: ER+, PR+, HER2-, low Ki67 (<20%); Luminal B-like (HER2 negative): ER+, PR+/-, HER2-, high Ki67 (≥20%) or ER+, PR-, HER2-, Ki67 any; Luminal B-like (HER2 positive): ER+, PR+/-, HER2+, any Ki67; Her2-positive (non luminal) ER-, PR-, HER2+, any Ki67; Triple-negative (ER-, PR-, HER2-, any Ki67) Abbreviations: HRT: Hormone Replacement Therapy; HER2: human epidermal growth factor receptor 2.

3.2. Trajectory Groups

3.2.1. GHS/QoL Trajectories

Taking into account statistical fit, class separation, parsimony and clinical interpretability, we selected the five-class LCGA solution as the optimal model (See Appendix B for detailed explanations). Figure 1a shows the resulting trajectories for the five-group model and the actual measurements of the patients that compose each class. An Excellent trajectory (n = 71, 13.2%) demonstrated very high baseline scores (≥90) that remained stable throughout the study period. Mean increase in slope is statistically significant (Table D1 in Appendix D) but overall changes are small (<5–10 points), indicating no clinically meaningful change over time, while individual trajectories cluster tightly at high values with limited variability. The largest class was the Good trajectory (n = 219, 40.7%), characterized by consistently high QoL scores primarily within the good range (≥70). The mean trajectory was essentially flat (Table D1 in Appendix D). This class exhibited greater individual variability than the Excellent class; however, most observations remained ≥70 throughout follow-up. The Moderate class (n = 169, 31.4%) exhibited baseline and follow-up scores predominantly within the moderate range (50–69). Substantial within-class variability was observed, with some individuals showing improvement over time while others remained stable or fluctuated across follow-up. Overall, the mean trajectory remained largely stable (Table D1 in Appendix D). A small Recovering class (n = 42, 7.8%) started with moderately reduced QoL (approximately 60–65) but exhibited a marked and sustained improvement exceeding clinically important change thresholds (≥10 points), reaching levels comparable to the higher-functioning classes by the end of follow-up. Most individuals in this class showed upward trajectories crossing the threshold for clinical importance. Finally, the Low deteriorating class (n = 37, 6.9%) showed low baseline QoL near or below the poor range (<50), with further statistically important (Table D1 in Appendix D) and clinically meaningful deterioration during the early months followed by only partial recovery. The early decline exceeded 10 points, representing clinically meaningful deterioration and subsequent improvement did not fully offset this loss. Although within-class variability was high, most individuals remained below 50–55 for much of the follow-up period.

3.2.2. HADS Depression Trajectories

We selected the four-class GMM solution as the optimal model (See Appendix C for detailed explanations). Figure 1b shows the resulting trajectories for the four-group model and the actual measurements of the patients that compose each class. The largest group was the Resilient trajectory (n = 321, 59.7%), characterized by persistently very low depressive symptom levels throughout follow-up. Within this class, variability was minimal, with overall or individual mean changes being negligible or very small (Table D2 in Appendix D). The Stable moderate/high trajectory (n = 136, 25.3%) exhibited consistently elevated depressive symptoms (mild to moderate range) that remained largely stable over time (Table D2 in Appendix D), indicating a chronic symptom burden. For most patients in this class, depressive symptom scores fluctuated considerably from one assessment to another (>0.2), yet the overall individual trajectories remained stable over time. A smaller Recovering trajectory (n = 54, 10.0%) showed higher baseline depressive symptoms followed by a sustained and clinically meaningful decrease over time (Table D2 in Appendix D), with most individuals improving to very low symptom levels (<0.4). Finally, the Delayed occurrence trajectory (n = 27, 5.0%) demonstrated low baseline depressive symptoms with a gradual, statistically significant (Table D2 in Appendix D) and clinically meaningful increase during follow-up, reaching levels >1, indicative of delayed onset of depressive symptomatology.

3.3. Predictors of C30 GHS/QoL Trajectories

3.3.1. Low Deteriorating QoL vs Rest

In clinical-site–adjusted univariable logistic regression analyses (Supplementary Figure S1), higher odds of membership in the low deteriorating QoL class were observed among unemployed/housewives (OR 2.49), individuals with ≥2 negative life events (OR 5.30), those with chronic disease (OR 1.80) or mental illness (OR 4.92), higher BMI (OR 1.08 per unit), adverse tumor characteristics including stage III disease (OR 2.99), grade III tumors (OR 3.74), Ki-67 ≥20% (OR 2.20), triple-negative subtype (OR 3.51), and receipt of neoadjuvant chemotherapy (OR 4.45). In contrast, middle income (OR 0.35), physical activity (low/moderate: OR 0.22; heavy: OR 0.10), moderate alcohol consumption (OR 0.42), estrogen receptor positivity (OR 0.29), endocrine therapy alone (OR 0.23), and combined chemo-endocrine therapy (OR 0.32) were associated with lower odds.

Both at baseline and month 3, the vast majority of considered psychological, psychosocial, and health-related quality-of-life scales differed significantly between patients in the low deteriorating QoL class and all other classes (Supplementary Figures S2, S3).

Baseline psychological distress was strongly associated with higher odds, including depression (HADS; OR 6.19), anxiety (HADS; OR 3.87), fear of cancer recurrence (FCRI; OR 2.60), distress (NCCN; OR 1.24), and negative affect (PANAS; OR 2.41), other blame (CERQ; OR 2,53), catastrophizing (CERQ; OR 1.45), overall negative cognitive emotion regulation score (CERQ; OR2.76), whereas optimism (OR 0.39), sense of coherence (ORs 0.82–0.89), acceptance (CERQ; OR 0.77) ability to cope with trauma (trauma-focus, forward-focus, flexibility, total PACT; ORs 0.60-0.78), mindfulness (MAAS; OR 0.48), resilience (CDRISC; OR 0.43), coping with cancer (CBI; OR 0.53), self-efficacy (OR 0.73), perceived social support (OR 0.69), and positive affect (PANAS; OR 0.45), put into perspective (CERQ; OR 0.62) were protective. Similarly, better baseline global health status/QoL (C30; OR 0.94), and functioning across almost all domains (physical, role, social, emotional and cognitive C30; ORs 0.95–0.98) were associated with lower odds, while higher financial impact (C30; OR 1.02) and symptom burden (fatigue, pain, arm symptoms, dyspnea, insomnia, appetite loss, constipation, diarrhea and systemic therapy side effects C30 or BR23; ORs 1.01-1.05) were associated with higher odds. Moreover, better body image (BR23; OR 0.98), future perspective (BR23; OR 0.99), and sexual enjoyment (BR23; OR 0.99) were associated with lower odds.

The scales not significantly associated with membership in the low deteriorating QoL class were: polarity (PACT), most CERQ scales (self-blame, rumination, positive refocusing, positive reappraisal and planning; although there is a tendency to significance with the exception of planning and self-blame), as well as nausea (C30), breast symptoms (BR23), sexual function (BR23) and upset by hair loss (BR23).

At month 3, patients in the low deteriorating QoL class showed substantially higher psychological distress compared with all other patients, including depression (HADS; OR 6.97), anxiety (HADS; OR 5.93), distress (NCCN; OR 1.29) negative affect (PANAS; OR 2.44), anxious preoccupation (MAC; OR 3.54) and helplessness (MAC; OR 3.92). Conversely, positive affect (PANAS; OR 0.37), perceived social support (OR 0.78), general self-efficacy (OR 0.71), adaptive family communication and cohesion (FARE; ORs 0.56–0.58), personal control beliefs (OR 0.85), and positive treatment beliefs (OR 0.67) were associated with lower odds. Among coping behaviors, exercising (OR 0.62) and looking at positive aspects (OR 0.76) were protective.

Health-related QoL at month 3 similarly differentiated trajectories. Better global health status/QoL (OR 0.95) and higher functioning (physical, role, emotional, cognitive, and social functioning C30; ORs 0.94–0.97) were associated with lower odds of low deteriorating QoL, whereas greater C30 symptom burden—including fatigue (OR 1.04), pain (OR 1.03), dyspnea (OR 1.03), insomnia (OR 1.02), appetite loss (OR 1.02), diarrhea (OR 1.02), financial impact (OR 1.02), and treatment-related side effects (OR 1.05)—was associated with higher odds. Breast cancer–specific domains further characterized this group, including poorer body image (OR 0.98), more breast and arm symptoms (ORs 1.03), reduced future perspective (OR 0.98), and lower sexual functioning and enjoyment (ORs 0.97–0.98).

Variables related to post-traumatic growth (PTGI; relating to others, new possibilities, personal strength, spiritual change, appreciation of life, total score), mental adjustment styles (MAC; fighting spirit, avoidance, fatalism, ) and certain copying behaviors (tried to relax, distracted yourself, prayed, etc.) showed no statistically significant difference between the groups.

To identify robust features associated with the low deteriorating QoL profile, we implemented a feature selection pipeline based on elastic net logistic regression across 30 multiply imputed datasets under penalized site and unpenalized site conditions. Predictor importance was quantified using selection frequency. Predictive model performance was evaluated using cross-validated predictions from elastic net logistic regression models.

A limited set of features emerged as robust in distinguishing the low-deteriorating class from the remaining patients (Table 3). At baseline, depressive symptoms, emotional functioning, global health status/quality of life (GHS/QoL), fatigue, pain, diarrhea, sense of manageability, coping with cancer, perceived social support, and other-blame coping were the most stable predictors, each selected with 100% or near-100% (≥97%) frequency under both penalized and unpenalized clinical site conditions. Receipt of neoadjuvant chemotherapy and triple-negative molecular profile followed, also demonstrating high selection stability. Negative life events showed borderline selection frequency, suggesting a weaker and less consistent contribution to the model.

At month 3, cognitive functioning, depressive symptoms, physical functioning, and treatment control beliefs emerged as the most robust predictors, each being selected in 100% of the 90 stability iterations under both penalized and unpenalized clinical site conditions. These were followed by anxiety symptoms and family communication and cohesion. In contrast, neoadjuvant chemotherapy showed unstable selection when the clinical site was forced into the model.

Both at baseline and month 3, the models demonstrated strong discriminative ability (ROC-AUC ~ 0.86) and high probabilistic accuracy (for chance performance log-loss=0.253 and Brier score=0.064) (Table 3).

3.3.2. Excellent QoL vs Rest

Compared with all other classes (Supplementary Figure S4), excellent QoL was associated with older age (OR 1.03), heavy physical activity (OR 2.60), moderate alcohol consumption (OR 2.24), lower BMI (OR 0.93), no negative life events (ORs 0.36–0.17), and absence of mental illness (OR 0.09). Favorable disease characteristics, including lower stage (ORs 0.58–0.20), lower grade (ORs 0.50–0.34), low Ki-67 (OR 0.53) and luminal A-like subtype (OR 2.13), were also associated with excellent QoL. In contrast, neoadjuvant chemotherapy was inversely associated (OR 0.15), whereas endocrine therapy alone was positively associated (OR 2.67) with excellent QoL.

At baseline and month 3, almost all considered psychological, psychosocial, and health-related quality-of-life scales differed significantly between patients in the excellent QoL class and all other classes, with consistent advantages observed across distress, coping, affect, functioning, and symptom domains (Supplementary Figures S5-S6).

At baseline, the Excellent QoL group was less likely to report negative psychological traits, including depressive symptoms (HADS; OR = 0.06), anxiety (HADS; OR = 0.14), fear of recurrence (FCRI; OR = 0.42), self-blame (CERQ; OR = 0.53), other-blame (CERQ; OR = 0.59), catastrophizing (CERQ; OR = 0.38), and overall negative cognitive emotion regulation strategies (CERQ; OR = 0.29).

Conversely, the Excellent QoL group was more likely to report positive psychological traits, including optimism (LOT; OR = 2.36); comprehensibility (SOC; OR = 1.16), manageability (SOC; OR = 1.15), and meaningfulness (SOC; OR = 1.18); positive coping strategies (PACT: forward focus OR = 1.56, trauma focus OR = 1.72, total coping OR = 1.36); adaptive emotion regulation strategies (CERQ: perspective-taking OR = 1.43, positive refocusing OR = 1.57, positive reappraisal OR = 1.37, acceptance OR = 1.41); mindfulness (MAAS; OR = 2.65); resilience (CDRISC; OR = 3.89); and positive affect (PANAS; OR = 2.97).

Moreover, better baseline global health status/QoL (C30; OR = 1.14) and better functioning across all C30 domains, including physical (OR = 1.09), role (OR = 1.05), emotional (OR = 1.07), cognitive (OR = 1.06), and social functioning (OR = 1.05), were associated with higher odds of belonging to the Excellent QoL class. In contrast, greater symptom burden, including fatigue (OR = 0.95), nausea (OR = 0.92), pain (OR = 0.95), dyspnea (OR = 0.97), insomnia (OR = 0.98), appetite loss (OR = 0.97), constipation (OR = 0.97), and diarrhea (OR = 0.96) (C30), as well as systemic therapy side effects (BR23; OR = 0.92), breast symptoms (OR = 0.98), and arm symptoms (OR = 0.96), was associated with lower odds of Excellent QoL. Moreover, better body image (BR23; OR = 1.04), future perspective (BR23; OR = 1.03), sexual function (BR23; OR = 1.01), and sexual enjoyment (BR23; OR = 1.01) were associated with higher odds of Excellent QoL.

Rumination (CERQ), planning (CERQ), polarity (PACT) and upset hair image (BR23) did not show a statistically significant difference between the two groups.

At month 3, individuals in the Excellent QoL group showed significant differences across a wide range of psychosocial and coping variables. They were significantly more likely to report strong family communication and cohesion (FARE; OR = 2.43), effective family coping (FARE; OR = 2.17), higher general self-efficacy (single item; OR = 1.76), better adherence to medical advice (OR = 1.56), greater perceived social support (single item; OR = 1.52), instrumental support (mMOS; OR = 1.53), emotional support (mMOS; OR = 2.12), and overall social support (mMOS; OR = 2.01). They also reported higher positive affect (PANAS; OR = 5.46), stronger personal control beliefs (OR = 1.28), treatment control beliefs (OR = 1.36), and greater use of exercise as a coping behaviour (OR = 1.34).

Conversely, the Excellent QoL group was significantly less likely to report depressive symptoms (HADS; OR = 0.03), anxiety (HADS; OR = 0.07), overall mental health distress (HADS total; OR = 0.02), general distress (NCCN Distress Thermometer; OR = 0.62), helplessness (MAC; OR = 0.15), anxious preoccupation (MAC; OR = 0.21), avoidance coping (MAC; OR = 0.56), negative affect (PANAS; OR = 0.15), and the coping behaviors crying (OR = 0.47), talking to the physician (OR = 0.62) or ask for help (OR 0.79).

All health-related QoL domains at Month 3 differentiated between trajectory groups. Better global health status/QoL and higher overall functioning were associated with a greater likelihood of belonging to the Excellent QoL class (ORs 1.05–1.16), whereas greater symptom burden and financial impact were associated with lower odds (ORs 0.93–0.98). Breast cancer–specific domains further contributed to this differentiation, with better body image, future perspective, and sexual functioning (ORs 1.01–1.04) and fewer treatment-related, breast, and arm symptoms (ORs 0.94–0.98) characterizing the Excellent QoL group.

Variables related to post-traumatic growth (PTGI; relating to others, new possibilities, personal strength, spiritual change, appreciation of life, and total PTGI score), specific mental adjustment styles (MAC; fighting spirit and fatalism), and certain coping behaviors (e.g., trying to relax, distraction, praying, and looking at positive sides and perceiving the situation as a challenge) showed no statistically significant differences between the groups.

Stability selection at baseline and month 3 identified a relatively broad set of factors distinguishing the excellent QoL class (Table 4).

At baseline, robust predictors were associated with psychological resources and lower distress, including lower anxiety and higher mindfulness, resilience, coping, self-efficacy, and perceived social support. Better global health status and functioning across physical, role, emotional, and cognitive domains, together with lower catastrophizing and symptom burden, further differentiated the excellent QoL class. In addition, treatment-related factors (receipt of endocrine therapy and absence of neoadjuvant chemotherapy), favorable molecular phenotypes (more Luminal A–like and less Luminal B–like [HER2+]), and selected sociodemographic characteristics (not being unemployed/housewife or following a vegetarian diet) contributed to class differentiation. Notably, the stability of pre-existing mental illness and positive affect dropped well below 60% when clinical site was forced into the model.

At month 3, robust predictors reflected lower psychological distress and treatment-related symptom impact, including lower anxiety, anxious preoccupation, fatigue, systemic therapy side effects, and distress, alongside greater positive affect, future perspective, and role and social functioning. These variables were selected in 100% or near-100% of stability iterations. Additional contributors included better physical and emotional functioning, stronger personal control beliefs, greater perceived social and emotional support, lower pain and arm symptoms, fewer depressive symptoms, better family communication and cohesion, and reduced exposure to multiple negative life events. In contrast, the stability of negative affect declined markedly when clinical site effects were accounted for.

Both at baseline and month 3, the models demonstrated strong discriminative ability and high probabilistic accuracy (for chance performance log-loss=0.390 and Brier score=0.115) (Table 4).

3.3.3. Recovery vs Moderate QoL

Among sociodemographic, lifestyle, clinical, and cancer-related factors, clinical-site–adjusted univariable logistic regression analyses (Supplementary Figure S7) identified exposure to negative life events as the only factor significantly associated with QoL recovery, with lower odds observed across increasing exposure levels (OR range = 0.16–0.29).

Several baseline (M0) factors were significantly associated with QoL recovery. Protective factors included higher optimism (LOT; OR = 2.72), greater resilience (CD-RISC; OR = 2.42), greater coping flexibility (PACT; OR range = 1.23–1.52), and greater use of adaptive cognitive emotion regulation strategies, such as perspective-taking (CERQ; OR = 1.57) and planning (CERQ; OR = 1.72), all of which were associated with increased odds of recovery. In contrast, risk factors included higher depressive symptoms (HADS; OR = 0.40), anxiety symptoms (HADS; OR = 0.46), and greater use of catastrophizing (CERQ; OR = 0.63), which were associated with lower odds of QoL recovery (Supplementary Figure S8). Regarding health-related QoL domains at Month 3, better global health status/QoL (C30; OR = 1.03) and higher functioning, including physical (OR = 1.03), role (OR = 1.02), emotional (OR = 1.02), and social functioning (all C30; OR = 1.03), were associated with a greater likelihood of belonging to the Recovery QoL class. In contrast, higher symptom burden, particularly fatigue (C30; OR = 0.97) and pain (C30; OR = 0.97), was associated with lower odds of recovery. Among breast cancer–specific domains, better sexual functioning (EORTC QLQ-BR23; OR = 1.02) was positively associated with QoL recovery, whereas most symptom and body image domains did not significantly differentiate recovery trajectories.

The psychological scales not significantly associated with membership in the Recovery QoL class included sense of comprehensibility, manageability, and meaningfulness (SOC); fear of recurrence (FCRI); trauma and total coping (PACT), with polarity coping showing only a borderline association; most cognitive emotion regulation strategies (CERQ), including self-blame, other-blame, rumination, positive reappraisal, acceptance, and negative overall CERQ (with borderline effects for positive refocusing and acceptance); distress (NCCN Distress Thermometer); and negative affect (PANAS).

At Month 3, site-adjusted models indicated that lower symptom burden and more adaptive coping remained associated with QoL recovery. Specifically, lower helplessness (MAC; OR = 0.13), lower anxious preoccupation (MAC; OR = 0.34), lower anxiety and depressive symptoms (HADS; ORs=0.12-0.21), lower distress (NCCN; OR = 0.86), lower negative affect (PANAS; OR=0.52), higher fighting spirit (MAC; OR = 3.93), higher positive affect (PANAS; OR = 3.46), greater general self-efficacy (OR = 1.74), stronger personal and treatment control beliefs (ORs = 1.39–1.78), higher family coping (FARE; OR=1.82) and greater emotional support (mMOS; OR = 2.01) were significant predictors of recovery. In addition, better global health status/QoL and higher functioning (physical, role, emotional, cognitive, and social; C30; ORs = 1.02–1.05), lower pain (C30; OR = 0.97), fatigue and nausea (C30; ORs = 0.97–0.98), and better sexual functioning, enjoyment, and future perspective (BR23; ORs = 1.02–1.03) were also associated with increased odds of QoL improvement (Supplementary Figure S9).

At Month 3, the scales that did not differ between the groups were post-traumatic growth (PTGI; relating to others, new possibilities, personal strength, spiritual change, appreciation of life, and total PTGI score), mental adjustment styles reflecting avoidance and fatalism (MAC), family communication and cohesion (FARE), adherence to medical advice, several coping behaviors (including trying to relax, distraction, praying/going to church, exercising, bursting into tears, talking to or asking help from someone important, and talking to the physician), and instrumental social support (mMOS). In terms of health-related QoL, dyspnea, insomnia, constipation, diarrhea, and financial impact (EORTC QLQ-C30), as well as breast symptoms, arm symptoms, and upset by hair loss (EORTC QLQ-BR23), did not significantly differ between groups.

A concise set of features emerged as robust in distinguishing the recovery class from the moderate QoL one (Table 5). At baseline, coping with cancer, mindfulness, optimism, perspective-taking emotion regulation, pain, positive affect, sexual functioning, and social functioning, together with middle and high income (vs low income), emerged as the most robust predictors, had a selection frequency of 100% or above 90%. These were followed by planning as a cognitive emotion regulation strategy. The stability of resilience, exposure to two or more negative life events, adherence to a special diet dropped markedly when clinical site was forced into the model. No psychological distress symptoms emerged as robust predictors.

At month 3, a similarly concise set of features robustly distinguished the recovery QoL class from the rest of the patients (Table 5). Helplessness, pain, positive affect, sexual functioning, personal control beliefs over illness, middle and high income (vs low income), coping by perceiving the situation as a challenge and general self-efficacy were the most stable predictors, each selected in 100% or nearly 100% of the stability iterations, even when clinical site was forced in the model. These were followed by social functioning, fighting spirit, depression and anxiety symptoms. The selection frequency of postmenopausal status, negative life events and and triple-negative disease dropped below the threshold when clinical site was forced into the model.

At baseline, the model showed fair discrimination and moderate probabilistic accuracy (Table 5; for chance performance log-loss=0.500 and Brier score=0.159). At month 3, performance improved, with good discrimination and better calibration and accuracy, indicating enhanced predictive ability over time. Together, these results suggest that incorporating baseline and month 3 information would substantially enhances the model’s ability to predict QoL trajectories compared with baseline alone.

3.4. Predictors of HADS Depression Trajectories

3.4.1. Stable Moderate/High vs Resilient

Compared with the Resilient class, several protective factors were associated with lower odds of belonging to the Stable Moderate/High Depression class, including residence in Finland (vs Portugal) (OR = 0.43), university education (OR = 0.52), being married or in a common-law relationship (vs. single) (OR = 0.47), higher income (OR = 0.38), postmenopausal status (OR = 0.59), engagement in physical activity (heavy: OR = 0.32), and moderate alcohol consumption (OR = 0.49). In contrast, risk factors associated with higher odds of stable moderate/high depression included residence in Italy (vs Portugal) (OR = 3.56), unemployment or being a housewife (OR = 2.59), exposure to two or more negative life events (OR = 2.38), history of mental illness (OR = 3.99), HER2-positive (non-luminal) tumor subtype (OR = 3.59), and receipt of neoadjuvant chemotherapy (OR = 2.20) (Supplementary Figure S10).

At baseline, the scales that did not show statistically significant differences between the Resilient and Stable Moderate/High Depression groups were self-blame (CERQ) and acceptance (CERQ), with other-blame (CERQ) showing only a borderline association. All other psychological, coping, and affective scales, as well as all health-related quality-of-life domains were significantly associated with group membership (Supplementary Figure S11).

At Month 3, the scales not significantly associated with group membership between the Resilient and Stable Moderate/High Depression classes were post-traumatic growth dimensions (PTGI; relating to others, new possibilities, personal strength, appreciation of life, and total PTGI score), with the exception of spiritual change, fatalism (MAC), and few coping behaviors (trying to relax, praying/going to church, talking to or asking for help from somebody important). In terms of health-related QoL, only diarrhea (C30) did not significantly differ between the groups. All other psychological, coping, social support, and QoL scales demonstrated statistically significant associations (Supplementary Figure S12).

The feature selection procedure identified a concise set of variables that robustly distinguished the stable moderate/high depression class from the resilient class (Table 6). At baseline, risk factors consistently selected across all stability iterations included higher anxiety and depressive symptoms, greater arm symptoms, higher financial impact, poorer future perspective, greater catastrophizing, lower sense of manageability and meaningfulness, lower optimism, poorer role functioning, and residence in Italy or Finland (vs Portugal). Additional variables distinguishing the stable moderate/high class were higher distress and lower resilience. Stability of coping with cancer, heavy exercise (vs none), unemployment or being a housewife (vs employment) dropped to 0% when clinical site was forced into the model.

At month 3, variables consistently selected across all stability iterations and penalized clinical site conditions included anxiety, distress, anxious preoccupation, negative affect, helplessness, emotional functioning, positive affect, greater emotional support and residence in Italy or Finland (vs Portugal). Additional selected variables included fatigue, pain, arm symptoms, spiritual change, sexual functioning and enjoyment and cognitive functioning. In contrast, the selection frequency of avoidance coping, radiotherapy, university education, use of exercise as a coping strategy and heavy physical activity dropped considerably, even to 0 when clinical site was forced to the model.

At baseline, the model demonstrated excellent discrimination and calibration performance, while at month 3, performance remained high, though slightly reduced, indicating robust predictive ability at both time points (Table 6; for chance performance log-loss=0.610 and Brier score=0.209).

3.4.2. Delayed occurrence vs Resilient

Regarding sociodemographic, lifestyle and clinical data, compared with the Resilient class, several risk factors were associated with a higher likelihood of belonging to the Delayed Occurrence Depression class. These included history of mental illness (OR = 42.84), HER2-positive tumors (OR = 2.99), triple-negative breast cancer subtype (OR = 3.63), and receipt of anti-HER2 therapy (OR = 2.91). Additional risk was observed for individuals with HER2-positive (non-luminal) tumors (OR = 4.80). In contrast, several protective factors were associated with lower odds of delayed depression, including residence in Finland (vs Portugal) (OR = 0.15), university education (OR = 0.28), middle and high income (vs low income) (ORs = 0.22–0.25), engagement in physical activity (low–moderate: OR = 0.32; heavy: OR = 0.20), never or former smoking (ORs = 0.21–0.31), and estrogen receptor–positive disease (OR = 0.29).

At baseline, poorer psychological resources and quality of life were significantly associated with higher odds of belonging to the Delayed Occurrence Depression class. Specifically, lower optimism (LOT; OR = 0.34), lower sense of coherence, including manageability (SOC; OR = 0.86) and meaningfulness (SOC; OR = 0.85), lower resilience (CD-RISC; OR = 0.36), and reduced use of positive cognitive emotion regulation strategies, including positive refocusing (CERQ; OR = 0.64) and lower overall positive CERQ (OR = 0.56), were associated with increased risk. In terms of health-related QoL, worse global health status/QoL (EORTC QLQ-C30; OR = 0.97), poorer physical and role functioning (ORs = 0.96–0.97), and greater symptom burden, namely fatigue (OR = 1.03), pain (OR = 1.05), insomnia (OR = 1.02), appetite loss (OR = 1.03), diarrhea (OR = 1.04), and financial impact (OR = 1.03), were also significantly associated with delayed depression, along with greater arm symptoms (BR23; OR = 1.04).

At baseline, core affective symptoms were not significantly associated with delayed depression, including depressive symptoms, anxiety symptoms, and overall mental health distress (HADS). Several coping styles and emotion regulation strategies also showed no significant differences, including forward-, trauma-, total, polarity, and flexibility coping (PACT), self-blame, other-blame, rumination, catastrophizing, perspective-taking, acceptance, planning, and overall negative CERQ, as well as mindfulness (MAAS), coping with cancer (CBI), general self-efficacy, and perceived social support. In addition, fear of cancer recurrence (FCRI), distress thermometer scores (NCCN), negative affect (PANAS), and most breast cancer–specific QoL domains, including body image, treatment side effects, breast symptoms, sexual function, sexual enjoyment, and upset by hair loss (EORTC QLQ-BR23), were not significantly associated. Among EORTC QLQ-C30 domains, cognitive functioning, dyspnea, nausea, and constipation did not differ significantly between groups.

At month 3, psychological distress and affective symptoms were strongly associated with higher odds of belonging to the Delayed Occurrence of Depression class, including depressive symptoms (HADS; OR = 9.50), anxiety symptoms (HADS; OR = 5.40), overall mental health distress (HADS total; OR = 17.51), and negative affect (PANAS; OR = 2.13). In contrast, lower positive affect (PANAS; OR = 0.34) and lower emotional and total social support (mMOS; ORs = 0.56–0.56) were associated with reduced odds of delayed depression. Maladaptive coping and behavioral responses were also relevant, with coping by bursting into tears (OR = 1.57) and talking to the physician (OR = 1.69) associated with higher odds.

In terms of health-related QoL, poorer global health status/QoL (EORTC QLQ-C30; OR = 0.97), lower physical, role, emotional, and cognitive functioning (ORs = 0.96–0.98), and greater symptom burden, including nausea (OR = 1.03), pain (OR = 1.02), appetite loss (OR = 1.02), diarrhea (OR = 1.03), and financial impact (OR = 1.02), were significantly associated with delayed depression. Breast cancer–specific domains further characterized this group, including greater arm symptoms (BR23; OR = 1.03) and poorer future perspective, sexual functioning, and sexual enjoyment (BR23; ORs = 0.97–0.98).

At month 3, several domains did not significantly differentiate the delayed depression and resilient groups, including post-traumatic growth (PTGI), most mental adjustment styles (MAC), family functioning (FARE; although it tends to significance), and general coping resources and behaviors (adherence to medical advice, general self-efficacy, perceived support, and most specific coping strategies). In terms of health-related QoL, social functioning, fatigue, dyspnea, insomnia, and constipation (EORTC QLQ-C30), as well as breast cancer–specific domains such as treatment side effects, breast symptoms, and upset by hair loss (BR23), were not significantly associated.

At baseline, a limited set of features emerged from the stability selection pipeline as robust discriminators between the delayed depression occurrence and resilient groups (Table 7). Among these, only diarrhea, pain, and role functioning demonstrated robust selection, each with 100% selection frequency under both penalized and unpenalized clinical site conditions. In contrast, the selection stability of sense of manageability and optimism declined substantially when clinical site was forced into the model.

At month 3, a broader range of factors differentiated the delayed depression trajectory. Variables showing high selection frequency under both penalized and unpenalized site conditions included diarrhea, emotional functioning, anxiety symptoms, pre-existing mental illness, university education, and coping through communication with the physician, followed by sexual functioning and triple-negative disease. Middle income (vs. low income) showed borderline predictive value, with selection frequency decreasing to approximately 40% after forcing clinical site into the model. Similarly, the stability of unemployment/housewife status and heavy exercise declined markedly once clinical site was unpenalized.

At baseline, the model showed good discrimination and calibration, while at month 3, performance remained stable, with similar probabilistic accuracy and a slightly lower discrimination (Table 7; for chance performance log-loss=0.273 and Brier score=0.072).

3.4.3. Recovery vs Stable Moderate/High

Compared with the Stable Moderate/High Depression class, predictors of belonging to the Recovery Depression class included high income (vs low income; OR = 3.18), HER2-positive disease (OR = 2.56), and luminal B–like (HER2-positive) tumor subtype (OR = 3.94) (Supplementary Figure S16).

At baseline (M0), predictors of belonging to the Recovery Depression class included higher optimism (OR = 2.44), greater sense of coherence–manageability (OR = 1.15), greater sense of coherence–meaningfulness (OR = 1.13), higher mindfulness (OR = 1.87), higher positive affect (OR = 1.90), and lower pain (OR = 0.98) (Supplementary Figure S17).

At Month 3 (M3), predictors of belonging to the Recovery Depression class included lower depressive symptoms (OR = 0.24), lower anxiety symptoms (OR = 0.17), better overall mental health (OR = 0.13), lower spiritual change (PTGI; OR = 0.70), lower helplessness (MAC; OR = 0.36), lower anxious preoccupation (MAC; OR = 0.47), lower distress (NCCN Distress Thermometer; OR = 0.86), higher perceived social support (OR = 1.31), greater treatment control beliefs (OR = 1.28), higher positive affect (OR = 2.19), lower negative affect (OR = 0.49), greater instrumental support (OR = 1.54), greater emotional support (OR = 2.04), greater total social support (OR = 1.90), better emotional functioning (OR = 1.04), better cognitive functioning (OR = 1.02), lower pain (OR = 0.98), lower insomnia (OR = 0.98), fewer arm symptoms (OR = 0.97), and a more positive future perspective (EORTC QLQ-BR23; OR = 1.01) (Supplementary Figure S18).

Based on variable selection method, retained variables at baseline when clinical site was penalized were: sense of manageability, optimism, clinical site (Italy vs Portugal), treatment type (endocrine therapy only vs chemotherapy ± anti-HER2), and high income (vs low income), all selected in 100% of the stability iterations, followed by residence in Finland (vs Portugal), selected in 60% of iterations.

At month 3, variables selected in 100% of the stability iterations under penalized site conditions included anxiety symptoms, clinical site (Italy vs Portugal), treatment type (endocrine therapy only vs chemotherapy ± anti-HER2), and high income (vs low income). Variables with high but slightly lower stability included spiritual change, following a special diet, and coping by talking to someone important (each selected in 90% of iterations), as well as emotional functioning. Additional contributors included upset by hair loss, metabolic diseases, residence in Finland (vs Portugal), negative affect, and emotional support (Table 8).

While optimism, income (baseline), endocrine only treatment (baseline and month 3), coping by talking to someone important (month 3), upset by hair image (month 3 and negative affect (month 3) initially demonstrated high selection stability (in the clinical site-penalized model, their selection frequencies dropped dramatically, when site was unpenalized (i.e. forced into the model), indicating that these predictors were site-dependent. In other words, if the models know which clinical site a patient is at, knowing their income or optimism level adds almost no extra predictive value to the specific model.

The selection frequency for Spiritual change (month 3) decreased from 90% (penalized site) to 40% (unpenalized site) and for high income (month 3) from 100% (penalized site) to 57% (unpenalized site). This indicates that while site variation explains a significant portion of the effect of these variables, a moderate, site-independent effect remains, suggesting this variable possesses a more robust, universal association with the outcome compared to other site-dependent variables.

Although stability analysis identified a consistent set of predictors, overall model discrimination at baseline and month 3 was modest (Table 8; for chance performance log-loss=0.596 and Brier score=0.204), indicating limited predictive capacity. Accordingly, stability-selected variables should be interpreted as robust correlates of the recovery Depression profile.

4. Discussion

Using latent class growth analysis (LCGA) and growth mixture modeling (GMM), we identified distinct trajectories of C30 Global Health Status/Quality of Life and HADS Depression among women with early breast cancer over an 18-month period following baseline, highlighting clinically meaningful heterogeneity in longitudinal patient-reported outcomes.

Five distinct trajectories of C30 Global Health Status/Quality of Life were identified. Two trajectories reflected stable high functioning, comprising the Good (40.7%) and Excellent (13.2%) groups, whereas nearly one-third of participants belonged to a Moderate trajectory (31.4%) characterized by persistently intermediate QoL. A small Recovering group (7.8%) demonstrated clear and clinically meaningful improvement over time. Finally, a small but clinically vulnerable subgroup (Low deteriorating, 6.9%) experienced sustained poor QoL with early deterioration and limited recovery.

Park et al. (2023) [48] identified three largely stable QoL trajectories over the first year following the end of primary treatment among breast cancer survivors (N = 124). The majority of women maintained relatively good QoL over time, while a smaller subgroup reported persistently low QoL following treatment. Di Meglio et al. (2022) [49] identified distinct long-term QoL trajectories among women with stage I–III breast cancer treated with adjuvant chemotherapy (N = 4,131), including stable excellent and very good patterns, as well as notably smaller persistently poor and deteriorating trajectories extending up to four years after diagnosis. Goyal et al. (2018) [50] examined quality of life over an 18-month period following baseline assessment (conducted within eight months after diagnosis) among women with newly diagnosed stage I–III breast cancer (N = 653) and identified six trajectories. These included persistently low or very low QoL trajectories, trajectories characterized by moderate or high QoL, as well as two improving trajectories leading to moderate or high QoL, similar to our study.

In the case of depressive symptoms, we identified four distinct trajectories, underscoring clinically meaningful heterogeneity in the longitudinal course of depression. The majority of participants followed a Resilient trajectory (59.7%), characterized by persistently low depressive symptom levels. A substantial proportion belonged to a Stable moderate/high trajectory (25.3%), marked by consistently elevated symptoms indicative of a chronic depressive burden. A smaller Recovering group (10.0%) demonstrated clear and clinically meaningful improvement over time, with symptoms declining to low levels. Finally, a small but clinically vulnerable Delayed occurrence subgroup (5.0%) experienced a gradual and clinically meaningful increase in depressive symptoms during follow-up, reflecting delayed onset of depression.

In a cohort of 4,803 women with stage I–III breast cancer, Charles et al. (2022) [51] examined trajectories of depressive symptoms measured using the Hospital Anxiety and Depression Scale (HADS) over the three years following diagnosis. Six distinct trajectory groups were identified, ranging from persistent non-cases to stable depression. Remission and delayed-onset patterns were identified that comprised only a small proportion of patients, consistent with the findings of the present study. Kant et al. (2018) [52] identified four distress trajectories among 181 newly diagnosed breast cancer patients over the first six months following surgery. A resilient trajectory comprised the majority of patients, whereas high-remitting, delayed, and chronic distress trajectories accounted for smaller proportions, similar to the pattern observed in our study. In 300 Chinese women with breast cancer, Li et al. (2022) [53] reported stable none/mild, stable low and high-decreasing depressive symptom trajectories over 18 months after discharge.

Across these studies, similar to our findings, stable good-QoL or depression trajectories typically comprised the largest proportion of patients, whereas deteriorating or persistently poor QoL trajectories consistently represented smaller but clinically vulnerable subgroups. Differences across studies in the cohort size, baseline timing (before or after the end of primary treatment), study duration, assessment intervals and modelling approach and model selection criteria likely explain the observed variation in the number and type of trajectories reported in the literature, as well as compared to our findings.

To identify parsimonious sets of predictors that jointly differentiate longitudinal quality-of-life and depression trajectory groups, we combined clinically interpretable logistic regression models with multivariable feature selection using elastic net regularization. Following trajectory clustering, individual trajectory classes were retained either to capture nuanced patterns of change in quality of life and depressive symptoms over time or to conduct clinically motivated binary contrasts that support risk stratification and interpretability. For quality of life, comparisons included low-deteriorating versus remaining trajectories, excellent versus remaining trajectories, and recovering versus moderate trajectories, the latter with comparable baseline GHS/QoL levels. For depression, analyses contrasted the resilient class with the stable moderate/high depression class, the delayed depression occurrence class with the resilient class, and the recovery class with the stable moderate/high depression class. Importantly, all analyses were conducted under two conditions, with the clinical site variable alternatively penalized and left unpenalized, to assess the robustness of predictor selection to potential site effects. Elastic net-penalized logistic regression enabled simultaneous variable selection and coefficient shrinkage, allowing identification of robust predictors while accounting for multicollinearity among candidate psychosocial, clinical, and sociodemographic factors.

Based on clinical-site adjusted univariate logistic regression, Low deteriorating QoL and excellent QoL are largely, but not perfectly mirror images.

Regarding low deteriorating QoL trajectory class, both at baseline and month 3 elastic net feature selection consistently prioritized a short list of variables. At baseline these reflect psychological distress (depressive symptoms, emotional functioning), specific physical symptoms (fatigue, pain, diarrhea) , and psychosocial resources (other blame, perceived support, manageability and coping with cancer, exposure to negative life events) irrespective of clinical site penalization, highlighting their robust and site-independent contribution. Clinical factors namely neoadjuvant chemotherapy and triple-negative molecular profile showed an importance however more moderate prediction power. At month 3, the selected predictors were related to functional status (cognitive and physical functioning), psychological symptoms (depressive and anxiety symptoms), believes about the control of treatment over the disease and family interpersonal relationships.

Regarding the excellent QoL trajectory, the analysis identified a broad set of robust predictors, strongly suggesting that stronger internal and interpersonal resources, proactive coping strategies, and lower mental distress are closely associated with achieving an excellent quality of life. Moreover, membership in the excellent QoL class reflected a global pattern of sustained superior functioning and lower symptom burden across all QoL domains, evident at both baseline and month 3. In addition, treatment-related factors—including receipt of endocrine therapy and absence of neoadjuvant chemotherapy—together with favorable molecular phenotypes (e.g. Luminal A–like) and selected sociodemographic characteristics (not being unemployed/housewife and not following a vegetarian diet) further contributed to differentiation of the excellent QoL class.

Together, these findings suggest that preventing deterioration is not simply the inverse of promoting excellence. Accordingly, interventions may require distinct emphases and targets, depending on whether the clinical goal is risk mitigation or the promotion of flourishing and optimal quality of life.

Recovery from a moderate QoL trajectory was characterized by a focused profile of adaptive psychological resources, positive affect, and preserved social and sexual functioning, rather than by low psychological distress. At both baseline and month 3, recovery was consistently associated with active and positive coping styles, optimism, mindfulness, personal control beliefs, and self-efficacy, alongside lower pain and better functioning. Socioeconomic factors (namely income) also played a role. Together, these findings suggest that QoL recovery reflects the mobilization of adaptive psychological and functional resources, highlighting targets that may be particularly relevant for interventions aimed at restoring well-being rather than merely preventing decline.

The contrast between the resilient and stable moderate/high depression trajectories was primarily defined by a consistent pattern of psychological vulnerability, functional impairment, and symptom burden. Across baseline and month 3, persistent depression was associated with higher anxiety, distress, negative affect, helplessness, and catastrophizing, alongside poorer emotional, role, and cognitive functioning and greater physical symptoms. In contrast, resilience appeared to be characterized by preserved meaning, optimism, manageability, and emotional functioning, rather than the absence of clinical or treatment-related exposures.

The comparison between the resilient and delayed depression occurrence trajectories suggests that delayed onset of depressive symptoms is less strongly rooted in baseline psychological vulnerability. At baseline, differences were minimal and primarily reflected physical symptoms and role functioning, whereas by month 3, delayed depression was characterized by worsening emotional functioning, increased anxiety, and greater engagement with illness-related concerns. Delayed depression may reflect a dynamic response to accumulating treatment- and disease-related stressors, rather than pre-existing psychological risk alone.

The contrast between the recovery and stable moderate/high depression trajectories revealed only a limited set of robust predictors once clinical site effects were taken into account. At baseline, sense of manageability emerged as the sole stable psychological factor associated with recovery. By month 3, recovery was primarily differentiated by lower anxiety, better emotional functioning, spiritual change, emotional support, income, and metabolic disease. The dominance of clinical site in these models suggests that country-level contextual factors, such as differences in healthcare systems, cultural norms, social support structures, and treatment pathways, account for a large proportion of the variance between recovery and persistent depression.

The current study extends the findings of a previous analysis by Karademas et.al. (2023) [54] stemming from BOUNCE dataset, which employed a shape-based clustering approach (advanced k-means using kmlShape in R [55]) combined with a random forest classifier to identify predictors of dichotomized outcomes. By adopting a more statistically rigorous probabilistic framework, namely Latent Class Growth Analysis (LCGA) and Growth Mixture Modeling (GMM), we moved from distance-based partitioning to model-based clustering, which explicitly accounts for the underlying data distribution and provides a formal basis for class selection through established fit indices. Furthermore, rather than collapsing trajectories into binary endpoints, a practice that can obscure clinically relevant temporal heterogeneity, we retained individual trajectory classes to capture nuanced patterns of change in quality of life and depressive symptoms over time, while conducting clinically motivated binary contrasts (e.g., excellent vs. non-excellent or low-deteriorating vs. remaining trajectories) to support risk stratification analyses. This strategy, combined with Elastic Net regularization, enabled a more stable and interpretable selection of predictors by effectively addressing multicollinearity and mitigating the risk of overfitting and unstable variable selection that can arise in non-parametric tree-based methods, particularly in settings with highly correlated predictors [45,56]. Finally, clinical site heterogeneity was explicitly addressed through secondary stability analyses in which site variables were forced into the model, allowing for the identification of robust, site-independent predictors.

The latent class trajectories identified in this study show strong conceptual alignment with the patterns reported in [54]. Both analyses consistently identified stable high-QoL/low-distress and chronic low-QoL/high-distress groups, as well as trajectories reflecting recovery or delayed response. Differences primarily relate to class composition and size. Recovery and delayed-response classes in the present study were substantially smaller (approximately one-half to one-third) than those reported by Vangelis et al. In addition, we did not identify a delayed-deterioration QoL trajectory (characterized by initially moderate-to-high levels followed by decline) or an “unstable good” V-shaped trajectory attributed to transient treatment effects, although V-shaped patterns were also evident within other trajectories. In the present analysis, mean trajectory shapes cannot exhibit sharp V-shaped patterns, as they are approximated by smooth (linear or quadratic) functions.

Differences primarily reflect the tendency of shape-based clustering to distinguish and group in one class timing-related subpatterns, while the model-based approach assumes groups are defined by shared growth parameters (like specific intercepts and slopes) and accounts for within-class variability.

The kmlShape algorithm is designed to identify complex trajectory shapes [55] and preserve it in the mean trajectory. Such an example is a sharp “V-shaped” (drop-then-recovery) patterns that may better reflect short-term treatment effects. However, this approach clusters individuals based on shape similarity alone, regardless of the timing of changes (e.g., whether a decline occurs at month 3 or month 9), and is insensitive to baseline (intercept) differences, as it focuses on relative movements along the path. (e.g. it can recognize that a patient dropping from 80 to 50 and back to 80 has the same V-shape as one dropping from 60 to 30 and back to 60).

In contrast, LCGA/GMM approach is a statistical tool that accounts for uncertainty and individual-level differences within the groups. They are better suited when the timing of changes is clinically meaningful (e.g. does dropping early predict worse outcomes than dropping late?) and due to their parametric formulation, are sensitive to intercept differences (e.g. a “High Baseline” group and a “Low Baseline” group are statistically different populations, even if they share the same slope or recovery shape). Moreover, because GMM typically fits smooth (e.g., quadratic) trajectories, very sharp V-shaped patterns may not be identified as distinct classes and can instead be treated as outliers around a stable underlying trend. Notably, in our results, the excellent QoL trajectory showed no treatment-related decline based on visual inspection of individual-level trajectories (i.e. the methodology distinguished a class that is not affected by treatment), whereas declines at treatment time points were observed across all other classes.

Differences between our findings and those reported by Karademas et al. (2023) likely reflect important methodological distinctions rather than true inconsistencies. In the previous analysis, trajectories were collapsed into dichotomous outcomes, whereas our approach retained multiple, distinct trajectory classes, allowing clinical factors—such as neoadjuvant chemotherapy and triple-negative disease—to emerge as relevant for specific QoL patterns rather than for broad outcome categories. As a result, predictors associated with particular clinical courses may have been obscured in the earlier binary framework. In addition, the use of random forests in the prior study may have influenced variable importance estimates in the presence of highly correlated predictors, where importance can be distributed across correlated features or masked entirely. In contrast, the elastic net framework explicitly addresses multicollinearity, enabling more stable identification of clinical predictors whose effects are conditional on specific trajectory memberships.

There is evidence that suggests that GMM may outperform alternative approaches, including k-means, in certain settings [57,58,59,60]. However, further methodological work is needed to clarify the conditions under which different longitudinal clustering methods converge or diverge. Hybrid or comparative frameworks ([57]; latrend R package) may offer promising avenues for improving the identification of clinically meaningful longitudinal phenotypes.

Future research must prioritize the early identification of QoL and distress profiles immediately following a breast cancer diagnosis. This period serves as a critical window for proactive risk stratification and the development of personalized supportive care plans. While the majority of patients maintain favorable QoL trajectories, a clinically vulnerable minority experiences persistent or deteriorating functional outcomes that require intensive, targeted intervention.

To improve the precision of these models, larger, well-characterized cohorts are necessary to capture the heterogeneity and temporal dynamics inherent in these high-risk subgroups. Advanced machine learning (ML) approaches, such as XGBoost or Random Survival Forests, offer a unique opportunity to integrate multidimensional factors, including clinical biomarkers, and psychosocial variables. However, researchers must employ robust techniques like Synthetic Minority Over-sampling (SMOTE) [61] or engineered up-sampling (ENUS) [62] to handle significant data imbalance, ensuring that models do not overlook the high-risk minority.

Furthermore, the transition to clinical practice requires interpretable and calibrated prediction frameworks [63]. When paired with interpretable frameworks like SHAP for clinical transparency [64] or Conformal Prediction [65,66,67] for calibrated risk estimates accompanied by rigorous confidence intervals, ML models can support shared decision-making without the typical opacity of black-box systems, thus, fostering trust in Clinical Decision Support (CDS) tools. By identifying key variables associated with trajectory membership, this analysis lays the foundation for clinically meaningful risk models that support shared decision-making and the efficient allocation of survivorship resources.

4.1. Strengths and Limitations

The methodological approach has several key strengths. Stability analysis across multiple imputations and repeats ensures findings are reliable, not accidental results from a single data version or data partition. Elastic Net (α=0.5) manages correlated clinical variables effectively. This prevents the model from arbitrarily picking one and ignoring the other, unless one is clearly a redundant proxy. Forcing clinical sites into the model is a critical “stress test” that identifies truly universal predictors by removing clinical site/country-specific bias and institutional confounding. The use of the λ=1se rule results in a simpler, more interpretable model that is less likely to overfit and more likely to generalize to new patient populations.

Several limitations should be considered when interpreting the findings. Prediagnosis levels of psychological variables were not available, limiting our ability to distinguish cancer-related changes from pre-existing conditions. Not all variables were assessed at baseline and month 3, which may have constrained the temporal resolution of the analyses. Differences across clinical sites, including variations in care pathways and cultural contexts, may have contributed to heterogeneity in patient-reported outcomes despite adjustment for site effects. In addition, the small sample sizes in specific trajectory classes may have constrained the stability of the model, increasing the risk of overfitting and limiting the generalizability of predictors identified within these smaller cohorts. Finally, some analyses were cross-sectional in nature, which precludes causal inference.

5. Conclusions

Findings suggest that heterogeneous subgroups of patients follow distinct adjustment pathways over time, with differences explained by variations in symptom burden, functional scales, coping styles and affect. These results provide insights into the multidimensional determinants of resilience and highlight potential intervention targets to support individualized care and long-term quality of life and mental health in women with breast cancer.