The Seek and Treat for Optimal Prevention of HIV/AIDS (STOP HIV/AIDS) cohort is a de-identified population-based cohort of all diagnosed PLWH in BC, followed between April 1^st, 1996, and March 31^st, 2020. The STOP HIV/AIDS cohort was formed through the annual linkage between the BC Centre for Excellence in HIV/AIDS’ Drug Treatment Program registry and other provincial administrative databases, as described in Supplementary Table S1 [23]. It is important to mention that, in BC, ART and related medical and laboratory monitoring are available free of charge to all PLWH residing in the province, which is likely to minimise any bias related to access to care in our analyses.

Eligible participants met the following criteria: (i) initiated ART between January 1, 2001, and December 31, 2018, in BC; (ii) were ≥18 years old at the first ART date; (iii) had >1 pVL measurement between the first ART date and the end of follow-up; (iv) were ART naïve at baseline; (v) had a pVL ≥50 copies/mL at the first ART date; (vi) had >1 year of follow-up after ART initiation; (vii) had ≥5 years of administrative data before the randomly selected index date to ascertain comorbidity status [24]; and (viii) had >1 CD4 measurement between 1.5 years before and 6 months after the index date. The index date itself was chosen at least one year after the participant’s ART initiation (i.e., between one-year post-ART and the end of follow-up). Index dates were chosen randomly so that results from the final model would not be relative to a specific clinical visit date. We excluded PLWH who had documented structured treatment interruptions or participated in blinded trials during follow-up (see Supplementary Figure S1). These criteria were used to enhance the reproducibility of our study and minimise potential biases.

We followed all eligible participants for >1 year from the index date until i) December 31^st, 2019; ii) the last contact date (i.e., the last filled ART prescription refill date, the last available laboratory test date or the date of the last interaction with the healthcare system); or iii) the date of death (all-cause), whichever came first. Baseline pVL and CD4 nadir (lowest CD4) values were obtained within the above-mentioned window period. We omitted anyone missing a CD4 nadir measurement within the window period. Supplementary Figure S1 outlines the step-by-step process for inclusion in the final analytic sample.

The outcome of interest was a one-year all-cause mortality probability. The mortality data was obtained from the BC Vital Statistics Agency database, available through STOP HIV/AIDS [25]. We chose chronic comorbidities: CVD, CLD, COPD, CKD, diabetes mellitus (DM), hypertension (HTN), substance use disorders (SUD), Alzheimer’s/dementia (AD/D), NADC, personality disorder (PD), schizophrenia (SCZ), mood/anxiety disorder (MAD), rheumatoid arthritis (RA), metastatic solid tumor (MST), osteoporosis (OP), asthma, peptic ulcer (PU), and osteoarthritis (OA), which are responsible for the excess mortality observed among PLWH despite the success of ART [7, 10]. These comorbidities were assumed to be irreversible once identified and considered present until the end of follow-up. In line with recent recommendations [26], we selected a broad list of 18 comorbidities that are most prevalent and clinically impactful among people living with HIV. We used administrative data definitions for each condition (see Supplementary Table S2), ensuring that we captured a wide spectrum of chronic comorbidities. We acknowledge that other conditions may also be relevant; future work with more comprehensive data could explore additional comorbidities.

The prevalence of each comorbidity was determined using a five-year lookback window from the index date [24], which is consistent with recommendations from Nanditha et al. [24], who demonstrated that longer windows substantially improve the accuracy of identifying chronic diseases in administrative data. Although a shorter window (e.g., two years) might have retained more participants, it risks missing previously diagnosed, but currently stable, conditions. Therefore, we opted for five years to minimise misclassification of comorbidities, aligning with other validated comorbidity indices. We acknowledge that this requirement reduced the eligible sample, but we believe it provides a more robust identification of chronic comorbidities. We ascertained diagnoses of these comorbidities based on the BC Ministry of Health (BC-MoH) case-finding algorithms or literature (where applicable) using the International Classification of Diseases (specifically the ninth revision with clinical modifications, and the Canadian tenth revision version) and Drug/Product Identification Numbers (Supplementary Table S2) [27]. Due to limitations inherent in administrative datasets, we were unable to capture the severity levels or staging of these chronic comorbidities. Consequently, comorbidities in our models were treated as binary variables (present vs. absent).

Other predictors included sex at birth (male, female), age (continuous), years on ART (i.e., years since ART initiation; continuous), pVL suppression (i.e., all pVL measurements <200 copies/mL; yes/no), CD4 nadir (<50, 50-199, 200-349, ≥350 cells/mm³; or continuous), and income assistance (yes/no), measured within the window period mentioned above. In multivariable analysis, we considered linear, quadratic, and cubic polynomial forms for age per 10-year increase, years on ART, and CD4 nadir per 100 cells/mm³, and we centred these variables at the median to control for multicollinearity. Due to a high proportion of missing data for ethnicity in our database, we used ‘income assistance’ as a proxy for socioeconomic status.

Categorical variables were expressed as counts and percentages while continuous variables were expressed as median values with 25^th to 75^th percentiles (Q1-Q3). In bivariable analyses, we compared categorical variables using the Chi-Square or Fisher’s exact test under the applicable conditions and continuous variables using the Wilcoxon rank-sum test [28]. We did not have an independent cohort for external validation. Thus, we performed internal validation using data-splitting (or split-sample) and a bootstrapping-based method and presented both ways. This method accounts for the over-optimism of using unadjusted bootstrapping results and adjusts the model for overfitting [29]. Cox proportional hazard regression models were used to derive the MRPi to predict one-year all-cause mortality probability. The coefficients of the multivariable models were presented as adjusted hazard ratios (aHRs) with 95% Wald Confidence Intervals (CI). Analyses were performed using SAS version 9.4 (SAS Institute, Inc. Cary, NC, USA) and R statistical software version 3.0.2 (R Foundation for Statistical Computing, Vienna, Austria). All p-values were double-sided, and significance was at the 5% level.

To improve readability and reduce redundancy, we provide a concise overview of the model-building steps here, while detailed derivation procedures are given in Supplementary Text S1.

We randomly split the sample into training and test datasets with a ratio 2:1 [29]. The training dataset was used for the model development and later validated on the test dataset. We built Cox proportional hazard regression models on the training dataset to assess the relationship between selected predictor variables and one-year all-cause mortality. We used a variable selection approach based on the Akaike information criterion and Type III p-value [30]. We considered the possibility that certain predictors (e.g., age, CD4 nadir, viral load suppression) might interact with comorbidities to influence one-year all-cause mortality. Accordingly, we tested these interaction terms in preliminary models. However, none of these interactions reached statistical significance or meaningfully improved model discrimination (Harrell’s c-statistic) or calibration. To maintain parsimony and interpretability, we, therefore, excluded interaction terms from the final model. We adopted a three-model approach to systematically assess the incremental predictive gain of each additional set of covariates. Model 1 included only comorbidities, Model 2 added age and sex at birth, and Model 3 further incorporated HIV-specific variables (e.g., plasma viral load suppression) and socioeconomic status (income assistance). This stepwise strategy allowed us to evaluate how much discrimination and calibration improved with each new set of predictors.

We used the coefficients derived from the chosen final predictive model (from the data-splitting method) to create the MRPi. First, we multiplied each coefficient by the value of the corresponding predictor to get a risk coefficient. Then, we summed the risk coefficients for each participant. Next, participants’ MRPi scores were calculated as shown in Supplementary Text S1. Scores range from 0 to 100 (higher values = higher outcome risk). Next, the predicted one-year all-cause mortality probability was computed using the baseline Cox survival function and the sum of risk coefficients [29]. To compute each participant’s one-year all-cause mortality probability, we applied the baseline Cox survival function to the sum of each individual’s risk coefficients. The full equation, including the baseline hazard and parameter estimates, is provided in Supplementary Text S1. Also, to facilitate clinical use, we developed both a web-based calculator and an Excel-based tool to estimate the MRPi and the corresponding one-year all-cause mortality probability. The web-based calculator provides a user-friendly interface, allowing clinicians to enter key patient characteristics to obtain an immediate mortality risk estimate. We have included the link to the web-based calculator and instructions for using the Excel tool in the Supplementary Data File 2 (Mortality Risk Prediction Calculator).

We assessed the model’s performance through discrimination and calibration. For predictive discrimination, we calculated Harrell’s concordance statistic (c-statistic) [31]. The c-statistic ranges from 0.5 to 1.0, with 0.5 considered discrimination by chance alone, 0.70-0.79 considered good, and ≥0.80 excellent [32]. We calculated the 95% CI for each c-statistic, assuming a normal distribution [28]. The final predictive model was the model with the best discriminative ability. For calibration, we assessed the agreement between the predicted mortality probability and observed mortality probability using a calibration plot [29].

The Cox model coefficients estimated from the training dataset were applied to the test dataset along with the estimation of the c-statistic. Calibration was performed as follows. We first predicted one-year all-cause mortality using our final Cox model derived from the training dataset [20]. Second, we binned the participants in the test dataset into ten subgroups based on a 5-point interval of their predicted mean MRPi [20]. Next, we estimated the observed one-year all-cause mortality probabilities using the Kaplan-Meier (KM) method and 95% CI for each subgroup. Finally, we graphically compared the predicted and the observed one-year all-cause mortality probabilities for the subgroups. The agreement between the predicted curve and the observed points indicates how well the model was calibrated.

The bootstrapping-based internal method estimates the future performance of the model on new participants by providing bias-corrected estimates when the model is applied to a new sample [29]. We carried out this intuitive approach to emphasise the internal validity of the MRPi. To perform the bootstrapping-based internal validation, we first built a Cox model using the overall dataset, and the best model was chosen based on its performance, as described above. We estimated c-statistics and predicted survival probabilities by applying the model to the full cohort. Then, we drew 400 bootstrap samples with replacement from the full cohort; for each bootstrap sample, we fitted a Cox model to measure its apparent performance. The optimism was estimated following the method by Harrell et al. [29, 31]. We then used the average of the optimism over all 400 bootstrap samples as a correction factor to estimate a validated performance (optimism-adjusted) measure by subtracting the estimated optimism from the apparent performance. For the calibration plot, we calculated and compared the observed and adjusted KM (y-axis) to the predicted (x-axis) one-year survival probabilities, where participants were binned into ten subgroups using deciles of their predicted one-year survival probabilities [29]. A 45-degree line with a slope equal to one and intercept equal to zero (perfect calibration) indicates how well the model was calibrated. Any deviation above or below this line implies the difference between the observed and predicted on year survival probabilities.

First, we compared participants included in this study with those excluded to address external validity. Next, we performed sex-based analyses to assess whether the performance of the MRPi varied by sex at birth (male vs. female). Specifically, we repeated the final model-building steps described above, stratifying by sex.

Table 1 shows the descriptive characteristics of the full cohort, comparing the training and test datasets at the index date. Among the 4,387 participants, 3,570 (81%) were male, 3,116 (71%) had a suppressed pVL, and 762 (18%) had a CD4 nadir <200 cells/mm³. The median age was 47 years (Q1-Q3: 38-54), follow-up time was 3.3 years (Q1-Q3: 1.4–6.3), and time on ART was 4.0 years (Q1-Q3: 2.0–6.0). MAD was the most prevalent comorbidity (2,154, 49%), followed by SUD (1,720, 39%), CLD (735, 17%), asthma (596, 14%), HTN (561, 13%), and the least prevalent, RA (38, 1%). Participants in the test dataset and training sets were very comparable, except for asthma prevalence at the index date where those in the test dataset were significantly more likely to have asthma than those in the training dataset (Table 1).

Overall, 714 (16%) deaths were observed with a mortality rate of 38.18 per 1,000 PYs (95% CI: 35.46-41.06 per 1,000 PYs) by the end of follow-up, and 54.09 per 1,000 PYs (95% CI: 47.22-61.70 per 1,000 PYs) in the first year. The KM survival curve from the index date to the end of the study period is shown in Supplementary Figure S2.

Table 2 summarises the multivariable-adjusted association of the predictors with one-year all-cause mortality and the c-statistics for the training and test datasets. Model 3 demonstrated superior discrimination (0.8355, 95% CI: 0.8187–0.8523 in the training dataset and 0.7965, 95% CI: 0.7664–0.8266 in the test dataset) of one-year all-cause mortality to Models 1 and 2, and it was chosen as the final predictive model. Among all predictors in Model 3, SUD (aHR 2.60, 95% CI: 2.08–3.27), MST (aHR 2.27, 95% CI: 1.47–3.50) and NADC (aHR 2.16, 95% CI: 1.63–2.86) had the strongest association with one-year all-cause mortality. Increasing CD4 nadir per 100 cells/mm³ and years on ART were associated with a reduced aHR (aHR 0.84 per 100 cells/mm³ increase, 95% CI: 0.80-0.88 and aHR 0.94 per one-year increase, 95% CI: 0.89–0.98, respectively) for one-year all-cause mortality. Other predictors significantly associated with higher one-year all-cause mortality risk are shown in Table 2. The unadjusted results for each variable in the training dataset are shown in Supplementary Table S3. The adjusted model (Supplementary Table S4) yielded bootstrapping-based adjusted c-statistic approximately the same as the unadjusted c-statistic and slightly lower than the c-statistic in the final predictive model, suggesting good discrimination (0.8227 [95% CI: 0.8077–0.8377] versus 0.8188 [95% CI: 0.8180–0.8196]).

The median MRPi in our population was 49.6 (Q1-Q3: 43.1-57.5). As the MRPi increased, the predicted one-year all-cause mortality also increased. For example, let us consider the hypothetical case of a male participant, 46 years old, on ART for four years with no chronic comorbidities, a CD4 nadir of 500 cells/mm³, and a suppressed pVL. This participant’s estimated MRPi was 41 and his one-year predicted all-cause mortality probability was 0.8%. If he was also diagnosed with SUD, his MRPi would increase to 49, and his one-year predicted all-cause mortality would rise to 2.0%. Alternatively, if his CD4 nadir dropped to 150 cells/mm³, his MRPi would be 47, and his predicted mortality would be 1.6%. Finally, if his pVL became unsuppressed, but his CD4 nadir remained the same, his MRPi would be 45, and his predicted mortality would be 1.2%.

Figure 1 depicts the calibration plot comparing predicted versus observed one-year all-cause mortality probabilities across the ten subgroups of predicted mean MRPi in the test dataset for the data-splitting internal validation method. When we applied the final Cox model to the test dataset, we found that the predicted and observed mortality were closely aligned (visually) (Supplementary Table S5) for MRPi <75 over the one-year observation period, which includes most of the data (1,437, 98%). The model overpredicted the probability of mortality for MRPi greater than or equal to 75. However, the less-than-perfect alignment was in the tail of the MRPi distribution, representing those individuals who were more likely to die. Therefore, the overall calibration plot suggests the model was well calibrated with good agreement between the predicted and observed mortality probabilities. Supplementary Table S5 summarises the data in Figure 1.

Figure 2 (Supplementary Table S6) shows the calibration plot comparing the mean predicted and observed KM (unadjusted and adjusted) one-year survival probabilities using the full dataset (bootstrap-based internal validation). When the survival probability was <95%, there was a non-significant, modest difference between the predicted and observed one-year survival probabilities with overlapping 95% CI. Though there was a closer correspondence between the predicted and observed one-year survival as the one-year survival probability increased above 95%, some 95% CI did not cross the ideal line implying an underestimation of the predicted one-year survival probability in these subgroups. These inconsistencies were likely because the predicted survival probabilities in some subgroups were higher or lower depending on the sample size distribution. Note that we did not have observed survival probabilities <75% in this dataset. Overall, there was good agreement between the observed KM and predicted one-year survival probabilities suggesting a good calibration.

The excluded participants were significantly more likely to be younger, and have longer median years on ART and lower mortality (Supplementary Table S7). Results of the sex-based model-building are shown in Supplementary Figures S7–S10 and Tables S8–S11. Females in the final analytic sample were younger than males. About 572 (70%) females versus 2084 (58%) males were <50 years (Supplementary Figure S3). Sex-based analyses indicated that females tended to have a higher comorbidity burden and slightly higher MRPi scores than males. However, the overall predictive performance of the MRPi (c-statistic) remained high in both subgroups. Multimorbidity was more common among females (602, 74%), and males (958, 27%) were more likely to have no comorbidity at the index date (Supplementary Figure S4). The most prevalent comorbidities were SUD, MAD and CLD for both males and females (Supplementary Figure S5). In addition, females experienced two times the burden of SUD and CLD compared with males (66% versus 33%, and 30% versus 14%, respectively). Supplementary Figure S6 shows the distribution of the MRPi by sex at birth. A total of 473 (13%) males versus 36 (4%) females had an MRPi between 30 and 39, while 225 (28%) females versus 448 (13%) males had an MRPi between 60 and 69.