Objective: The primary objective of this retrospective observational study was to describe the evolution of CD4 cell counts over time, for patients on antiretroviral therapy (ART) in Mildmay Uganda. It was also aimed at determining whether the evolution depends on selected patient characteristics.

Methods: Since repeated measures were taken from the patients, correlation has to be taken into account when analyzing these data. As such, this study illustrates the application of Linear Mixed Models to describe the changes in CD4 counts over time, for patients on ART. Patients with baseline age 15 years were considered in this study.

Findings: Results show, that the evolution could be described by a cubic function of time. It was also noted that the effect of gender on the evolution depends on the patient's baseline CD4 category, also on the NNRTI drug that the patient uses.

Conclusion: The benefit of early treatment was shown in this study i.e. patients who started ART at higher baseline CD4 counts evolved higher than those who started at lower CD4 counts. The effect of gender on the evolution depended on the patient's baseline CD4 category, also on the NNRTI drug used.

Keywords: CD4 count, Longitudinal analysis, linear mixed models, anti retroviral therapy (ART), HIV/AIDS

ART, anti retroviral therapy; HIV, human immunodeficiency virus; AIDS, acquired immunodeficiency syndrome; MARPs, most at risk persons; NRTI, transcriptase Inhibitor; AZT, zidovudine; d4T, stavudine; TDF, tenofovir; NNRTI, non-nucleoside reverse transcriptase inhibitor; REML, restricted maximum likelihood; ML, maximum likelihood; EB, empirical bayes; SAS, statistical analysis system

The Human Immunodeficiency Virus (HIV) causes Acquired Immunodeficiency Syndrome (AIDS) by reducing a person's ability to fight infection. HIV attacks an immune cell called the CD4 cell which is responsible for the body's immune response to infectious agents.¹ As such, the number of CD4 cells per cubic millimeter of blood is widely used as an important biomarker for progression to AIDS when studying the efficacy of drugs to treat HIV-infected patients.² It is often measured repeatedly over follow-up periods in large-scale studies.

Since the beginning of the epidemic, almost 75 million people have been infected with the HIV virus, about 36 million people have died, 35.3 million [32.2-38.8 million] people were living with HIV at the end of 2012. An estimated 0.8% of adults aged 15-49 years worldwide are living with HIV, Sub-Saharan Africa remains most severely affected with 71% of the people living with HIV worldwide and nearly 1 in every 20 adults is living with HIV.³

2012 HIV/AIDS statistics showed that 7.2 percent of Uganda's population was living with HIV. This amounted to an estimated 1.4 million people, which included 190,000 children. An estimated 62,000 people died from AIDS in 2011 and 1.1 million children had been orphaned by Uganda's devastating epidemic. HIV prevalence has been rising since its lowest rate of 6.4 percent in 2006. New infections are diagnosed in 150,000 people a year, of whom 20,600 are children.⁴

Antiretroviral Therapy (ART) services have been available to HIV patients in Uganda for a while now. In May 2013, Uganda began the process of updating the ART guidelines to include recommendations from the 2013 WHO ART Guidelines which are aimed to dramatically increase access to ART and in turn provide new opportunities to save lives, improve clinical outcomes and reduce HIV incidence. The new guidelines recommend initiation of ART in all adults and adolescents with HIV at CD4 cut of <500 cells/mm³ regardless of clinical stage. The guidelines also recommend immediate ART initiation for patients with: HIV and active TB disease, HIV and HBV co-infection with evidence of severe chronic liver disease, HIV positive partner in a sero-discordant sexual relation-ship, Most at risk Persons (MARPs) in hot-spots (fisher folks, Commercial sex workers, long distance truck operators), Pregnant and Lactating Mothers and for all HIV positive children below 15 years of age.⁵ It is worth noting that these guidelines have been changing over time, for instance, the 2010 recommendations had a threshold of 350 CD4 cells/mm³.³ Prior to 2010, there was a threshold of 250 CD4 cells/mm³.

Few studies have been conducted about longitudinal analysis of CD4 cell count data for patients on ART especially in Sub-Saharan Africa. For instance, Adams and Luguterah¹ report that the pattern of growth in CD4 cell counts was not linear. Picat et al.,⁶ use nonlinear mixed-effects models, and conclude that higher long-term CD4 counts were predicted for children starting ART younger, and with higher CD4 counts. Reda et al.⁷ use mixed models regression and report a substantial increment in weight and CD4 lymphocyte count among the patients who were taking ART in eastern Ethiopia.

This study was aimed at describing the evolution of CD4 cell counts, for patients with baseline age 15 years, who initiated ART between January 2009 and December 2012 at Mildmay Uganda. It was also aimed at determining whether the evolution depends on selected patient characteristics. By considering changes over time, the longitudinal approach has the added advantage of observing changes more accurately, by increasing the power and validity of measuring the change in CD4 cell counts.

Study Population

Mildmay Uganda (MUg) opened in 1998 to provide palliative outpatient care for people living with HIV/AIDS, and to act as a teaching and training centre for HIV/AIDS health care personnel in Uganda. MUg is primarily a treatment centre for persons living with HIV/AIDS and their families. It currently offers family centered care and support to approximately 24,000 clients. MUg has one central clinic located in Lweza, 12km out-side Kampala, the capital of Uganda. MUg is also involved in a district health systems strengthening programme to build HIV/AIDS health care capacity in existing health centres in 18 local districts.⁸ Against this background, this study was conducted using routine data for patients with baseline age 15 years, who started ART between 2009 and 2012 at MUg, and were on first line regimen drugs. CD4 counts are normally taken for patients on ART on routine basis every 6 months. The criteria for starting ART in Uganda followed the WHO guidelines. First-line ART comprised a Nucleoside Reverse Transcriptase Inhibitor (NRTI) backbone i.e. Lamivudine (3TC) plus one of Zidovudine (AZT), Stavudine (d4T) or Tenofovir (TDF), plus a Non-nucleoside Reverse Transcriptase Inhibitor (NNRTI) drug i.e. either Nevirapine or Efavirenz.

The study was approved by the institutional review board at Mildmay Uganda. Informed consent was not obtained from participants since this was a retrospective, observational study and data were analyzed anonymously. It is worth mentioning however, that patients in care at Mildmay Uganda, consent and assent to use of their data for scientific research.

The Data

The response variable (CD4 count) was log transformed because it was heavily skewed. Other variables that were considered in the study included: the baseline CD4 category ( 200, 201-350, 351-500, >500), baseline age categories (15-24, 25-34, 35-49, 50) Gender of patient (0=Female, 1=Male), NRTI backbone used (d4t+3TC, TDF+3TC and AZT+3TC), NNRTI drug used (Nevirapine, Efavirenz), Did the patient change drug combination (Yes/No), Year of starting ART (2009, 2010, 2011, 2012).

Exploratory data analysis (EDA)

EDA was done in order to gain insight into the data prior to analysis. Descriptive statistics were generated and tabulated while histograms were constructed to show the distribution of number of measurements per patient. 2 tests were used to compare baseline characteristics (in form of categorical variables) between patients treated with either Efavirenz or Nevirapine.

Individual pro le plots were then constructed to gain some rough picture about how subjects evolve as well as to provide indications in terms of between and within subjects variability. Moreover, this exploration provided ideas about what random effects to start with.

Average pro le plots were constructed to describe the mean evolution of CD4 counts, overall and according to different subgroups. From such exploration, indications were obtained, about the functional form of the evolution and also whether the evolution depends on given covariates. Smoothing using the Loess method was applied because the measurements were not equally spaced across the different subjects.

The variance structure was also explored. This was aimed at getting insight about how the variance evolves over time. This gave indications in terms of how the variance was to be modeled so that valid inferences could be made.⁹

The Linear mixed effects model

The Mixed Model is one of the statistical methods used to analyze repeated measurements in continuous longitudinal data in an easy, valid and flexible manner.¹⁰ The Linear mixed model handles longitudinal data analysis in the continuous case. Moreover, since it assumes a linear regression for each cluster separately, it can be used for data with unequal number of measurements per cluster. Against this background a linear mixed-effects model was fitted. According to Verbeke and Molenberghs,⁹ the linear mixed-effects model is defined as,

where $Y_i$  is the $n_i$ dimensional response vector for the $i^{th}$ subject, 1≤ I ≤ N, N is the number of subjects, $X_i$ and $Z_i$ are $\left(n_i\times p\right)$ and $\left(n_i\times q\right)$ dimensional matrices of known covariates, is a p-dimensional vector containing the fixed effects, $b_i$ is a q-dimensional vector containing the random effects, $b_i$ is q-dimensional vector of residual components, and D is a covariance matrix of random effects.

Using the general guidelines for model building proposed by Verbeke and Molenberghs,⁹ the preliminary mean structure was first stated using a model with plausible combinations of covariates considering the cubic evolution of CD4 counts of patients, with random intercepts and random slopes. This was followed by the investigation of serial correlation, reduction of random effects, and reduction of mean structure consecutively.

Including serial correlation, if present, is far more important than correctly specifying it.⁹ Hence to check for the appropriate serial correlation, likelihood based tests were used. Restricted Maximum Likelihood (REML) was preferred to Maximum Likelihood (ML) testing here, because it reduces the well-known finite sample bias in the estimation of the covariance.¹¹

The covariance estimates for the random effects were assessed using mixtures of chi-square distributions due to the boundary problem. This problem arises because these tests involve variances which cannot be negative. Thus a likelihood ratio test of the null hypothesis that a variance is zero is testing a null hypothesis that is on\the boundary of the parameter space" for a variance. One consequence is that the usual null distribution for the likelihood ratio test is no longer valid; instead, the null distribution is a mixture of chi-squared distributions.¹¹

Finally, to discover the most parsimonious mean structure, F-tests and Likelihood Ratio Tests were employed under ML estimation.

Since the random effects ( $b_i$ 's) of the linear mixed effects model are stochastic, Bayesian methods were used to estimate them. The estimates, known as Empirical Bayes (EB) estimates, were presented in scatter plots to help identify outlying pro les.

Statistical software

Statistical Analysis System (SAS) version 9.3 was used for statistical analysis and graphics. Some graphics were produced in R version 2.15.3. For statistical tests, a significance level of 5% was used.

Exploratory data analysis

Descriptive statistics: Table 1 below shows the baseline characteristics of the population studied with Age 15 years (n=3241), at the point of ART initiation. Overall, there were more females (68%) than males (32%). Patients starting nevirapine were predominantly female, more likely to use zidovudine (than Tenofovir or Stavudine) in their NRTI backbone and had lower median CD4 counts (195 versus 209 cells/mm³). The few patients who used Stavudine in their NRTI backbone were not considered in the subsequent analysis, even so because that drug has been phased out of use.

Figure 1 below shows the distribution of the number of measurements per patient. The largest number of patients (575) had two measurements, followed by those with four measurements (492) and least of all was one patient with 12 measurements.

Figure 1 Frequency of number of measurements per patient.

Tables 2 & Table 3 below show the changes that happened within the first line regimen drugs used. It can be seen that 89 patients changed from nevirapine to efavirenz based ART (see Table 2) while 76 patients changed from zidovudine to Tenofovir use (see Table 3).

Exploratory plots

Individual Pro les: Figure 2 depicts pro les for all patients. What can be observed is that generally there is evidence of between subjects variability as well as within subject variability. The subjects have largely variable CD4 values at the start and also possibly different evolutions over time, this suggests that perhaps linear mixed models with random intercepts and slopes could be plausible starting points.

Figure 2 Individual pro les for all patients.

Mean structure: Figure 3 shows the average evolution of CD4 counts over time. We observe a sharp linear increase in the mean CD4 counts up to month 4, then a little dip at month 6 and from month 9 up to about month 40 it appears to flatten out, followed by a more gentle increase for the later time points. In this regard, we might have to consider a higher order polynomial for modeling the mean structure.

Figure 3 N=57; Epidemiological distribution of the pathological fractures, traumatic fractures, and nonunion.

In terms of gender (see Figure 9), females seem to have higher average CD4 counts than males at all time points. In terms of age categorization (see Figure 10), apart from the first 10 months in which the subjects with (35-49) years seemed to evolve differently from the others, at the other time points all age categories seemed to evolve in a similar manner. It is worth noting though, that the patients above 50 years had the lowest CD4 counts on average at most of the time points.

In terms of the NRTI backbone (see Figure 11), patients using AZT seemed to evolve in a similar way to those using TDF. The patients using d4T evolved in a different way, however, since these were few, they were not included in further analysis.

In terms of the NNRTI drug (see Figure 12), for the first 7 months, patients on nevirapine and those on efavirenz seemed to evolve in a similar way, between 7 to 20 months, those on efavirenz seemed to have higher CD4 counts while beyond 20 months those on nevirapine had higher counts.

In terms of Baseline CD4 categorization (see Figure 13), it is evident that the evolution differs depending on the baseline CD4 categorization. It is also seen that, patients who start in higher CD4 categories remain with higher CD4 counts at almost all time points.

From Figure 14, we observe that there are minor variations in the average pro les between the patients who changed treatment and those who did not change. Whether these variations are significant will have to be tested using the fitted models.

Variance Structure: Figure 15 shows the observed variance function, which is changing over time, indicating the need for other random effects in addition to the random intercepts.

Selection of preliminary mean and random effects structures

Following from the observations in the exploratory analysis, good models that best describe the observed average trends and also reflect the observed variance structures were sought. Verbeke and Molenberghs⁹ suggest that in case of highly unbalanced datasets with many covariates, it is necessary to use the most elaborate model one is prepared to consider for the mean structure, making it as extensive to answer the research question. In this case, all main effects and their two-way interactions were deemed sufficient to answer the research question. Several polynomials were then fitted on the same graphs with the observed average trends and the polynomials (functions of time) that seemed to best reflect the observed average trends were chosen. Figure 4 below shows the various fits. The decision on which polynomial to use was however taken after exploration of subject specific regression models elucidated in the next subsection.

Figure 4 Smoothed average evolution of CD4 counts by gender.

The choice of random effects structures was based on the observed variance functions as highlighted in the exploratory data part.

Exploration of subject-specific regression models

To further explore subject-specific evolutions, subject-specific regression models were fitted. Polynomial models with increasing powers were fitted to capture the varied trends of the individual evolutions. Table 4 shows that Rmeta² for the regression model with the linear effect of time was low at 0.4196, and it radically increased when the time covariate was raised to higher powers. Tests for model extension were conducted and p-values are also shown. Models with higher powers of time proved to be significant improvements from models with lower powers as shown by the significant p-values for the tests. However, with the aim of generating a model which is not only flexible but also parsimonious, the cubic polynomial model was chosen.

The Linear mixed effects model

A linear mixed effects model was fitted, using the repeated measurements of the CD4 counts of each person as a response. They were fitted against a cubic function of time, including all the 2-way interactions of all the covariates. On top of that, random effects for the intercept and linear slope were included to account for the variability between the different subjects. Extending the random effects structure to the quadratic and cubic slopes led to failure of convergence; because of this, robust inference was considered in the final model to circumvent problems of possible misspecification.

The need for reduction of the random effects structure was assessed using a mixture of chi-square distributions as mentioned in the methodology. Hence, the p-value under REML estimation for the comparison of Model 1 (with both random intercepts and slopes) versus Model 2 (with only random intercepts) (see Table 5) can then be calculated as:

From the above result it was concluded that the variance structure should not be simplified by dropping the random linear slope effect from the model.

Finally, the most parsimonious model in terms of the mean structure was assessed using the likelihood ratio test. The model was found to have a significant t based on the null model likelihood ratio test, implying that it is indeed necessary to model the covariance structure of the data. Summary of test results are shown in Table 6.

The final model with only significant effects was thus formulated as: $\begin{array}{l}\log \left(CD4coun{t}_{ij}\right)={\beta }_0+b_{0i}+{\beta }_{1}Gender+{\beta }_{2}cd4ca{t}_i+{\beta }_{3}chang{e}_i+{\beta }_{4}yea{r}_i\\ +{\beta }_{5}nnrt{i}_i+\left({\beta }_6+b_{1i}\right)tim{e}_{ij}+{\beta }_{7}tim{e}_{ij}{}^2+{\beta }_{8}tim{e}_{ij}{}^3\\ +{\beta }_{9}Gende{r}_{i^{*}}tim{e}_{{}_{ij}}+{\beta }_{10}cd4ca{t}_{i^{*}}tim{e}_{ij}+{\beta }_{11}nnrt{i}_{i^{*}}tim{e}_{ij}\\ +{\beta }_{12}nnrt{i}_{i^{*}}tim{e}_{i{j}^2}+{\beta }_{13}Gende{r}_{i^{*}}cd4ca{t}_i+{\beta }_{14}Gende{r}_{i^{*}}nnrt{i}_i\end{array}$

Where: $\log \left(CD4coun{t}_{ij}\right)$ is the natural log of the CD4 count of the ith patient at the jth time, $tim{e}_{ij}$ is the time (in months) at which the jth CD4 measurement is taken for the ith patient and ${\beta }_0,{\beta }_1,\dots \dots ,{\beta }_{14}$ are regression coefficients while $b_{0i}$ and $b_{1i}$ denote subject specific deviations from the common intercept and linear slope. The latter are assumed to follow a multivariate normal.

Table 7 below shows the covariance parameter estimates from the final model. It shows that the highest variability came from the random intercepts. It also shows that the variance of the random intercepts was higher than that of the random slopes, pointing to higher between patient variability than the within patient variability. The covariance of random effects is negative, implying that patients with high baseline CD4 counts tend to have smaller slopes than patients with low CD4 counts at baseline.

Results from the final model (see Table 8) show that there was a significant cubic evolution of the CD4 counts. The effect of gender on the evolution depends on the patient's baseline CD4 category (see Figure 5), also on the NNRTI drug that the patient uses (see Table 9).

Figure 5 Smoothed average evolution of CD4 counts by Age group.

Figure 6 Smoothed average evolution by NRTI backbone.

Table 8 also shows that patients who started ART between the years 2010 to 2012 had higher intercepts (CD4 starting values) than those who started ART in 2009. It also shows that patients who changed treatment drugs had lower intercepts (CD4 starting values) than those who didn't change.

Figure 8 in the appendix shows the observed and predicted pro les plotted against each other. It shows that our model fits relatively well to the observed data.Empirical Bayes estimates for random effects were also calculated in order to check whether there were outliers or not. The scatter plot of random intercepts and slopes (see Figure 7-15), showed that there were some potential outliers which may need to be investigated further.

Figure 7 Smoothed average evolution by NNRTI drug used.

Figure 8 Smoothed average evolution by Baseline CD4 category.

Figure 9 Smoothed average evolution by change of treatment.

Figure 10 Smoothed variance function.

Figure 11 Observed and fitted mean functions.

Figure 12 Predicted CD4 count levels showing interaction of gender and baseline CD4 category.

Figure 13 Predicted CD4 count levels showing interaction of gender and NNRTI drug.

Figure 14 Observed (dotted lines) and Predicted connected lines) pro les for females (left pannel) and males (right pannel) disaggregated by baseline CD4 categories.

Figure 15 Scatter plot of random intercepts and slopes.

None.

Authors declare that there are no conflicts of interests.

1. Adams M, Luguterah A. Longitudinal analysis of change in CD4+ cell counts of HIV-1 patients on antiretroviral therapy (ART) in the Builsa district hospital. Eouropean Scientific Journal.2013;9(33):1857‒7881.
2. Guo X, Carlin B. Separate and joint modelling of longitudinal and time-to-event data using standard computer packages. The American Statistician. 2004;58(1):16‒24.
3. World Health Organization; 2014.
4. HIV and Aids in Uganda. Avert. 2017.
5. Ario AR. The National antiretroviral treatment guidelines for Uganda 2013. Atic newsletter. 2014;11(1):1‒3.
6. Picat M Q, Lewis J, Musiime V, et al. Predicting Patterns of Long-Term CD4 Reconstitution in HIV-Infected Children Starting Antiretroviral Therapy in Sub-Saharan Africa: A Cohort-Based Mod-elling Study. PLoS Med. 2013;10(10):e1001542.
7. Reda A A, Biadgilign S, Deribew A, et al. Predictors of Change in CD4 Lymphocyte Count and Weight among HIV Infected Patients on Anti-Retroviral Treatment in Ethiopia: A Retrospective Longitudinal Study. PLoS ONE. 2013;8(4):e58595.
8. Funk A, Kanters S, Nansubuga M, et al. Cohort Pro le: The MUg Observational Cohort. Int J Epidemiol. 2012;41(16):1594-1594.
9. Verbeke G, Molenberghs G. Linear Mixed Models For Longitudinal Data. Springer-Verlag: New York, USA; 2000.
10. Burton P, Gurrin L, Sly P. Extending the simple linear regression model to account for correlated responses: an introduction to generalized estimating equations and multi-level mixed modelling. Stat Med. 1998;17(11):1261‒1291.
11. Fitzmaurice G M, Laird N M, Ware J H. Applied Longitudinal Analysis. John wiley & Sons; 2004.
12. Wools-Kaloustian K, Kimaiyo S, Diero L, et al. Viability and e ectiveness of large-scale HIV treatment initiatives in sub-Saharan Africa: experience from western Kenya. AIDS. 2006;20(1):41‒48.