SubSaharan Africa’s population is among the most HIV infected in the world. In 2014, there were 1.4 million (1.2 - 1.5 million) new infections and around 790,000 (690,000 - 990,000) people died of AIDS-related illnesses . These alarming statistics indicate a generalised HIV epidemic. On the one hand, HIV/AIDS infection is correlated with high risk behaviors such as occasional unprotected sex. On the other hand, the burden of the disease is correlated with poverty, weakness of health system, discrimination toward girls and women, low education level, malnutrition, migration, stigmatization of persons living with HIV/AIDS (PLWHA) and low volume of funding for the fight against HIV/AIDS  . In addition, PLWHA are sometimes less observant during the treatment period or co-infected by Tuberculosis and then may face failure of treatment  .
In Burkina Faso, the overall HIV prevalence was estimated to be 1.8% in 2003 . The antiretroviral therapy (ART) program started in 1999 with the financial support of international funding institutions and NGOs . The Global Fund to Fight against AIDS, Malaria and Tuberculosis was the leading funding program between 2003 and 2007 . Its intervention in Burkina Faso was the subject of an evaluation at the end of the five-year period 2003-2007. It aimed at measuring the intervention effects on health service access as well as HIV/AIDS related morbidity and mortality. Among other questions that deserve additional attention, there is that of the effectiveness of care through the analysis of the response to treatment of people who benefited from it. Viral load and CD4 cells count are relevant indicators that can be considered for such an investigation. CD4 cells count is low cost, simple to measure and a good predictor of the HIV dynamics during treatment . Since it is unclear how health conditions at the beginning of treatment and the demographic characteristics are correlated to the temporal variations of CD4 cells count, we propose to address this gap in this paper. The aim of this study is to model temporal variations of CD4 cell counts in a sample of ART patients with a longitudinal tree regression approach.
Data come from the Global fund to fight AIDS, Tuberculosis and Malaria Five year evaluation survey database. The sample included both male and female older than 15 that had initiated ART between 1st Jan 2003 and 31st Dec 2007. Data consisted of socio-demographic characteristics and semesterly CD4 count evaluations during the follow-up period. More details on patient sampling and data collection can be found in . Finally, we applied the following patient selection criteria: at least two CD4 counts available of which one is the baseline CD4 count; first line ART regimen, and baseline WHO clinical stage available.
We also applied two operations to complete cases. First, we created 3 new variables: body mass index (BMI) and the CD4 count evaluation clinical visit index (Time). Then we corrected the positive asymmetry of the distribution of CD4 count by applying a square root transformation.
Baseline Data Description
Our study included 3459 ART patients from 14 ART centres. Most centres (98%) were urban. For most patients (95%), the regimen consisted of two nucleoside reverse transcriptase inhibitors (NRTI) plus one non-nucleoside reverse transcriptase inhibitor (NNRTI). Table 1 shows that 72% were female, 70% were married, 95% were infected by HIV1, 81% initiated ART at WHO clinical stages between 3 and 4, and 83% had CD4 count lower than 200. The median age at ART initiation was 35 (30; 40).
Table 1. Patients’ characteristics.
Note: 1) NGO: non-governmental organization; CDT: Centre de Dépistage et de traitement de la Tuberculose; CTA: Centre de Traitement Ambulatoire; CDVA: Centre de Dépistage Volontaire et Anonyme. 2) We used Kruskall-Wallis and Mann-Whitney tests to compareCD4counts distribution.
3. Statistical Methodology
3.1. Statistical Model
The statistical model considered in this study for the data analysis belongs to the class of the varying-coefficients regression models as introduced by  after . In this setting the regression coefficients are allowed to vary with respect to some covariates called effect modifiers or moderators. When dealing with the framework of the linear model, the dependence of the regression coefficients of a predictor on effects modifiers should be understood as the expression of an interaction between X and the effects modifiers. Let’s consider a response variable Y conjointly with covariates (X, Z) where and . Let’s denote and respectively the spaces of the potential values of X and Z and with x in and z in ; the varying-coefficients model will take the following form: . For the present study the functions are modeled through a binary tree as follows:
where is rectangular partition of and is the indicator function. This choice offers is a flexible way to consider interactions between the covariates Z and the covariates X without a tight specification of such interactions at the model statement step. Moreover this framework allows the identification of subgroups of individuals with specific shape of CD4 variation over time. Finally the model that will be considered for the data analysis will be expressed as follows:
and is the indicator function. This is a special case of varying-coefficient regression model .
3.2. Model Fitting Methods
We chose a generalized linear mixed model (GLMM) conditional inference tree algorithm for tree construction . It helps to overcome some instance instability and bias problems that are common in recursive binary tree models fitting  . In brief, at each iteration, the algorithm combines the fixed-effect estimation by a model-based partitioning algorithm (MOB)  with a linear mixed model random effect prediction. MOB is based on parameter instability tests . Instability refers to a significant difference in coefficient estimates for two subsets of the dataset. The algorithm builds the tree after repetitions of the following steps:
GLMM tree algorithm stopping criterion is based on the linear mixed model log-likelihood denoted . One strategy to fit the model to data is outlined in the following algorithm :
We used R package glmertree  for model fitting.
We have looked for an optimal joint value of a minimal terminal node size (minsize) and a numeric significance level by using Bayesian Information Criterion (BIC). That model selection strategy avoids to choose arbitrary values for minsize and . In our application, the linear mixed model includes two fixed-effects (Treatment and clinical visit index) and one random effect (patient ID). The covariates used to defined the subgroups are Age, Gender, Marital status, baseline CD4 count, WHO clinical stage, the mode of entry into the active list, the history of opportunistic infection at ART initiation, the HIV type and Body Mass Index.
3.3. Stability Analysis and Model Diagnostics
Stability is an essential property of a fitted regression tree model. It ensures that the potential instability of the model is minimized and then the model use, as for prediction task. A model stability can be assessed by fitting the model on data obtained by bootstrap resampling of the training dataset. Variables and cutpoints that were not selected by the original tree may be selected by replicate samples. Metrics for stability assessment include the relative variable selection frequency, the mean frequency of the variable selections per tree and the frequency of each cutpoint over the trees :
• The relative variable selection frequency for a partitioning covariate , equals the total number of replicate trees that have selected at least once, divided by the total number of replicate trees.
• The mean frequency of the covariate selections per tree for is the total number of times is selected for partitioning by a replicate tree over the repetitions, divided by the total number of replicate trees.
• The relative frequency of a cutpoint equals the total number of replicate trees that have selected to split the variable , divided by the total number of replicate trees.
A covariate selection is stable if its average split count is close to its number of selections in the original tree and its frequency of selection is close to 100%. Graphical methods are used to analyze variable cutpoints’ variability. A histogram is used to illustrate the cutpoint variability when the partitioning variable is numerical. It is expected that the cutpoints selected in the original tree have the highest frequencies (one or more peaks in the histogram). For an ordered categorical variable, a barplot is used to show the frequency of all possible split points. For an unordered categorical variable, a specific plot is used to visualize the partitions’ variability over the replicates. The same color is used for categories that belong to the same node. The combination of categories that corresponds to a partition observed in the original tree is marked on the right side of the plot by a solid red line. In addition, two dashed lines enclose the area representing the partition. Number(s) on the right side of the area indicate the level(s) of the corresponding split(s) in the original tree. In conclusion, a split point is stable if it is selected by most replicate trees.
We have implemented a non parametric bootstrap method for generalized linear mixed models  . Mainly, the method consists of three steps. First, for all , we computed scalars and by centering and scaling the predicted random variable and the predicted random variable respectively. Consider and the empirical means of and . Denote
the empirical variances of and of respectively. The marginal residuals is also centered and scaled. we compute the predicted values by replacing by and by . Note that the two random variables follow a standard normal distribution. Thirdly, to obtain each bootstrap dataset, the response y is replaced by a bootstrap of the predicted response .
To study the plausibility of the model assumptions, we proceed in two steps. First, we computed the least confounded residuals  and assessed the normality assumption using quantile-quantile plot. Secondly, we evaluated the homoscedasticity assumption by plotting predicted values against standardized conditional residuals. For this purpose, we wrote a R code for parametric bootstrap that account for fitted mixed model as add-on for the package stable learner  to realize the stability analysis of linear mixed model based recursive binary tree.
4.1. Identified Difference in Temporal Variation of CD4 Cell Counts with Respect to Interactions between Covariates
Our model highlights seven subgroups of patients with different temporal variations of CD4 counts (Figure 1). The set of patients at an advanced disease stage (CD4 < 225) is the most splitted. It consists of at least four subgroups of chronological increase of CD4 count. Baseline CD4 count is selected for the first split. It is the most correlated with the temporal variation of CD4 count. Baseline CD4 count is a predictor of sustained virologic response . In the subgroup of patients with baseline CD4 count ≤ 73 cells/μl, Gender is the most correlated with the CD4 count variation. Female patients have a faster chronological increase of CD4 counts than male patients ( and respectively).
Figure 1. Conditional inference tree grown on the study dataset. CD4 count is the study outcome. Baseline CD4 count (in CD4), Gender (Genre) and Age were selected as splitting variables. P-value in each edge refers to the smallest p-value found in parameter instability tests.
Uptake and adherence to HIV services may explain this difference. In the subgroup of patient with baseline CD4 counts between 73 cells/μl and 129 cells/μl as well as in the subgroup of patients with baseline CD4 counts between 129 cells/μl and 225 cells/μl, Age is correlated with CD4 count temporal variation. In these subgroups, patients with aged 34 years or below have a slightly higher chronological CD4 count increase compared with older patients. The variable Age selection might be explained by what is called discordant immunological response, likely to happen among older patients . Conversely, the group of patients with baseline CD4 count ≥ 225 cells/μl is homogeneous according to partioning rules.
Table 2 shows a chronological increase of CD4 count on average in all subroups ( and ). This finding corroborates the fact that ART treatment aims to reduce the viral load and progressively to restore the immune system. The increase is faster for patients with an initially weakened immune system. On average the coefficient is small when baseline CD4 count is large. In addition, the ratio . This shows that individual random effect is important in the analysis of CD4 cell counts time variation.
Subgroups also differ by the difference in treatment response between patients treated with 2NRTI + 1NNRTI regimen and those treated with other ART regimens (2 NRTI + IDV, 2 NRTI + LPV/r, 2 NRTI + NFV). It is commonly admitted that there is no optimal treatment for all patients. Chronological CD4 count levels are higher for patients treated with 2NRTI + 1NNRTI than for those treated with other ART regimens group in the subgroup of patients with baseline CD4 count ≤ 129 and aged 34 and over as well as in the subgroups of patients with baseline CD4 counts between 129 cells/μl and 225 cells/μl (Figure 2). The difference in chronological CD4 count levels between ART regimens is highest
Table 2. Summary of fixed-effects estimation.
Figure 2. Chronological increase of CD4 count in the subgroups identified by the model. Time is in 6-months time scale. (a)-(g) Time (semester).
in patients aged 34 or younger and with baseline CD4 counts between 129 cells/μl and 225 cells/μl. In the remaining subgroups, we found no significant difference ( ).
4.2. Conditional Inference Tree Stability Assessment and Model Validation
We performed the stability analysis with 500 bootstrap samples. Relative frequencies of selecting baseline CD4 count, Gender, Age are all equal to 100% (Table 3 and Figure 3). Their mean frequency of selections per tree are respectively equal to their selection frequency in the initial tree (Table 3). In addition, Figure 4 shows that all the bootstrap trees have splitted baseline CD4 count on the same levels as in the initial tree. Age is splitted at age 34 on level 3 as in the original tree. Similarly, Gender is splitted on level 3. Thus, baseline CD4 count, Gender and Age can be considered as stable. Most of the bootstrap trees (97.6%) and the original tree are identical (Figure 3). Finally, variable EntryMod need not be retained. It is only selected by 2.4% of bootstrap trees. To sum up, the fitted tree is stable. Finally, Figure 5 and Figure 6 show that there is no evidence against normality and homoscedasticity assumptions respectively.
Figure 3. Frequencies of the different trees built over the repetitions. Dashed horizontal red lines mark the frequency of the original tree. It is enclosed by a solid vertical red line at the right of the plot.
Figure 4. Stability of the cutpoints selection for partioning. Dashed vertical lines mark the original tree cutpoints. The number above a dashed vertical red line indicates the level at which the split occurred in the tree. For example, Age is splitted twice on level 3. For the categorical variable Gender, the split occurred between Male and Female on the level 3 in the original tree.
Figure 5. Quantile-quantile plot of standardized least confounded residuals.
Figure 6. Homoscedasticity assumption evaluation plot.
Table 3. Variable selection overview.
5. Discussion and Concluding Remarks
In this study, we have analyzed data from a retrospective cohort of adults who started Antiretroviral Therapy at different centres between 2003 and 2007 in Burkina Faso. The results showed that the population of ART patients in Burkina Faso may be split into seven subgroups, according to the CD4 count variation over time. The number of the subgroups indicates that this population is heterogeneous as regards treatment response. This feature is particularly pronounced among advanced infected patients (CD4 < 225). The high proportion of patients (83%) in this condition may partly explain this finding. Gezie et al.  found Age as a predictor of CD4 counts over time in Ethiopia. Younger age was positively correlated to CD4 count increase during the course of treatment. But no cutpoint was provided to allow a comparison with our finding.
On average, we have observed an increase of CD4 counts in all subgroups. This finding confirms the benefit of ART treatment whatever the patients’ health conditions. However, CD4 counts do not recover a normal value when ART is started at a low baseline CD4 count  . In addition, we found that the increase is lower among patients with greater baseline CD4 counts. Garcia et al. reported that the lower increase was not observed when it was adjusted by viral load at last determination . This study reveals that female patients have a faster increase than male patients. This can be explained by representations of masculinity in Burkina Faso that could lead to late presentation and poor compliance with treatment  .
Another finding of this study is that there are differences in treatment response between patients treated with 2 NRTI + 1 NNRTI regimen and those treated with other ART regimens (2 NRTI + IDV, 2 NRTI + LPV/r, 2 NRTI + NFV) in the groups of patients with baseline CD4 between 73 and 225. The differences could be related to a better adherence to ART regimen in the subgroups. Protease inhibitors are expected to be more effective than 1 NNRTI. But tolerance to treatment may explain why they are less effective.
Abrogoua et al.  have also shown a heterogeneity in the population of patients under ART in Cote d’Ivoire. Their study identified four and fourteen CD4 counts variation subgroups in two samples of ARV-treated patients enrolled in a clinical trial between 2006 and 2007 and followed for two years in Cote d’Ivoire. The samples included 87 nevirapine treated patients and 164 efavirenz treated patients respectively. Age, sex, weight, Karnofsky’s score, haemoglobin, occurrence of an opportunistic infection were used for CD4 cell counts modeling. Unlike Abrogoua et al. , we studied 3459 HIV infected patients followed up to five years in routine clinical practice. In addition, we included the treatment regimen in the model.
A part of this work was carried out during a stay of the first author as a doctoral student at the Laboratoire de Mathématiques et de leurs Applications de Pau (LMAP), at the Université de Pau et des Pays de l’Adour (UPPA). He is grateful for the welcome given to him by the members of LMAP and for having benefited from its resources. Authors thank Dr. Stephen Centrella and Mrs Marie Henriette Somda for their kind assistance in the English writing.
Simon Tiendrébéogo: Analyzed data, wrote analysis tools and wrote the paper.
Séni Kouanda: Contributed to write the paper and provide guidance on the clinical interpretation of the findings.
Blaise Somé: Contributed to write the paper and reviewed the manuscript.
Simplice Dossou-gbété: Contributed to the data analysis and to the writing of both analysis tools and the paper.
This article is original and contains unpublished material. The corresponding author confirms that the other authors have read and approved the manuscript and no ethical issues involved.