Data source
Data were obtained from the China Health and Retirement Longitudinal Study (CHARLS), a prospective cohort study of adults aged 45 years and above and their partners in China. The study used a multistage probability sampling strategy to select participants. Participants were enrolled in 2011 and were followed up in 2013, 2015, and 2018. Detailed study design of the CHARLS is available elsewhere [25].
A total of 17,708 individuals completed the baseline survey. Of these, 15,186 were surveyed at wave 2 (2013), 13,565 at wave 3 (2015), and 11,988 at wave 4 (2018). All participants with one or more missing values in chronic conditions (N = 3,164) were excluded from the data analysis. In addition, all participants who were under the age of 45 years (N = 244) or had multimorbidity (N = 3,032) at baseline were also excluded. Thus, we included 5,548 participants in the current analysis. Further, in analyzing the association between multimorbidity trajectories and incident ADL disability, additional 576 participants were excluded for the following reasons: having missing values in the baseline ADL scale variable, or having ADL disability at baseline. Detailed procedure of sample selection is shown in Fig. 1. The biomedical ethics committee of Peking University approved the CHARLS and written informed consent was obtained from all participants.
Measures
Chronic conditions and multimorbidity
We ascertained the history of 14 chronic diseases by asking “Have you been diagnosed with the following conditions by a doctor”: hypertension, dyslipidemia, diabetes, cancer, chronic lung diseases, liver disease, heart disease, stroke, kidney disease, stomach disease, emotional problems, memory-related disease, arthritis, and asthma. These diseases are the leading causes of death in Chinese and are generally irreversible [26]. The self-reported data on chronic diseases are in high agreement with that based on medical records [27]. All diseases and conditions were defined as a binary variable (1 = present, 0 = absent). Multimorbidity was defined as having two or more of any aforementioned diseases. Participants who were qualified as having multimorbidity during follow-up would be defined as having a new onset of multimorbidity, with all participants free of multimorbidity at baseline.
Physical function assessment
We used the Chinese version of the Katz Activities of Daily Living (ADL) scale to measure disability [28]. This scale assessed perceived difficulties in the six ADLs, including dressing, bathing, eating, getting into and out of bed, toileting, and controlling urination and defecation. Following previous studies [29], participants were classified as having ADL disability if they reported having any degree of difficulties in performing at least one ADL. Otherwise, they were considered as having no ADL disability. Participants who had no disability at baseline but were qualified as having ADL disability at any wave of follow-up were considered as developing an incident ADL disability.
Covariates
We included the following baseline variables as covariates: age (45–59, 60 + years), sex (male, female), self-rated health (poor, fair, good) [30], education (no formal education, primary school, middle school or above), body mass index (BMI) (underweight: < 18.5 kg/m2, normal: 18.5–22.9 kg/m2, overweight: 23.0–27.4 kg/m2, obesity: ≥ 27.5 kg/m2), [31] occupation (agricultural, non-agricultural, unemployed/retired), ever smoking (no, yes), alcohol drinking (no, yes), sleep duration (< 7 h/day, 7–8 h/day, > 8 h/day), marital status (married/cohabitation, single: divorced, separated, widowed, or never married), annual household expenditure (≤ 2800 yuan, 2801–4846 yuan, 4847–8325 yuan, > 8325 yuan), social participation (whether the respondent participated in any social activities), location of residence (rural, urban) and public health insurance status (no, yes).
Statistical analyses
The multimorbidity trajectory groups were derived according to the longitudinal data of chronic conditions. But due to the high dimensionality of the longitudinal data of chronic conditions, we cannot directly derive the trajectory groups. To address this challenge, we first used exploratory factor analysis (EFA) to construct multimorbidity patterns to reduce the dimensionality of the 14 chronic conditions. We then calculated the corresponding factor scores of the multimorbidity patterns and used the factor scores as the outcome variables in deriving trajectory of the multimorbidity using the group-based multi-trajectory modeling (GBMTM).
We used EFA to determine the latent multimorbidity patterns underlying the 14 chronic conditions in the baseline sample. The patterns were determined based on their interpretability. We used the weighted least squares means and variance (WLSMV) estimator to estimate the EFA model [32]. We evaluate the goodness of fit by the comparative fit index (CFI), standardized root mean square residual (SRMR), root mean squared error of approximation (RMSEA), and tucker and lewis index (TLI) [33,34,35]. For better interpretations, an oblique rotation of factor loading matrices was performed, with each resulting factor loading representing the strength of association between each condition and the latent multimorbidity patterns. A factor loading of ≥ 0.40 indicates a strong association. The multimorbidity patterns were named according to the conditions that were most strongly associated with them. After identifying the latent multimorbidity patterns at baseline, we estimated the factor scores for each pattern over the follow-up using factor loadings of the multimorbidity patterns at baseline. The multimorbidity pattern score ranged between -0.31 and 2.68, with a higher score suggesting a greater number of conditions belonging to a specific multimorbidity pattern. We excluded conditions with a prevalence of < 1.0% in the baseline sample to achieve better robustness [36].
Based on the longitudinal factor scores of multimorbidity patterns, we identified subgroups of participants with similar joint trajectories of the multidimensional scores using the GBMTM. In this model, the outcome variables were the factor scores of each multimorbidity pattern over the follow-up, and these variables were modeled using the censored normal distribution by setting the lowest score (-0.31) as the censored minimum and the highest score (2.68) as the censored maximum. Based on previous research on multimorbidity trajectories in Korean older adults [37], we hypothesized that there would be 2–6 distinct trajectories of multimorbidity. Model fitting proceeded iteratively by comparing models with a varying number of groups (2–6 groups) and shapes of trajectories (linear, quadratic, and cubic). Model selection was based on the Bayesian information criterion (BIC) and Akaike’s information criterion (AIC) value, in which the model with the lowest BIC and AIC value was preferred. In addition, an ideal model is one in which the proportion assigned to each trajectory group (based on the maximum posterior probability rule) is greater than 5%; the average posterior probability of group membership is at least 0.7. Finally, the final models should have sufficient clinical relevance and interpretability [38].
We used the generalized estimating equation (GEE) model with an independent correlation matrix to evaluate the association between multimorbidity trajectories and incident disability. The independent variable was multimorbidity trajectory group, with participants without multimorbidity as the referent group. The dependent variable was disability status (yes, no) which was modeled with a logit link function. All variables included in this study were repeatedly measured in 2011, 2013, 2015, and 2018. Three models were fitted consecutively: model 1 was adjusted for time, model 2 was additionally adjusted for sociodemographic characteristics (age, gender, marital status, living area, occupation, annual household expenditure, social participation, and public health insurance status), and model 3 was further adjusted for health-related characteristics (self-rated health, BMI, smoking, alcohol drinking, and sleep duration). We reported the odds ratios (ORs) with 95% confidence intervals (CIs) of the association between multimorbidity trajectory groups and disability. We tested the interaction between multimorbidity trajectory groups and time in the GEE model, but the interaction was not statistically significant, therefore the interaction term was not included in the models.
The EFA analyses were performed with Mplus version 8.0, the GBMTM was conducted using the “Traj” program in Stata version 16.0 (StataCorp, College Station, TX) [39], and the GEE model was fitted using R 4.1.2. A two-sided P value < 0.05 was considered statistically significant.
Sensitivity analysis
We conducted a set of sensitivity analyses: First, we performed three EFA using the cross-sectional data from 2013, 2015, and 2018 waves of CHARLS to assess the robustness of the EFA, where we included the same chronic conditions as those at baseline to make the results comparable with those of the main analysis. Second, to evaluate the influence of incomplete data on the trajectory analysis, we conducted sensitivity analysis by restricted participants with complete data for at least 3 waves (N = 5,888). Third, to evaluate the effects of the excluded sample, the characteristics from samples with and without participants who were excluded during sample selection were compared and fitted the GEE model regarding the association of multimorbidity and incident disability among those excluded from the current analysis. Fourth, we conducted an analysis that only excluded participants with multimorbidity at baseline (N = 6,630) to assess the impact of the excluded sample on the study results. Finally, we divided participants without multimorbidity into no morbidity and single morbidity groups to examine the impact of single morbidity on the association between multimorbidity trajectories and incident disability using the GEE model.