Cox Proportional Hazards Modelling of Time-to-Event Outcomes Among HIV Patients on Antiretroviral Therapy at Specialist Hospital Yola, Adamawa State, Nigeria

The evaluation of HIV and AIDS has shifted from simple case counting to the modelling of complex and time-dependent clinical outcomes (L​u​v​a​n​d​a​ ​e​t​ ​a​l​.​,​ ​2​0​2​3). Globally, the expansion of highly active antiretroviral therapy (HAART) has driven public health systems toward longitudinal patient tracking (M​o​n​t​a​n​e​r​ ​e​t​ ​a​l​.​,​ ​2​0​1​4). While traditional epidemiological indicators such as prevalence provide a useful cross-sectional snapshot of the proportion of infected individuals at a given point in time, they do not capture temporal variations in disease progression, patient retention, or the occurrence of clinical events (W​o​r​l​d​ ​H​e​a​l​t​h​ ​O​r​g​a​n​i​z​a​t​i​o​n​,​ ​2​0​2​3). As a result, modern biostatistics increasingly relies on survival analysis techniques, particularly the Cox proportional hazards regression model, to examine the relationship between risk factors and the time until specific clinical events occur (D​z​i​n​z​a​ ​&​ ​N​g​w​i​r​a​,​ ​2​0​2​2).

In Adamawa State, surveillance records from the Adamawa State Agency for the Control of AIDS (ADSACA) indicate that the prevalence of HIV declined from 2.5% in 2014 to 1.1% (L​e​a​d​e​r​s​h​i​p​ ​N​e​w​s​,​ ​2​0​2​4). As of 2023, 36,137 people living with HIV were reported to be on antiretroviral medication in Adamawa State (S​o​l​a​c​e​B​a​s​e​,​ ​2​0​2​3).

Historically, most epidemiological research in Adamawa State relied on static methods such as log-linear models, frequency counts, and binary logistic regression. For example, A​k​i​n​r​e​f​o​n​ ​e​t​ ​a​l​.​ ​(​2​0​2​3​) applied log-linear models to examine associations among gender, age, viral load, and marital status among 3,779 HIV/AIDS patients attending Specialist Hospital Yola, between 2019 and 2020. While such methods are useful for identifying demographic associations with HIV status, they are limited when the research interest involves long-term clinical outcomes and time-to-event data (D​z​i​n​z​a​ ​&​ ​N​g​w​i​r​a​,​ ​2​0​2​2). Traditional regression approaches typically collapse time-dependent outcomes into simple binary categories, thereby ignoring the temporal dimension of patient care and the presence of censorship arising from treatment default, relocation, or study closure (D​z​i​n​z​a​ ​&​ ​N​g​w​i​r​a​,​ ​2​0​2​2).

As standard prevalence metrics and cross-sectional analyses do not account for time-to-event dynamics, public health administrators cannot accurately determine when patients face the highest risk of adverse outcomes. This methodological gap also limits the ability of clinicians to quantify how clinical markers, such as viral load, interact with demographic characteristics over time to influence patient trajectories (D​z​i​n​z​a​ ​&​ ​N​g​w​i​r​a​,​ ​2​0​2​2).

In a neighbouring North-Eastern state, A​h​m​a​d​u​ ​e​t​ ​a​l​.​ ​(​2​0​2​4​) applied the Cox proportional hazards model to HIV patient records from Federal Medical Centre Jalingo, Taraba State. Their analysis identified CD4 count stage and marital status as significant predictors of the hazard of adverse outcomes, with patients in advanced clinical stages (stage IV) having approximately four times the risk of the event compared with those in earlier stages, and single patients having twice the risk of married patients. These findings from a contiguous state reinforce the value of survival modelling for understanding time-to-event dynamics among people living with HIV in the region and provide a useful comparative context for the present study at Specialist Hospital Yola.

A​p​e​a​g​e​e​ ​e​t​ ​a​l​.​ ​(​2​0​2​0​) applied the Cox proportional hazards model to time-to-event data of HIV/AIDS patients in Nigeria. Their study identified sex, age, and patients’ health status as significant predictors of mortality, while higher CD4 cell counts were associated with a reduced risk of death. The findings provide further methodological support for the use of the Cox model in analysing survival outcomes among HIV patients in Nigerian treatment settings.

In the local setting, A​k​i​n​r​e​f​o​n​ ​e​t​ ​a​l​.​ ​(​2​0​2​3​) used log-linear models to analyze HIV/AIDS data from 3,779 patients at Specialist Hospital Yola. Although the study successfully identified associations involving gender, age, viral load, and marital status, its reliance on static categorical analysis underscored the need for time-to-event methods tailored to the same facility and population.

I​k​p​e​ ​e​t​ ​a​l​.​ ​(​2​0​2​5​) conducted a retrospective cohort analysis of the Nigeria national HIV data repository to identify predictors of time to first treatment interruption among patients on antiretroviral therapy (ART) across 2,226 public facilities. Using Kaplan–Meier survival estimates, log-rank tests, and multivariable Cox proportional hazards regression, the study found that age grouping, ART regimen class, recorded viral load, CD4 cell count, World Health Organization (WHO) clinical stage, and facility location were predictive of treatment interruption. These findings underscored the value of national-level survival modelling for ART programme evaluation in Nigeria and provided a useful benchmark for facility-specific investigations such as the present study.

D​z​i​n​z​a​ ​&​ ​N​g​w​i​r​a​ ​(​2​0​2​2​) compared parametric and Cox regression models using HIV survival data from Malawi. Their findings demonstrated the advantage of the Cox model in handling censored observations in resource-limited African settings, providing methodological support for its application in similar contexts such as Adamawa State.

Beyond West and Southern Africa, C​h​i​r​n​e​t​ ​e​t​ ​a​l​.​ ​(​2​0​2​4​) examined time to viral load suppression among HIV patients receiving ART in the Gebi Resu zone, Afar Region, Ethiopia. Applying Cox proportional hazards regression, the authors confirmed that the proportional hazards assumption was satisfied and identified baseline WHO clinical stage, cotrimoxazole prophylaxis therapy, and drug adherence as significant predictors of delayed viral suppression. This study reinforced the broader methodological validity of Cox-based survival modelling for monitoring ART outcomes in resource-limited settings across diverse African contexts.

Collectively, these studies showed that while static models remained common in parts of Nigeria, survival analysis techniques, particularly the Cox proportional hazards model, offer superior insight into temporal risk factors. However, there remains a notable scarcity of Cox-based survival studies specific to Adamawa State. The present study addressed this gap by applying a multivariable Cox proportional hazards regression model to examine the influence of socio-demographic and clinical factors on time to treatment interruption or death among HIV patients on antiretroviral therapy at Specialist Hospital Yola, Adamawa State, Nigeria.

This study was conducted using HIV/AIDS patient data obtained from Specialist Hospital Yola, Adamawa State, Nigeria. The hospital is a major referral health facility in the state and provides comprehensive healthcare services, including HIV/AIDS diagnosis, treatment, monitoring, and ART management. It serves patients from different parts of Adamawa State and the neighbouring regions, rendering these areas a core centre for HIV/AIDS care and clinical research.

The data used in this study were obtained from the medical records of HIV patients on antiretroviral therapy at Specialist Hospital Yola, Adamawa State. These records formed the basis of an earlier analysis using log-linear models that examined associations between socio-demographic factors and HIV/AIDS prevalence (A​k​i​n​r​e​f​o​n​ ​e​t​ ​a​l​.​,​ ​2​0​2​3). The present study re-analyzed the same facility-based dataset with the addition of follow-up time and event status coding for the application of Cox proportional hazards models.

A retrospective survival analysis design was adopted using time-to-event data obtained from HIV/AIDS patient records. Survival analysis was considered appropriate because the study focused on the time until the occurrence of treatment interruption or death among patients on ART. The Cox proportional hazards regression model was applied to evaluate the effects of demographic and clinical factors on the survival outcome.

The analysis included 2,887 HIV patients on ART who had documented follow-up time and outcome information. Patients with incomplete data on the key variables required for the survival analysis (age, sex, marital status, and survival time/event status) were excluded. Two participants had missing or invalid viral load values; they were retained in analyses that did not involve viral load and were excluded only from models and tables that included viral load status. Consequently, the Cox regression models that included viral load were based on 2,885 participants.

The dependent variable was survival time, defined as the duration in months from the date of ART initiation to the date of the recorded Current ART Status (or last clinical contact). The primary outcome was a composite event of interruption in treatment (IIT) or death. Current ART Status was coded as the failure indicators as follows: Event (failure = 1) = IIT or Dead; Censored (0) = Active, Active-Restart, Active-Transfer In, or Transferred Out.

The independent variables examined were age group (categorized as ≤29 years, 30–39 years, 40–49 years, and ≥50 years), sex (female or male), marital status (married, single, widowed, divorced, or separated), and viral load status. Viral load was classified as suppressed (<1000 copies/ml) or unsuppressed (≥1000 copies/ml) according to the World Health Organization threshold.

The data used in this study were obtained from the medical records of HIV patients receiving antiretroviral therapy at Specialist Hospital Yola, Adamawa State. These records formed the basis of an earlier analysis using log-linear models that examined associations between socio-demographic factors and HIV/AIDS prevalence (A​k​i​n​r​e​f​o​n​ ​e​t​ ​a​l​.​,​ ​2​0​2​3). Permission to access the anonymized patient records was granted by the management of Specialist Hospital Yola for academic research purposes.

The Cox Proportional Hazards Regression Model was used to assess the effects of the independent variables on the hazard of treatment interruption or death. The model proposed by C​o​x​ ​(​1​9​7​2​) is given as:

where, h(t/X) is the hazard function at time t given the covariate vector X, h₀(t) is the baseline hazard function, and β is the vector of regression coefficients. In the present analysis, the covariate vector X contained indicator variables for age group (three indicators), sex (one indicator), viral load status (one indicator), and marital status (four indicators), giving a total of nine coefficients that were estimated.

The hazard ratio (HR) was used to interpret the effect of explanatory variables on survival outcomes and is defined as:

where, HR > 1 indicates an increased hazard (risk); HR < 1 indicates a reduced hazard (protective effect); HR = 1 indicates no effect.

The regression parameters of the Cox model were estimated using the partial likelihood estimation method. The partial likelihood function is expressed as:

where, the coefficients were obtained by maximizing this partial likelihood function.

Data cleaning, coding, and analysis were carried out using Stata statistical software. Descriptive statistics were first used to summarize the baseline characteristics of the study participants. Categorical variables were presented as frequencies and percentages, while continuous variables were summarized using means, standard deviations, and medians. Differences between participants who experienced treatment interruption or death and those who did not were examined using the chi-square test.

Survival probabilities were estimated using the Kaplan-Meier method, and differences in survival experience across groups were compared using the Log-rank test. Univariate Cox proportional hazards regression was performed to obtain crude HRs and 95% confidence intervals. All variables were subsequently entered into a multivariate Cox model to estimate adjusted HRs while controlling for potential confounding.

The proportional hazards assumption was assessed both graphically and statistically. Log-minus-log plots and Schoenfeld residual plots were explored for the covariates included in the final model. In addition, the global Schoenfeld residual test was conducted. A non-significant result p > 0.05 was interpreted as evidence that the proportional hazards assumption was reasonably satisfied. A p-value of less than 0.05 was considered statistically significant.

Figure 1 presents the Schoenfeld residual plots for the covariates in the final multivariate model. The plots showed no strong systematic trend over time. The smoothed residual lines remained roughly horizontal and close to zero. These graphical findings, together with the non-significant global Schoenfeld test (p = 0.669), confirmed that the proportional hazards assumption was reasonably satisfied.

Figure 2 shows the log-minus-log plots for viral load status, sex, age group, and marital status. For viral load status, the curves were approximately parallel, hence supporting the proportional hazards assumption. Plots for sex and age group also showed roughly parallel lines. Although some deviation was noted in the smaller marital status categories, the overall pattern did not indicate major violation of the proportional hazard assumption. These graphical findings are consistent with the formal Schoenfeld residual test.

The Kaplan-Meier survival curve (Figure 3) stratified by viral load status revealed a clear separation between the two groups. Participants with unsuppressed viral load (≥1000 copies/ml) had consistently lower survival probability throughout the follow-up period than those with suppressed viral load (<1000 copies/ml). This difference was highly statistically significant (Log-rank test, p < 0.001), indicating a strong association between unsuppressed viral load and increased risk of treatment interruption or death.

In contrast, the Kaplan-Meier curves by sex showed almost complete overlap between males and females across the entire observation period. The Log-rank test confirmed that there was no significant difference in survival experience between the two sexes (p = 0.609).

Similarly, when survival was examined according to age group (≤29, 30–39, 40–49, and ≥50 years), the curves remained close to one another with no distinct pattern of separation. The difference across age groups was not statistically significant (Log-rank test, p = 0.635), suggesting that age did not substantially influence time to treatment interruption or death in this cohort.

The Kaplan-Meier curves by marital status displayed minor variations, particularly among the smaller categories (widowed, divorced, and separated). Although the overall Log-rank test reached statistical significance (p = 0.028), the visual separation between the major groups (married and single) was limited, indicating that the practical difference in survival by marital status was modest. Note that the chi-square tests in Table 1 assess the cross-sectional association between each covariate and the occurrence of treatment interruption or death, whereas the Log-rank tests evaluate differences across the entire survival curves.

From the baseline characteristics (Table 1), a total of 2,887 participants were included in the analysis. The majority were female (70.1%) and married (58.4%). The largest age groups were 30–39 years (33.9%) and 40–49 years (32.4%). Most participants had suppressed viral load (<1000 copies/ml) (95%)while only 5% had unsuppressed viral load (≥1000 copies/ml).

There was no statistically significant difference in the cross-sectional distribution of sex (Pearson χ² p ≈ 0.54) Unsuppressed viral load was significantly more common among participants who experienced the event (Pearson χ² p ≈ 0.032). The corresponding Log-rank tests that evaluate differences across the entire survival curves are reported in Section 3.2 (sex p = 0.609; age group p = 0.635; viral load p < 0.001; marital status p = 0.028). The median follow-up time was 36 months.

In the univariate Cox regression (Table 2), each variable was examined individually. Age group, sex, and marital status showed no statistically significant association with treatment interruption or death. Participants aged 30–39, 40–49 years, and ≥50 years had crude HRs close to 1 when compared with those aged ≤29, and none of these associations reached significance. Similarly, male participants had a slightly elevated but non-significant hazard compared with females (Crude HR = 1.05; 95% confidence interval (CI): 0.88–1.24; p = 0.611).

In contrast, viral load status emerged as a strong predictor. Participants with unsuppressed viral load (≥1000 copies/ml) had nearly twice the hazard of treatment interruption or death compared with those who had suppressed viral load (Crude HR = 1.98; 95% CI: 1.44–2.73; p < 0.001).

All variables were subsequently entered into a multivariate Cox model to obtain adjusted estimates (Table 2). After controlling for age group, sex, and marital status, the association between unsuppressed viral load and treatment interruption or death remained strong and statistically significant (adjusted HR = 1.96; 95% CI: 1.42–2.70; p < 0.001). The HRs for age group, sex, and marital status changed only slightly and remained non-significant.

The detailed parameters of the final multivariate Cox model are presented in Table 3.

Based on these parameters, the model can be expressed as:

where, the coefficient for unsuppressed viral load (β = 0.6726) was the only statistically significant parameter (p < 0.001). Exponentiating this coefficient yielded an adjusted HR of 1.96, confirming that participants with unsuppressed viral load had a 96% higher hazard of the event compared with those with suppressed viral load, after adjustment for all other variables in the model.

This study examined factors associated with treatment interruption or death among HIV patients receiving care on ART at Specialist Hospital Yola, Adamawa State. The most important finding was that unsuppressed viral load (≥1000 copies/ml) was a strong and independent predictor of treatment interruption or death. After adjusting for age group, sex, and marital status, HIV patients on ART with unsuppressed viral load had nearly twice the hazard of treatment interruption or death compared with those who had suppressed viral load (adjusted HR = 1.96; 95% CI: 1.42–2.70; p < 0.001).

This finding is consistent with evidence from large Nigerian cohorts of HIV patients on ART. C​h​u​n​ ​e​t​ ​a​l​.​ ​(​2​0​2​2​), in a retrospective longitudinal study involving over HIV patients on ART across 18 Nigerian states, demonstrated that detectable viraemia significantly increased the risk of subsequent virological non-suppression and failure. Similarly, T​o​m​e​s​c​u​ ​e​t​ ​a​l​.​ ​(​2​0​2​3​), in a large cross-sectional analysis of 585,632 HIV patients on ART in Nigeria, reported that viral non-suppression (≥1000 copies/ml) remained an important programmatic challenge, with several clinical and demographic factors influencing its occurrence.

In the present study, socio-demographic variables (age group, sex, and marital status) were not significantly associated with the outcome in either univariate or multivariate models. This contrasts findings by T​o​m​e​s​c​u​ ​e​t​ ​a​l​.​ ​(​2​0​2​3​), who reported higher likelihood of viral non-suppression among males and younger age groups in Nigeria. A study conducted at the same facility (Specialist Hospital Yola) by S​t​e​p​h​e​n​ ​e​t​ ​a​l​.​ ​(​2​0​2​4​) found an adherence rate of 88.4% among HIV patients on ART, with forgetfulness being the most common reason for missed doses, thus highlighting the importance of adherence support in this setting.

The Kaplan-Meier curves in Figure 3 showed a clear separation between participants with suppressed and unsuppressed viral load. Participants with unsuppressed viral load (≥1000 copies/ml) had significantly poorer survival experience compared with those with suppressed viral load (Log-rank test p < 0.001). Although the Log-rank test for marital status reached statistical significance, the practical differences between the major categories were limited, to be consistent with the non-significant adjusted HRs for marital status.

The proportional hazards assumption was satisfied, as confirmed by Schoenfeld residual plots (Figure 1), log-minus-log plots (Figure 2), and the global Schoenfeld test (p = 0.669). This supports the validity of the Cox model used in the analysis.

This study has several limitations. It was a single-centre retrospective analysis based on routine hospital records, so residual confounding could not be eradicated. Important clinical covariates such as CD4 cell count, ART regimen, adherence measures, nutritional status, and co-infections were not available or were incompletely recorded, thus not being able to be included in the models. In addition, the quality and completeness of routine programme data may have introduced some misclassification of event or censoring status. These limitations should be considered when interpreting the findings. Future prospective multi-centre studies with richer clinical data are recommended.

This study among HIV patients on antiretroviral therapy at Specialist Hospital Yola demonstrated that unsuppressed viral load (≥1000 copies/ml) was the only independent predictor of treatment interruption or death. Patients with unsuppressed viral load had nearly twice the hazard of experiencing treatment interruption or death compared with those who had suppressed viral load, after adjusting for age group, sex, and marital status. Socio-demographic factors such as age, sex, and marital status were not significantly associated with the outcome.

These findings highlighted the critical importance of achieving and sustaining viral suppression among HIV patients on ART. Strengthening routine viral load monitoring, adherence counselling, and prompt clinical management of patients with unsuppressed viral load should be prioritized at Specialist Hospital Yola and similar treatment centres to improve patient treatment outcomes and to attain the broader goal of HIV epidemic control in Nigeria.