TABLE 1 | Distribution of explanatory variables for the 56,308 subjects aged between 17 and 65 years with marital status at entry to the study.

lies in its flexibility when modeling time-varying covariates We consider subjects between ages 17 and 65. Most cases in the African context for demography usually consider age for marriage to start at 15 years (Manda and Meyer, 2005). However, for this study, since there were no events occurring before 17 years of age, we opted to use a starting age of 17 years. It is then subdivided into q = 24 unit intervals between 17 and 65 years (which do not necessarily need to be equal). A subject's marital status is assumed to change not more than once within a period of 2 years, hence we chose to use 2-year time intervals in the analysis. Thus, we have intervals [17, 19], [19, 21], . . .[64, 65] which we denote as t = 1, 2, . . . , 24.

For assessing factors impacting the state occupied by a subject who possesses a number of characteristics at a particular age interval, a multinomial regression model is used (Grilli and Rampichini, 2007). For single transitions between marital states, into and out of a marriage, separate binary logistic regression models were utilized (Allison, 1982; Manda and Meyer, 2005; Steele, 2011; Clark et al., 2013). Separate response variables were created for the separate transition types, namely; marital dissolution (i.e., a dummy variable of 0 if subject is married and 1 if widowed or separated) or re-marriage (i.e., 0 if still widowed or separated and 1 if re-married). For analyzing ways of dissolving a marriage, a competing risks model (Allison, 1982; Jenkins, 1995; Steele, 2011) was used. We considered transitions out of a Married state, where Separated and Widowed are the competing events. A random variable, Yit, is created such that it is coded 0 whenever subject i is at risk of leaving a Married state, 1 when a transition is made into a Separated state and 2 when a transition is made into a Widowed state (the competing event).

Additionally, since observations are done on each individual repeatedly, there is possible correlation among observations TABLE 2 | Number of transitions between marital states from 1 January 2004 to 31 December 2016 for the 56,308 subjects aged between 17 and 65 years.

within the same subject. There may exist other unobserved subject-specific factors which may influence marital state occupation or transition between marital states. For example, one person may have reasons for not remarrying after the death of a partner, while another may easily remarry shortly. It is, therefore, important to account for these variations by including a subject-specific random effect in the analysis. In the binary logistic regression models, there is a separate random effect for each type of transition and these are allowed to be normally distributed with a mean of 0 and a variance of σ 2 . Moreover, since the marital states are possibly recurrent in a subject, the competing risks regression model controls for this using a multilevel model. The correlations will enable us to understand how one subject who experienced one type of marital transition is likely to experience the other. Tomita et al. (2017) used the same data set on HIV, accounting for its multilevel nature. However, they used continuous-time survival methods. Using discrete- time metrics, we fitted these in STATA15.1 using its inbuilt standard tools **xtlogit** for random effects logit model and **gsem** for multinomial logit.

## 3. RESULTS

## 3.1. Descriptive Statistics

The explanatory variables considered in this study include gender, if income is earned, if subject is employed, highest education level attained by subject, age at first sexual intercourse and age a first marriage for subjects who have ever been married.  **Table 1** shows the distribution of some variables of interest at entry into the study where applicable. In all categories of the variables, never married subjects had the highest proportion. At entry into the study, 13.6% of the female participants were married and 10.3% of the males were married. Most participants did not earn an income and of those who did, 24.4% were married and 69.5% were never married. Of those who did not earn an income, 5.9% were widowed.

Number of transitions (in person-years) into the different marital states are represented in **Table 2**, where most subjects remained in the Never Married state and the most common transition type was Married to Widowed with 1, 561 transitions followed by the Never Married to Married transition with 1, 274

transitions. During the follow-up period, most subjects remained in the Never Married state with 267, 284 person-years. The median age at first marriage was found to be 34 years. Of the 17, 124 who ever married, 29%, 38%, 28%, and 5% had their first marriage at ages ≤ 22, 23 − 30, 31 − 40, and ≥ 41 years groups, respectively.

**Figure 1** displays the proportion at each age, of subjects who were occupying each state. It is clear that below 46 years of age, the biggest part of the population was constituted by subjects who were never married. Moreover, as shown in **Figure 2**, the hazards of entering a first marriage were very low (close to 0) in the younger ages, but increased just slightly with age. **Figure 3** also displays the trends over time for marital state occupation. Over the whole study period, the proportion of never married subjects was the highest. It decreased until 2008 and then stabilized afterwards. Proportion of married subjects was considerably low but slightly increased with time, then declined after 2009. The rates of marital separation were almost constant over time while those of widowhood started to increase slightly after 2008.

## 3.2. Results for Marital State Occupation

We began by considering factors leading to staying in a marital state as shown in **Table 3** below. The odds ratios are also displayed in **Figure 4**. All confidence intervals reported were at 95%. Compared to never married subjects, we assessed the odds of being married. Males had a significantly lower likelihood of being married than females [Odds Ratio (OR) for gender was 0.84(CI = 0.76; 0.93)]. Similarly, the odds ratio for primary, high school and tertiary education are 1.68(CI = 1.43; 1.98), 2.42(CI = 2.04; 2.86), and 6.43(CI = 4.89; 8.47), respectively, thus the likelihood of being married relative to never married were highest among subjects with tertiary education compared to those who never went to school. Those who did not earn income and those who were not employed had a lower likelihood of being in a Married state relative to Never Married [OR= 0.93(CI = 0.87; 1.01) and 0.93(CI = 0.86; 0.99), respectively]. Between the years 2006 and 2007, a subject was more likely to be married than never married when compared

to years 2004 − 2005. After 2007, subjects were less likely to be married than never married and the trend reduced over the years up to 2016 [with OR = 0.86, CI = (0.78; 0.94), OR = 0.71, CI = (0.64; 0.78), OR = 0.60, CI = (0.54; 0.67), OR = 0.50, CI = (0.44; 0.56), and OR = 0.43, CI = (0.38; 0.48)] for periods 2008 − 2009, 2010 − 2011, 2012 − 2013, 2014 − 2015 and 2016, respectively. However, the variation in the likelihood of being married, between subjects was quite substantial [which was allowed to vary with a standard deviation of 0.22 (CI=0.21; 0.23)].

For Separated vs. Never Married state occupation, although there is heterogeneity among subjects (with a standard deviation of the heterogeneity component of 0.09(CI = 0.08; 0.10)), males and those without employment were less likely to be in a separated state than females and those who were employed

\**Statistically significant variable.*

[OR = 0.42(CI = 0.24; 0.73) and OR = 0.52(CI = 0.35; 0.78) for gender and employment, respectively]. When compared to those who never went to school, the more educated a subject became, the higher the likelihood of being in a Separated state [OR = 1.62(CI = 0.59; 4.46), 4.50(CI = 1.52; 13.34) and 14.86(CI = 4.14; 53.31) for primary, secondary, and tertiary education, respectively].

When comparing subjects who occupied a Widowed state relative to those who never married, males were less likely to be widowed compared to females [OR = 0.09(CI = 0.07; 0.11)]. When compared to subjects without any education, the likelihood of occupying a widowed state significantly increased with increasing level of education [OR = 1.34(CI = 1.11; 1.62), 1.62(CI = 1.29; 2.02), and 2.96(CI = 1.96; 4.46) for primary, secondary, and tertiary education, respectively]. There existed some significant subject-to-subject variation in the likelihood of the Widowed vs. Never Married state occupation [standard deviation of the heterogeneity component, OR = 0.431(CI = 0.231; 0.631)].

## 3.3. Results for Single Marital State Transitions

Results of the binary transitions are displayed in **Table 4** below. These transitions include marital dissolution, re-marrying after a dissolved marriage (Separated or Widowed →Married), marital dissolution due to a separation as well as a dissolution as a result of death of a partner. The hazards of transitions over the years are displayed in **Figure 5**. For each transition type, a separate baseline hazard was allowed. Using the smallest BIC, they all had a linear baseline effect of age except for the married→separated transition whose baseline was constant. These have been represented

TABLE 4 | Results for single transitions for the different family dynamics for subjects aged between 17 and 65 years between 2004 and 2016 for the 56,308 subjects.

\**Statistically significant variable.*

graphically in **Figure 2**. Additionally, there was homogeneity among subjects' re-marriage rates, since the standard deviation of the heterogeneity component was estimated as SD = 0. Increasing age significantly reduced the rate of re-marriage [OR = 0.92(0.88; 0.97)]. Males were significantly less likely to experience a marital dissolution than females [OR = 0.27, CI = (0.17; 0.45)]. Raising age at first sex was associated with a lower rate of marital dissolution [OR = 0.97(CI = 0.95; 0.99)]. The likelihood of having a marital dissolution was significantly associated with recent periods, with the period 2008 − 2009 having the highest effect [OR = 24; CI = (5.53 − 108.37)] [standard deviation of 1.44, CI = (1.04, 2.01)].

Males were significantly less likely to experience a Married to Widowed transition, when compared to females [OR = 0.20, CI = (0.12; 0.36)]. The rest of the covariates did not significantly affect transition into widowhood. However, there existed some variation between subjects which is unaccounted for in the analysis [standard deviation of 1.50(CI = 1.07; 2.09)]. Lastly, higher ages at first marriage were associated with lower rates experiencing a Married→Separated transition for those who entered first marriages before the age of 40. The likelihood of a Married→Separated transition slightly decreased with increasing age at first marriage [OR = 0.55, CI = (0.10; 2.94) and OR = 0.38, CI = (0.06; 2.56) for age groups 23 − 29

and 30 − 40, respectively]. The standard deviation for random effects was 3.43(CI = 2.50; 4.72) which implies that there could be other factors associated with transitions from Married to

TABLE 5 | Results for exiting a marriage for subjects aged between 17 and 65 years since 2004–2016 from a population of 56,308 subjects.

\**Statistically significant variable.*

Separated state which were not captured by the model that causes the variation.

## 3.4. Results for Termination of a Marriage

Results of the competing risks model using the multivariate binary response are displayed in **Table 5** below. The hazards of transitions over the years are displayed in **Figure 6**. These concern transitions out of a Married state, whose competing destination marital states were Separated and Widowed. The Married→Separated transition had a constant baseline age effect while the Married→Widowed had a linear age effect. Marrying at later ages was associated with a lower likelihood of transition from married to widowhood. [OR = 0.65, CI = (0.43; 0.99) and OR = 0.40, CI = (0.19; 0.88) for those who entered first marriage between 30 − 40 and ≥ 41 years, respectively]. Males were significantly less likely to experience a Married→Widowed transition when compared to women [OR = 0.20, CI = (0.12; 0.36)]. Substantial heterogeneity [OR of SD = 3.56(CI = 2.66, 4.77) and SD = 1.50(CI = 1.08, 2.09)] between subjects on the Married → Separated and Married → Widowed transitions, respectively, were observed.

## 4. DISCUSSION

Marriage and marriage dissolution are important life events that affect adults in a society. We used the rich data from one of Africa's largest population-based cohorts based in KwaZulu-Natal South Africa to model family formation and dissolution using a series of multilevel discrete time event models. In line with Garenne (2004) and Hosegood et al. (2009), our results demonstrate that marital rates are low (13.6% females and 10.3% males) and that age at first marriage is remarkably high in this population (median of 34 years). In addition, the results of the multilevel discrete time to event models revealed that marrying at later ages had a clear association with a low rate of marital dissolution, whilst more educated subjects and later age of sexual debut were associated with a lower likelihood of experiencing a marital dissolution. The role of these factors on marriage warrants further discussion and may be investigated in future research which can be supported by data.

The study has several strengths including the use of rich data from one of Africa's largest population-based cohorts and stateof-the-art statistical methods. However, even more advanced statistical models can still be used to model the transitions in this process simultaneously. As such, in addition to modeling the hazards and covariate effects, other methods may then be used to determine more useful statistics, such as expected waiting times, sojourn times, and transition probabilities which may then be used for prediction. Additionally, the random effects in this study were modeled using a normal distribution but a different distribution (such as gamma mixture, by Jenkins, 1997) may also be considered where different assumptions are made about the subjects. Different approaches may also be used to handle missing data.

Although marriages may get terminated as a result of poor (psychological) health (Wang and Amato, 2000), marital dissolution might lead to undesired health outcomes on both partners and children, such as stress, and high blood pressure. We found significant factors which contribute to a marital dissolution which may help with policy decision making. Many prevention programs in HIV encourage delay in sexual debut in order to reduce risk of HIV acquisition in the highly vulnerable youth (Karim et al., 2017). Thus, programs which encourage delay in sexual debut may have additional benefits in terms of reducing rates of marital dissolution thus ultimately improving psychological and physical health.

## DATA AVAILABILITY STATEMENT

The datasets for this study can be found upon request from Africa Health Institute whose website is https://www.ahri.org/. Data sets with restricted access require a data access agreement to be completed. A request is then submitted to the applicable data custodian for specific data sets on the repository. For requests of ad hoc data sets, beyond the data sets archived on the Africa Centre data repository, these have to be directed to Africa Centre's Helpdesk (help@africacentre.ac.za).

## ETHICS STATEMENT

The studies involving human participants were reviewed and approved by Ethical approval for all the data collected by AHRI was obtained from the University of KwaZulu-Natal's Ethics Committee (BE 169/15). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

## AUTHOR CONTRIBUTIONS

JB performed the literature review, data management, statistical analysis, and wrote the initial draft of the manuscript. SM conceived, designed, and developed the methodology of the study. HM helped with the design and methodology of the study. SM and HM reviewed the statistical analysis and helped with revision of the manuscript. FT contributed to the contextualization of the study findings and substantively helped with the revision of the manuscript. JB, SM, HM, and FT read and approved the final version of the manuscript.

## FUNDING

This work was funded and supported through the DELTAS Africa Initiative [SSACAB]. The DELTAS Africa Initiative was an independent funding scheme of the African Academy of Sciences (AAS)'s Alliance for Accelerating Excellence in Science in Africa (AESA) and supported by the New Partnership for Africa's Development Planning and Coordinating Agency (NEPAD Agency) with funding from the Wellcome Trust [grant 107754/Z/15/Z-DELTAS Africa Sub-Saharan Africa Consortium for Advanced Biostatistics (SSACAB) programme] and the UK government. The views expressed in this publication are those of the author(s) and not necessarily those of AAS, NEPAD Agency, Wellcome Trust or the UK government. We also received funding from University of KwaZulu- Natal's University Capacity Development Programme (UCDP) and the Teaching development grant [APP-TDG-226], both from South African Department of Higher Education and Technology (DHET) through University of Pretoria. We acknowledge the Africa Health Research Institute's Demographic Surveillance Information System and Population Intervention Programme which was funded by the Wellcome Trust (201433/Z/16/Z) and the South Africa Population Research Infrastructure Network (funded by the South African Department of Science and Technology and hosted by the South African Medical Research Council).

## ACKNOWLEDGMENTS

We gratefully acknowledge the AHRI for providing with the datasets. The authors were also grateful to the study participants and the work and support of the fieldwork and database teams at AHRI. The Biostatistics Research Unit of the South Africa Medical Research Council provided research visits to their Pretoria offices. We also acknowledge the support we received from the University of KwaZulu-Natal through the provision of research facilities.

## REFERENCES

death: a rural South African population-based surveillance study. PLoS Med. 10:e1001409. doi: 10.1371/journal.pmed.1001409

Hashemi, L., and Homayuni, H. (2017). Emotional divorce: child's well-being. J. Divorce Remarriage 58, 631–644. doi: 10.1080/10502556.2016.1160483

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Batidzirai, Manda, Mwambi and Tanser. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Progression of HIV Disease Among Patients on ART in Ethiopia: Application of Longitudinal Count Models

### Belay Desyebelew Andualem1 \* and Birhanu Teshome Ayele2,3

*1 Department of Statistics, Dire Dawa University, Dire Dawa, Ethiopia, 2 Department of Statistics, Addis Ababa University, Addis Ababa, Ethiopia, 3 Division of Epidemiology and Biostatistics Unit, Faculty of Medicine and Health Sciences, Stellenbosch University, Stellenbosch, South Africa*

Although the world has been fighting HIV disease in unity and patients are getting antiretroviral therapy treatment, HIV disease continues to be a serious health issue for some parts of the world. A large number of AIDS-related deaths and co-morbidities are registered every year in resource-limited countries like Ethiopia. Most studies that have assessed the progression of the disease have used models that required a continuous response. The main objective of this study was to make use of appropriate statistical models to analyze routinely collected HIV data and identify risk factors associated with the progression of the CD4+ cell count of patients under ART treatment in Debre Markos Referral Hospital, Ethiopia. In this longitudinal retrospective study, routine data of 445 HIV patients registered for ART treatment in the Hospital were used. As overdispersion was detected in the data, and Poisson-Gamma, Poisson-Normal, and Poisson-Gamma-Normal models were applied to account for overdispersion and correlation in the data. The Poisson-Gamma-Normal model with a random intercept was selected as the best model to fit the data. The findings of the study revealed the time on treatment, sex of patients, baseline WHO stage, and baseline CD4+ cell count as significant factors for the progression of the CD4+ cell count.

Keywords: HIV/AIDS, CD4 count, longitudinal data, Poisson-Normal model, Poisson-Gamma-Normal model, antiretroviral therapy (ART), Ethiopia, Debre Markos

## 1. INTRODUCTION

HIV disease continues to be a serious health issue for resource-limited countries like Ethiopia. According to the UNAIDS (2016) fact sheet, there were about 2.1 million new cases of HIV in 2015 globally (1). About 36.7 million people were living with HIV around the world, and, as of June 2016, 18.2 million people living with HIV were receiving medicine to treat HIV, called antiretroviral therapy (ART). An estimated 1.1 million people died from AIDS-related illnesses in 2015, and 35 million people have died from AIDS-related illnesses since the start of the epidemic. CD4+ cell counts are the primary targets of HIV. The relentless destruction of CD4+ cell counts by HIV, either directly or indirectly, results in the loss of HIV-specific immune responses and, finally, non-specific immune response in the AIDS stage. The estimation of peripheral CD4+ cell counts has been used as a tool for monitoring disease progression and the effectiveness of antiretroviral treatment (ART) (2). The changes in the CD4+ cell counts are important indicators of the response to ART.

### Edited by:

*Jim Todd, University of London, United Kingdom*

### Reviewed by:

*Edson Mollel, Kilimanjaro Christian Medical University College, Tanzania Kathryn Risher, University of London, United Kingdom*

### \*Correspondence:

*Belay Desyebelew Andualem<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Digital Public Health, a section of the journal Frontiers in Public Health*

Received: *20 May 2019* Accepted: *23 December 2019* Published: *19 February 2020*

### Citation:

*Andualem BD and Ayele BT (2020) Progression of HIV Disease Among Patients on ART in Ethiopia: Application of Longitudinal Count Models. Front. Public Health 7:415. doi: 10.3389/fpubh.2019.00415* Initial CD4+ cell count, age, gender, smoking, unemployment, WHO stage, hospital, opportunistic infections, body mass index, changing doctors during outpatient follow up, use of alcohol and drugs, and duration of treatment (in months) are some of the significant determinants that affect CD4+ cell count progression of patients on ART (3–5).

Most studies conducted in the area fitted statistical models that require (multivariate) Normal distribution by considering CD4+ cell counts as continuous variable. When this assumption is violated, even after transformation, considering Poissonrelated models is a natural choice. One of the common problems one can be faced with in analyzing count data like the CD4+ cell count is overdispersion. A Negative Binomial model can be considered to overcome this problem. Trindade et al. applied Poisson and Negative Binomial models using the multilevel (ML) approach and the generalized estimations equations (GEE) to model CD4+ cell counts of 587 HIV seropositive patients, and they stated that the best marginal model to fit the data was the Negative Binomial (NB) with an exchangeable correlation structure (6). Tekle et al. also employed different count data analysis methods starting from the ordinary Poisson regression model to study CD4+ cell counts of 222 HIV positive patients, and they found that Poisson-Normal-Gamma is the best model to fit their data (7). In this study, we applied various count data models to study the progression of the CD4+ cell count of HIV patients and identified risk factors for progression of patients' CD4+ cell count in Debre Markos Referral Hospital, Ethiopia.

## 2. MATERIALS AND METHODS

In practice, it is common to have response variables of a count type-like number of the CD4+ cell count in a cubic milliliter of blood. Some data analysts treat the CD4+ cell count as a continuous measure and apply the linear mixed effects model. But that practice ignores two facts: the data are really discrete, and the distributions of count variables are usually skewed. For these reasons, the use of models that assume (multivariate) normality might not be efficient (8). Even if the data is transformed and these models are applied, the interpretation might not be straightforward. In scenarios like this, it is better to apply statistical models that account for the nature of the data.

Our data includes 445 HIV-positive patients who started ART treatment between December 2005 and July 2014 in Debre Markos Referral Hospital, Ethiopia. The minimum number of measurements was two and the maximum was seven. Patients

with less than two measurements and age of <15 years were excluded from the study.

For our data, the assumption of multivariate normality failed, and this suggested that use of a linear mixed model was not appropriate (**Table 1**). The Poisson regression model with normal random effects and models that account for both correlation between repeated measures and overdispersion simultaneously were thus considered in line with Booth et al. (9) and Molenberghs et al. (10, 11).

## 2.1. Variables in the Study

### 2.1.1. Dependent Variable

The dependent variable of this study was the CD4+ cell count per cubic millimeter of blood of HIV-infected patients who are under ART treatment.

### 2.1.2. Independent Variables

The independent variables considered in this study were selected based on related literature (5, 7). These include the sex of patients, age of patients (age at the initiation of the treatment), baseline CD4+ cell count (the CD4+ cell count of the patients at the start of the treatment), WHO clinical stage at baseline (stage I, stage II, stage III, and stage IV), marital status at baseline, baseline weight, level of education at baseline, functional status at baseline, TB status at baseline, and time in months. Functional status was defined as WHO categories: Ambulatory and Working. Patients who are able to perform activities of daily living but not able to work or play are classified as ambulatory and the who are able to perform usual work in or out of the house, harvest, go to school or for children, normal activities, or playing were classified as working.

## 2.2. Poisson Model

Let Yi be the ith CD4+ cell count and is Poisson distributed with mean λi . The density function of Yi can then be written as

$$f(Y\_i = y\_i | \lambda\_i) = \frac{e^{-\lambda\_i}\lambda\_i^{y\_i}}{y\_i!} = \exp\{y\_i \ln \lambda\_i - \lambda\_i - \ln y\_i!\}, \quad \text{(1)}$$

The Poisson distribution belongs to the exponential family, with natural parameter θi equal to ln λi , scale parameter φ = 1, and variance function v(λi) = λi (12). The logarithm is the natural link function, leading to the classical Poisson regression model Yi ∼ Poisson(λi), with log(λi) = X T i β.

## 2.3. Poisson-Gamma Model

The standard Poisson distribution requires the mean and variance to be equal. When this assumption fails, the Poisson-Gamma model should be used to fit the data. Assume that Yi |θi ∼ Poi(θiλi), where θi denotes an independent and identically distributed (iid) sample of unit mean Gamma random variables with shape parameter α (9). Conditional on θi , the CD4 count of the ith patient follows a Poisson distribution with mean θiλi . The counts are then marginally independent Poisson-Gamma random variables [Yi ∼NB(α, λi)] with mean λi and variance λi + λ 2 i /α. Hence, the parameter α quantifies the amount of overdispersion with α = ∞ corresponding to no overdispersion

### TABLE 2 | Summary of CD4+ cell count at different time points.

TABLE 3 | Summary of CD4+ cell count progression for some categorical covariates.

with respect to the Poisson distribution. The mass function of the Poisson-Gamma random variables is given by

$$Pr(Y\_i = \boldsymbol{\wp}; \alpha, \lambda\_i) = \frac{\Gamma(\boldsymbol{\wp} + \alpha)}{\Gamma(\alpha)\boldsymbol{\wp}!} (\frac{\alpha}{\lambda\_i + \alpha})^{\alpha} (\frac{\lambda\_i}{\lambda\_i + \alpha})^{\boldsymbol{\wp}} \tag{2}$$

The Poisson-Gamma model (also known as the Negative Binomial model) is given by log(λi) = X T i β.

## 2.4. Poisson-Normal Model

For µij =E(Yij|bi) and known link function η(.), the generalized linear mixed model can be expressed as:

$$\eta(\mu\_{i\bar{j}}) = \eta[E(Y\_{i\bar{j}}|b\_i)] = X\_{i\bar{j}}^T \beta + Z\_{i\bar{j}}^T b\_i \tag{3}$$

where Yij is the CD4+ cell count of the ith patient at jth visit (measurement). β= a p-dimensional vector of unknown fixed regression coefficients. bi = a q-dimensional vector of unknown random regression coefficients for the ith individual, and these are often assumed to be drawn independently from the N(0, D), and D is the variance-covariance matrix of the random effects. Xij and Zij are p-dimensional and q-dimensional vectors of known covariate values, respectively (10). The generalized mixed Poisson model with normal random effects (Poisson-Normal model) becomes

$$\ln(\lambda\_{i\bar{j}}) = X\_{i\bar{j}}^T \beta + Z\_{i\bar{j}}^T b \mathbf{i} \tag{4}$$

This model is referred to as the Poisson-Normal model because it assumes Poisson distribution for the counts and normal distribution for the random effects bi (10, 11).

## 2.5. Poisson-Gamma-Normal Model

According to Molenberghs et al. (10, 11), a model combining the ideas from the Poisson-Normal and overdispersion models for repeated Poisson data with overdispersion can be specified as follows Yij ∼ poi(θijλij)

$$\lambda\_{i\dot{\jmath}} = \exp(X\_{i\dot{\jmath}}^T \beta + Z\_{i\dot{\jmath}}^T b\_i) \tag{5}$$

where θij capture overdispersion and denote an independent and identically distributed (iid) sample of unit mean gamma random variables with shape parameter α and scale parameter β=1/α, and where bi ∼ N(0, D) and θij ∼ Gamma(α, β). This model is called the Poisson-Gamma-Normal (combined) model because it includes both Normal (bi) and Gamma (θij) random effects to account for correlation and overdispersion, respectively.

## 2.6. Methods of Parameter Estimation

In this study, we used glmer and glmer.nb functions in R under packages MASS and lme4. A Laplace approximation was used to obtain parameter estimates. The R code used to fit the models is available in **Supplementary Material**.

## 2.7. Model Comparison

To select the important variables, first the main effect, main effect by time interaction, and plausible main effect by main effect interactions were incorporated to the initial candidate models,

the non-significant interaction effects were then removed, and the models were refitted again and so on. The best model that can fit the data was selected using various information criteria (AIC, BIC, and −2loglikelihood) (**Table 7**). The model with smallest values of information criteria was selected as the final model.

## 3. RESULTS AND DISCUSSION

## 3.1. Descriptive Analysis

In this section, CD4+ cell count data obtained from 445 HIV patients on ART treatment in Debre Markos Referral Hospital were summarized. The majority of the HIV patients [347 (78.0%)] started antiretroviral treatment with CD4+ cell counts <200 cells/mm3 . At the start of the treatment, the median CD4+ cell count of the patients was 145 CD4+ cells/mm3 of blood with IQR of 107.00 CD4+ cells/mm3 of blood. The minimum and maximum baseline CD4+ cell counts were three and 971 CD4+ cell cells/mm3 of blood, respectively.

The summary of CD4+ cell counts at different time points is given in **Table 2**. As can be seen in **Table 2**, the median CD4+ cell count increased over time. The IQR of CD4+ cell counts increased at some points and then started to decrease after the 24th month. The number of patients decreased at some points and increased at others, which implies the presence of intermittent missingness in the data. That means some patients were falling out of care and then re-engaging, or they did not have CD4+ cell counts that were spaced perfectly every 6 months.

Data on demographic and clinical characteristics of the patients was collected at the start of antiretroviral treatment. Among the 445 patients, 280 (62.9%) were females. The male patients had a 134.84 mean baseline CD4+ cell count, while the female patients had a mean baseline CD4+ cell count of 168.91. On average, female patients started ART treatment at a relatively higher CD4+ cell count. The difference in mean CD4+ cell count of the two groups increases as time increases. The average CD4+ cell count of females was higher than males at all time points and the difference increases over time.

WHO stage III had a higher number of patients [282 (63.4%)] as compared to the other three stages. WHO stage II took second place in number of patients [77 (17.3%)], and WHO stage IV had the smallest number of patients [27 (6.1%)]. As expected, patients on WHO stage I had a higher CD4+ cell count at all time points as compared to patients of the other three stages of the disease.

Patients with a working functional status have a higher mean CD4+ cell count at all time points than that of patients with ambulatory functional status. Among the 445 HIV patients included in this study, 346 (77.8%) were patients with working functional status and 99 (22.3%) were ambulatory (**Table 3**). About 27.3% of the 445 HIV patients were TB positive at baseline. TB negative patients had a higher mean CD4+ cell count at

all time points compared to TB positive patients, implying the impact of the HIV-TB coinfection.

## 3.2. Exploratory Data Analysis

**Figure 1** depicts the individual profile plot of the CD4+ cell count of HIV-infected patients included in the study. The plot provides some information on the between patients' CD4+ cell count variability and illustrates the over-time change in patients' CD4+ cell count. Some individuals have an erratic CD4+ cell count and others have a CD4+ cell count that slowly increases over time. As one can see from the graph, there is a considerably large difference in the intercepts of individual trajectories. Similarly, some trajectories are steeper, while others were almost horizontal, indicating the possible variability in the slope of CD4+ cell counts. Therefore, because of the variability in the intercept and slope of trajectories, using a mixed model could fit the data very well. The overall mean profile plot of the CD4+ cell count shows somehow a linear increasing pattern of CD4+ cell count over time (**Figure 2**), suggesting that a linear time effect seems reasonable. The mean CD4+ cell count

TABLE 5 | Poisson-Normal model.

increases at a high rate from baseline till the 6th month and then starts to increase slowly from 6 to 24th month and decreases at month 30.

## 3.3. Model Results

**Table 4** summarizes the parameter estimates of Poisson and Poisson-Gamma regression models employed on the CD4+ cell count. All parameters included in the Poisson regression model are significant at 5% level of significance. For this study, the data were overdispersed, as the sample variance of CD4+ cell count at all time points was greater than its corresponding sample means (**Table 2**). A likelihood ratio (LR) test was used to test the null hypothesis that the restriction in the Poisson model was true. The test revealed that the null hypothesis was rejected, implying the presence of overdispersion in our data. The Poisson-Gamma model better fits the data as compared to the Poisson model with smallest AIC value. The Poisson-Normal model with both random intercept and slope was found to be the best fit since it has smaller information criteria values compared with the only random intercept model. The parameter

estimates of this model are displayed in **Table 5**. Depending on this model time, the WHO stage and initial CD4+ cell count were found to be significant factors of patients' CD4+ cell count progression.

An improvement in both the Poisson-Gamma and Poisson-Normal models as compared with the Poisson model in fitting the data is an indication of the occurrence of both correlation and overdispersion in the data. The Poisson-Gamma-Normal (Negative Binomial log-linear mixed) model proposed by Booth et al. (9) and Molenberghs et al. (10, 11) was fitted to overcome this problem of correlated and overdispersed count data, and the random intercept Poisson-Gamma-Normal Model is a much better fit because of its lower AIC (27,379.9), BIC (27,488.4), and −2loglikelihood (27,342) values as compared to the Poisson-Normal models (**Table 6**). Therefore, the final model to fit our data was the random intercept Poisson-Gamma-Normal model. We have also tried the Poison-Gamma-Normal model with different (random) linear slopes for a time, but we found that the Poison-Gamma-Normal with random intercept was better based on information criteria (AIC and BIC).

Based on the results obtained from the Poisson-Gamma-Normal model, time in months, sex, and baseline CD4+ cell count were found to be significant factors of the CD4+ cell count of a patient (**Table 8**). For a given patient, keeping the random intercept and other covariates constant, one more month on ART increased the CD4+ cell count by a multiplicative factor of e 0.0243 = 1.0246.

A female patient had a CD4+ cell count of 1.1215 times that of a male patient, adjusting for other covariates and random intercept. A unit change in baseline CD4+ cell count increased the CD4+ cell count of a patient by a factor of 1.0034, fixing the values of the other covariates and the random intercept constant.

The dispersion parameter (1/α) has been estimated, in the final model, as 7.7009, and the Gamma (overdispersion) random effects are assumed to follow a Gamma distribution with unit mean and shape parameter α (0.130).

TABLE 6 | Poisson-Gamma-Normal model.

TABLE 7 | Summary of information criteria of different models.

TABLE 8 | Poisson-Gamma-Normal model.

## 3.4. Discussion

The effects of demographic and clinical factors on the progression of CD4+ cell counts over time of HIV patients taking ART treatment in Debre Markos Referral Hospital were assessed using Poisson longitudinal models since the response variable of interest CD4+ cell count is a count variable.

The results of the summary statistics revealed that the value of IQR is high at all time points, which might be an indication for high variation among the patients' CD4+ cell count at baseline as well as at different time points after the initiation of ART treatment. This variation might have been caused by the year at which the patients started ART treatment, as there have been different WHO's CD4+ cell count cut-off points to initiate ART treatment at different times. Although most of the patients included in our study started with lower CD4+ cell counts (<200 cells/mm3 ), there were patients who had higher baseline CD4+ cell counts (971 cells/mm3 ). Despite the continuous effort to initiate early, some patients still presented with lower CD4+ cell counts, which might be due to patients' lack of willingness to get tested (13, 14) or difficulties to provide treatments to all patients in lower-income countries including Ethiopia. Hence, we believe that our result could be generalizable. The final model also indicated that initial CD4+ cell count (CD4+ cell count at the start of the treatment) significantly affects CD4 count progression. Therefore, based on our findings we recommend patients to start the treatment early as of the WHO's "treat all" recommendation.

The sign of the parameter estimate of WHO stage III is positive, which implies that a patient with WHO stage III has a higher CD4+ cell count as compared with a patient of WHO stage II. It might be because the number of patients with WHO stage III are much higher (non-comparable) than patients with WHO stage II. The relationship between CD4+ cell count and WHO stage III might also be explained by the baseline CD4+ cell count. Duration of treatment also have a positive effect on the CD4+ cell count progression of HIV patients. This means patients with longer time on ART treatment have good recovery of CD4+ cell count than that of patients with short duration on the treatment.

## 4. CONCLUSION

An analysis of CD4+ cell count data using conventional models like linear mixed models might be inadequate as the data were highly skewed and may not satisfy normality (multivariate) assumption as demonstrated in our data.

In this study, CD4+ cell count data of 445 HIV patients under ART in Debre Markos Referral Hospital was analyzed using different longitudinal count models, and the Poisson-Gamma-Normal model was selected as the final model to fit the data based on different selection criteria. The Poisson-Gamma-Normal model handles overdispersion and correlation simultaneously.

The duration on ART treatment (time in months), sex of patients, and baseline CD4+ cell count were all identified as potential risk factors of CD4+ cell count progression. Having a good CD4+ cell count at baseline had a positive impact on CD4+ cell count evolution over time.

Although good CD4+ cell count progress in response to ART was observed, most of the patients (78.0%) were at decreased CD4+ cell counts (<200 cells/mm3 ) when enrolled for ART treatment, which might have contributed to low CD4+ count recovery in some patients.

## 5. LIMITATIONS AND RECOMMENDATION

In our study, we only considered patients from one hospital. The likelihood inference of the models considered in this study are valid under MCAR (missing completely at random). In the current study, we did not carry out a sensitivity analysis, and we only considered linear slopes models, although a different linear slope for different time periods seems reasonable. Hence, we recommend that researchers consider sensitivity analysis and data obtained from different Hospitals. The age and weight of patients might have a non-linear relationship with the CD4+ cell count. We recommend smoothing techniques like splines to be explored for further studies. The assumption of multivariate normality that is assumed by most statistical models used in longitudinal data analysis should be checked before analysis. Efficient methods like the ones used in this study could be considered if the assumption is violated.

## DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

## ETHICS STATEMENT

Before data collection, a letter of support written by the Statistics Department of Addis Ababa University was submitted to Debre Markos Hospital and permission to collect anonymized data was obtained. The data was extracted by trained data clerks in the ART Clinic and none of the researchers had access to original cards of patients. Written informed consent for participation

## REFERENCES

was not required for this study in accordance with the national legislation and the institutional requirements.

## AUTHOR CONTRIBUTIONS

BAn conceived the idea, performed the data cleaning and analysis, interpreted the ensuing results, and drafted the manuscript. BAy supervised the study, contributed to the conception, and revised the manuscript. Both the authors read and approved the final draft.

## ACKNOWLEDGMENTS

We acknowledge the ART case unit information center of Debre Markos Referral Hospital, Ethiopia for the data they supplied. The manuscript was prepared on the basis of my Masters' thesis.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpubh. 2019.00415/full#supplementary-material

R codes used to fit the models are attached as supplemental data.

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Andualem and Ayele. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Survival of Children Living With HIV on Art in Zambia: A 13-Years Retrospective Cohort Analysis

Tendai Munthali 1,2 \*, Charles Michelo1 , Paul Mee3 and Jim Todd4

*1 School of Public Health, University of Zambia, Lusaka, Zambia, 2 Department of Public Health, Ministry of Health, Lusaka, Zambia, 3 MeSH Consortium, Department of Public Health Environments and Society, Faculty of Public Health and Policy, London School of Hygiene and Tropical Medicine, London, United Kingdom, 4 Department of Population Health, London School of Hygiene and Tropical Medicine, London, United Kingdom*

Background: Research conducted before the introduction of anti-retroviral therapy (ART), showed that the majority of children living with HIV (CLHIV) would die before their second birthday. In Zambia, ART was rolled out to the public health system in 2004 with subsequent improved survival in CLHIV. However, the survival rates of CLHIV on ART in Zambia since 2004 have not been extensively documented. We assessed survival experiences and the factors associated with survival in CLHIV on ART in Zambia.

### Edited by:

*Philip Ayieko, University of London, United Kingdom*

### Reviewed by:

*Samuel Manda, South African Medical Research Council, South Africa Lucas Malla, KEMRI Wellcome Trust Research Programme, Kenya*

> \*Correspondence: *Tendai Munthali<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Digital Public Health, a section of the journal Frontiers in Public Health*

Received: *10 May 2019* Accepted: *09 March 2020* Published: *31 March 2020*

### Citation:

*Munthali T, Michelo C, Mee P and Todd J (2020) Survival of Children Living With HIV on Art in Zambia: A 13-Years Retrospective Cohort Analysis. Front. Public Health 8:96. doi: 10.3389/fpubh.2020.00096* Methods: We conducted a retrospective cohort analysis of CLHIV (aged up to 15 years) using routinely collected data from health facilities across Zambia, over 13 years to ascertain mortality rates. We explored survival factors using Cox regression giving adjusted hazard ratios (AHR) and 95% confidence intervals (95% CI). Nelson Aalen estimates were used to show the cumulative hazards of mortality for different levels of explanatory factors.

Results: A total of 65,448 eligible children, were initiated on ART between 2005 and 2018, of which 33,483 (51%) where female. They contributed a total survival time of 275,715-person years at risk during which 3,265 children died which translated into an incidence rate of 1.1 deaths per 100 person-years during the review period. Mortality rates were highest in children in the first year of life (Mortality rate 2.24; 95% CI = 2.08–2.42) and during the first year on ART (Mortality rate 3.82 95% CI = 3.67–3.98). Over 50% of the children had been on ART for 5–10 years by 2018, and they had the lowest risk of mortality compared to children who had been on ART for <5 years.

Conclusions: Children with HIV in Zambia are surviving much longer than was predicted before ART was introduced 14 years ago. This key finding adds to the literature on analysis of survival in CLHIV in low income settings like Zambia. However, this survival is dependent on the age at which ART is initiated and the time on ART highlighting the need to increase investments in early infant diagnosis (EID) to ensure timely HIV testing and ART initiation for CLHIV.

Keywords: survival, HIV- infected, children, Zambia, anti-retroviral treatment (ART), initiation

## INTRODUCTION

UNAIDS has a target to end all new HIV infections by 2030, but in 2018 there were still 1.7 million children living with HIV (CLHIV) (1). Worldwide, every year there are 110,000 deaths in CLHIV under the age of 15 years, and 180,000 new HIV infections from mother to child transmission (1, 2). In Zambia, in 2018 there were ∼62,000 children living with HIV, 79% of them were on ART and there were 3,000 AIDS related annual deaths (3, 4). Without ART, more than 50% of children infected with HIV will die before their second birthday (5, 6). From its inception ART has been shown to reduce mortality among adults and children living with HIV more than 50% since the year 2000 (7–9). With the roll out of HIV treatment services, the initiation of CLHIV on ART in Zambia increased steadily from 24,000 in 2010 to 49,116 in 2018 (2, 3). This translated into a reduction in HIV-related mortality from 3,600 in 2016 to 3,000 in 2018 (3, 4).

However, studies show that mortality among CLHIV on ART continues to be high, especially among infants with advanced disease, and in the first 12 months following ART initiation despite the use of prevention of mother-to-child transmission of HIV (PMTCT) drugs (10, 11). Studies have shown that initiating children on ART early (birth-7 weeks) reduced mortality by 76% compared to postponing ART until progression to advanced HIV disease or CD4 count thresholds are reached (11, 12). However, other studies have shown that WHO clinical stage 3 or 4 and young age at the time of ART initiation were associated with a 5-fold increased risk of death, while absolute CD4 counts of <350 cells/mm3 at ART initiation and underweight (weightfor-age Z-score <-2) were associated with 3.6 and 3,5 times increase in mortality among CLHIV (5, 13–15). In addition, overall mortality rates among CLHIV on ART vary from country to country depending on several factors. In a 7-year study in South Africa, a mortality rate of 4.7 per 100 person-years was reported (16), whereas in Ethiopia a mortality rate of 1.2 per 100 person-years was reported in Addis Ababa and 2.1 per 100 person-years in Eastern Ethiopia after a 2 year follow up (14). Furthermore, a mortality rate of 3.0 per 100 person-years was reported in a 1 year follow up study in Nigeria (17).

Despite the varied mortality rates, CLHIV are surviving longer on ART. Studies in Asia and Africa have shown survival rates that ranged from 84 to 97% after 12 months of ART for children initiated on ART at a median age of 5–7 years. Another multi country study in Africa also reported above 90% (93%) survival rates after 2 years of follow up on ART (18). There is limited global data on long term survival of children on ART. In Zambia the long-term survival of CLHIV and clinical factors that are associated with patterns of survival have not been extensively documented. We therefore assessed the probability of survival and factors associated with mortality among CLHIV in Zambia.

## MATERIALS AND METHODS

## Study Population and Sample

This study was a quantitative retrospective cohort study of CLHIV with records in the Zambian SmartCare data. SmartCare is an electronic patient monitoring system which is used in about 600 government health facilities in all districts in Zambia (19). The system records patient characteristics with a unique identification number, at the first contact with the clinic, and then records clinical information at each subsequent health facility visit. A total sampling from the SmartCare records of all children with a positive HIV result, aged 15 years and younger at the time of ART initiation, regardless of missing data during the period under review, was conducted for this analysis.

## Clinical Procedures

In the Zambian health care system, all children with a positive HIV test result and diagnosis are eligible to receive HIV care according to national ART prevention and treatment guidelines. This includes documenting medical history, physical examination, anthropometric measurements, socio demographic information, and WHO clinical stage. In addition, measurement of CD4+ T-cell counts or percentages, hemoglobin levels (HB), renal, liver function tests, and HIV viral loads are taken. Eligible children are then treated with a first-line regimen, either Nevirapine based (NVP), Efavirenz based (EFV), or Nucleoside reverse transcriptase inhibitors (NRTI). Children are then asked to return for clinical evaluation every 3 months (20, 21).

## Data Definitions and Statistical Methods

Follow up of children started at the date they initiated ART, and continued until death, or the date of censoring. The event of interest, death was coded as 1 for those who died, and 0 for those who were censored. The date of censoring (end date) contained the date of death for those that died, and the last clinic date for those who were lost to follow up, with all other children censored on 1st January 2018, or the date of their 15th birthday, whichever came first. The variable age at ART initiation was obtained by subtracting the child's date of birth from the date of ART initiation and was categorized into 5 age bands (birth to <1 year, 1–2, 3– 5, 6–9, and 10–15 years). Time on ART was divided into 4 points that is at 1, at 2, at 5, and at 10 years on ART. Children's weight was used to calculate the weight for age Z-scores (WAZ) for children in the dataset using a WHO standard (22).

Categorical variables were summarized across all children as frequencies and percentages. Mortality rates were defined as the number of deaths divided by the person-years of followup with 95% confidence intervals (95% CI) for mortality rates defined for each subgroup of children. Nelson Aalen cumulative Hazards were used to graphically show the differences in mortality by subgroups of children. A modified Cox regression model, adjusting for possible heterogeneity at the facility level, was used to obtain hazard ratios (HR) and 95% confidence intervals (95% CI) for factors associated with mortality, with subsequent adjustment for sex, age at ART initiation, and province to obtain adjusted HR (aHR) and 95% CI. In order to model the unobserved covariates at health facility level, an additional parameter, theta, drawn from a Gamma distribution, was estimated in the modified Cox regression model, representing the frailty or variation in the log of the survival function across health facilities.

For socio-demographic variables, such as age, sex, year, health facility and province, a complete case analysis was conducted, using all the children, assuming the small amount of missing data were missing at random. The results present the crude analysis, the age-adjusted analysis and the frailty model that adjusts for the effect of clustering. For variables where a large number of children had missing data, a subset analysis on children

TABLE 1 | Background characteristics of CLHIV on ART in Zambia from 2004 to 2018.

\**Tested using Cox proportional Hazards ratio.*

\*\**Adjusted for adjusted for age at ART start, sex and province.*

with valid data was performed to assess the mortality risks for those variables.

## RESULTS

A total of 70,718 children were recorded as having been initiated on ART from 496 health facilities and 71 districts across Zambia during the period under review. Of these only 65,448 were analyzed in survival analysis translating to 33,483 (51%) female, and 31,965 (49%) male.

The largest proportion of children analyzed were between 10 and 15 years and these accounted for about 25% of the sample while the smallest proportion was accounted for by children under 1 year of age at 10%. Mortality was also highest among children that initiated ART under 2 years, and lowest among children that initiated ART at older ages (10–15 years). Children that initiated ART during the period under review, were largely initiated on NRTIs (47%) and NVP based regimen (29%) (**Table 1**).

## Health Facility and Provincial Characteristics of CLHIV on ART

Out of the three levels of health facilities that were analyzed, hospitals and health centers accounted for the highest number of person years of follow-up (140,270 and 133,904, respectively) with health posts having the lowest person years of follow up (402). Out of the 3,265 children that died, hospitals in the country accounted for about 1,542 deaths translating to about 47% of deaths recorded during the reporting period. On the other hand, facilities in Lusaka (812, 24%), Copperbelt (583, 18%) and Southern (467, 14%) provinces had the highest number of deaths amongst all the ten provinces while facilities in Muchinga province had the lowest number of deaths with 81 deaths (2.5%). Children in Muchinga, Luapula and North western provinces also had the lowest chances of survival compared to children in Central province (North western, aHR 2.05 95% CI = 1.69–2.4; Muchinga, aHR1.95 95% CI = 1.50–2.50; Luapula, aHR 1.60 95% CI = 1.30–1.95) (**Table 2**).

## Survival of CLHIV on ART and Associated Factors

A total of 65,448 children were analyzed using proportional hazards survival analysis, with a total of 275,715 person-years at risk. The individual children had follow-up times that ranged from 1 month to 13.5 years. Survival of children on ART was highest among children who initiated ART between 10 and 15 years, with these children having 72% lower hazard ratio compared to children initiating ART between birth and 1 year (aHR 0.28,95% CI = 0.25–0.32). There was considerable variability in survival amongst the children initiating ART from the different provinces (**Table 2**). Factors associated with mortality included age when initiating ART and time on ART although these two variables were correlated as only children who initiated ART at younger ages could have a longer duration on ART. Compared to children who were diagnosed and initiated ART in the first year of life, children initiating ART between

\**Tested using Cox proportional Hazards ratio.*

\*\**Hazards adjusted for age at ART start and sex.*

\*\*\**Frailty model adjusting for clustering at facility level.*

*The variability of survival between the health facilities 2.18 (95% CI 1.86–2.63).*

the first and the second years of life were 32% less likely to die (aHR 0.68 95% CI = 0.61–0.75), with lower hazard ratios in each older age group (**Table 3**). Compared to WHO clinical stage 1, children with WHO clinical stages 3 (aHR 10.6; 95% CI = 9.70–11.74) and 4 (aHR 16.49 95% CI = 14.56–18.67) had increased mortality (**Table 3**). The frailty coefficient was >0, showing that the survival rates differed by a factor of 2.18 across the health facilities.

Three separate models were used restricting the analysis to subsets of children with data for CD4 counts at initiation, WAZ and the first line drug regime used at ART initiation. All models were adjusted for sex, WHO clinical stage, and time on ART, and the estimates for these variables were similar to the full model although the significance of the effects was greatly reduced (**Table 4**). Higher mortality was observed in those who were severely malnourished (aHR 4.51, 95% CI 2.57–7.92) compared to those who were not severely malnourished. Higher values of frailty for the model including malnourishment indicates the additional clustering of malnourishment in health facilities, increasing the variability of the survival rates by 3.22 across the health facilities. In the same model, those who initiated ART with a CD4 count of <500 cells/µl (aHR 1.76, 95% CI 1.22–2.54) had higher mortality than those who initiated ART with a CD4 count >500 cells/µl (**Table 4**).

## Estimation of Cumulative Hazards Curves of Death Among CLHIV on ART

Nelson Aalen cumulative Hazard estimates were analyzed for factors affecting survival in CLHIV on ART in Zambia and cumulative Hazard curves were derived. Factors analyzed included WHO staging, age at ART initiation, ART drug regimen and weight for age Z-scores. Overall cumulative hazards increased rapidly between the first and third year of follow up time across all factors that were estimated. The cumulative proportions of children who died on NVP based drug regimen was 15% higher than children who were on NRTI and EFV based regimen in the first 5 years of follow up time. Similar increases in mortality were also evident among children that initiated ART in WHO clinical stage 4, where steady increases of up to 30% were noticed for close to 10 years of follow up see **Figure 1**.

## DISCUSSION

Our study outlines survival of children on ART over a 13 year period (2004–2017) in Zambia. Studies in the past showed that without ART, more than 50% of CLHIV would die before their second birthday (5). For all children, mortality rates at ART initiation were high but lower than those reported in CLHIV prior to ART availability. Overall the SmartCare database showed 3,265 children on ART died during the 13 years, which gave an overall estimate of 1.1 deaths per 100-person years.

Mortality rates in our study decreased with longer duration on ART. Mortality rates in children who had been on ART for 5 years or more were only slightly higher than those reported for HIV negative children in Zambia (23). The overall mortality rate of 1.1 per 100 person-years reported in this study changed with duration on ART. In the first year on ART, the mortality rate was 12 per 100 person-years, which reduced to 8 per 100 person-years between the first and second year, 1.6 per 100 between the second and fifth year, and 0.62 after the fifth year. Our mortality rates over the first 2 years on ART were higher than those reported in a 2 year follow up study in Addis Ababa, Ethiopia (1.2 per 100 person-years), a 1 year follow up study in Makurdi, Nigeria (3.0 per 100 person-years), and a 2 year follow up study in Eastern Ethiopia (2.1 per 100 years). However, our overall mortality rate was lower than those reported in a 7 year South African study (4.7 per 100 person-years) and an eight-year study in Jos, Nigeria (14, 16, 17). The variation in mortality rates could be due to larger sample size in our study compared to other studies used in the comparison. Another reason could be the year in which the study was conducted as WHO guidelines being used at that time could have affected time to ART initiation, drug regimen and in turn mortality (10, 14).

Children that initiated ART in Lusaka, Copperbelt, and Southern provinces as well as children that initiated ART in hospitals and health centers had higher rates of mortality in our study. Male children on ART in our study had increased risk of mortality in contrast to findings by Zanoni et al. (15), in a South African study where female child had increased chances of mortality compared to their male counterparts. Nonetheless, our

TABLE 4 | Multi level adjusted and unadjusted hazards ratios for factors associated with mortality among CLHIV on ART in Zambia clustered at facility level.

*All models adjusted for sex and WHO clinical stage clustering at health facility level.*

findings are in line with findings by (16) where higher mortality rates was experienced among males (66.7%) than females (33.3%) as well as in a cohort of HIV uninfected children where male children had 2% higher mortality than female children sub Saharan Africa (24).

This study used a modified Cox regression to allow for differences in survival due to unobserved covariates at facility level. In these models, the estimated variance of health facility frailty terms was significantly >0, indicating large variability in the survival rates across facilities. The coefficient of variation for mortality across health facilities, showed mortality rates were more than twice as high (coefficient of variation = 2.18) in some facilities than the average. This difference could be due to varying resources in different facilities, training and availability of staff, and the quality of services. This high inequality needs to be addressed by the Ministry of Health to ensure equitable services to reduce HIV mortality for all Zambian CLHIV.

We also found that the risk of mortality was highest among infants and younger children initiated on ART, which is the period for highest risk of death for HIV negative children also. The risk of death reduced with increasing age. This could be due to non-adherence, sub optimal doses or exposure to perinatal antiretroviral therapy, which may lead to rapid disease progression and later death in infants and younger children (25). Another reason could also be that diagnosing HIV infection in infants is intricate and this can delay ART initiation in this age group (26). However, the proportion of children dying in this study did not account for over 50% of children analyzed in any age group showing improvements in survival compared to studies by Mulugeta et al. (5) and Tariro (6). Our study was in line with Studies in Asia and Africa that showed survival rates above 50% (84% to 97%) after 12 months of ART for children initiated on ART at a median age of 5–7 years, as well as another multi country study in Africa that also reported above 90% (93%) survival rates after 2 years of follow up on ART (18).

We found late initiation on to ART from time of HIV testing was predictive of reduced survival. Studies have shown that late presentation for early infant diagnosis (EID) services leads to poor prognosis (27, 28). This is also evidenced by findings of Bolton-Moore et al. (29) in earlier studies in Zambia. Our findings showed high mortality rates among children that initiated ART between 1 week and 1 month after HIV diagnosis. Current WHO guidelines call for same day ART initiation after HIV diagnosis which gives children the greatest benefit in survival (30). Our findings are in line with findings in a study done in South Africa that reported that 64% of deaths occurred within 3 months of ART initiation (15); a Nigerian study, where 81.3% and 84.4% of deaths occurred in the first 6 months and within 1 year, respectively (16) and an Ethiopian study that reported a high number of deaths within the first sixth months

of starting ART (31, 32). This pattern of increased mortality could suggest that HIV-infected children presented with poor immunological or clinical states at ART initiation or that CLHIV require some time before the benefits of ART are fully realized, especially when ART is initiated early (13, 18, 33). However, our study showed reducing mortality with increase in time to ART initiation above 1 months. This is similar to findings in a 5-year cohort study in Thailand where children that initiated ART in infancy had higher risks of mortality compared to children that initiated ART at older ages (18). The cohorts of children initiating ART at older ages, probably reflects those who survived with HIV, and therefore reflect a "survivor effect" of the older children (18).

The SmartCare database showed 3,265 children on ART died over the 13 years, but some deaths could have been missed and not reported to the health facility staff. In the SmartCare data, CD4 counts, weight and ART regime were not routinely recorded for every child but may only be recorded when triggered by a clinical event. We did not attempt to impute these data but instead carried out subgroup analyses in those children where the data are recorded to assess whether risk factors differed in the subgroups. The measures of advanced HIV diseases, low CD4 count, severe acute malnutrition, and WHO clinical stages 3 and 4 were all associated with reduced survival and this is seen in our results (13, 15). These results are consistent with earlier an study in Zambia which showed that mortality was higher in children with lower WAZ scores (29).

This study was not without limitations. The main purpose of the SmartCare data is for clinical care, treatment and follow up of CLHIV in health facilities where staff are overburdened. In these busy health care facilities vital information such as laboratory results and weight were not collected at every visit, resulting in a large number of missing data, which meant these factors could not be included in our main analysis The results for these variables should be interpreted with caution, as they did not represent all children. However, the consistency of the demographic risk factors across the subgroup analysis is encouraging. The SmartCare database may have missed some deaths when they did not occur at the health facility, as these might not have been reported, or entered into the data base, which could underestimate the mortality in children on ART. Moreover, loss to follow up of children initiated on ART is common among children on ART in Zambia and has been documented in earlier studies (20, 34).

## CONCLUSIONS

Children with HIV in Zambia are surviving much longer than was predicted before ART was introduced 14 years ago. This key finding adds to the literature on analysis of survival in CLHIV in low income settings like Zambia. However, survival was dependent on the age at which ART was initiated, the time on ART, type of health care facility used at ART initiation and the province where services were accessed. This highlights the need to increase investments in early infant diagnosis (EID) and make services more equitable across Zambia to ensure timely HIV testing and ART initiation for CLHIV.

## DATA AVAILABILITY STATEMENT

The Zambian Ministry of Health has the sole authority over the dataset that was used in this study, thus, it cannot be shared online. Any further information or interest to use the data should be addressed to the Zambian Ministry of Health (www.moh.gov.zm).

## REFERENCES

## AUTHOR CONTRIBUTIONS

TM, JT, and CM participated in the conception of the study, co-ordination, acquisition of data and drafted the manuscript. TM and JT carried out the statistical analysis. CM, PM, JT, and TM reviewed all the content related to interpretation of the findings and participated in the critical review and editing of the manuscript drafts for scientific merit and depth. All authors read and approved the final manuscript.

## FUNDING

This study was supported by the SEARCH (Sustainable Evaluation through Analysis of Routinely Collected HIV data) Project. We acknowledge funding by the Bill & Melinda Gates Foundation grant [OPP1084472] entitled 'Using routinely collected public facility data for program improvement in Tanzania, Malawi and Zambia'.

## ACKNOWLEDGMENTS

We would like to thank the Ministry of Health for granting access to the data that was analyzed. We would also like to thank the Bill & Melinda Gates Foundation for providing funding that made this study possible.

initiation in Southern Africa: the IeDEA-SA collaboration. PLoS ONE. (2013) 8:e81037. doi: 10.1371/journal.pone.0081037

urban Zambia: a retrospective cohort study. BMC Pediatrics. (2013) 10:54. doi: 10.1186/1471-2431-10-54

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a shared affiliation, though no other collaboration, with several of the authors JT and PM within the last 2 years, at time of review.

Copyright © 2020 Munthali, Michelo, Mee and Todd. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Misreporting of Patient Outcomes in the South African National HIV Treatment Database: Consequences for Programme Planning, Monitoring, and Evaluation

David Etoori 1 \*, Alison Wringe1 , Chodziwadziwa Whiteson Kabudula2 , Jenny Renju1,3 , Brian Rice4 , F. Xavier Gomez-Olive2 and Georges Reniers 1,2

*1 Department of Population Health, London School of Hygiene and Tropical Medicine, London, United Kingdom, 2 MRC/WITS Rural Public Health Transitions Research Unit (Agincourt), Faculty of Health Sciences, School of Public Health, University of Witwatersrand, Johannesburg, South Africa, 3 Kilimanjaro Christian Medical University College, Moshi, Tanzania, 4 MeSH Consortium, Department of Public Health Environments and Society, Faculty of Public Health and Policy, London School of Hygiene and Tropical Medicine, London, United Kingdom*

### Edited by:

*Marc Jean Struelens, European Centre for Disease Prevention and Control (ECDC), Sweden*

### Reviewed by:

*Philip Ayieko, University of London, United Kingdom Christen Ann Said, University of California, San Francisco, United States Deborah K. Glencross, University of the Witwatersrand, South Africa*

> \*Correspondence: *David Etoori<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Infectious Diseases – Surveillance, Prevention and Treatment, a section of the journal Frontiers in Public Health*

> Received: *25 April 2019* Accepted: *12 March 2020* Published: *07 April 2020*

### Citation:

*Etoori D, Wringe A, Kabudula CW, Renju J, Rice B, Gomez-Olive FX and Reniers G (2020) Misreporting of Patient Outcomes in the South African National HIV Treatment Database: Consequences for Programme Planning, Monitoring, and Evaluation. Front. Public Health 8:100. doi: 10.3389/fpubh.2020.00100* Background: Monitoring progress toward global treatment targets using HIV programme data in sub-Saharan Africa has proved challenging. Constraints in routine data collection and reporting can lead to biased estimates of treatment outcomes. In 2010, South Africa introduced an electronic patient monitoring system for HIV patient visits, TIER.Net. We compare treatment status and outcomes recorded in TIER.Net to outcomes ascertained through detailed record review and tracing in order to assess discrepancies and biases in retention and mortality rates.

Methods: The Agincourt Health and Demographic Surveillance System (HDSS) in north-eastern South Africa is served by eight public primary healthcare facilities. Since 2014, HIV patient visits are logged electronically at these clinics, with patient records individually linked to their HDSS record. These data were used to generate a list of patients >90 days late for their last scheduled clinic visit and deemed lost to follow-up (LTFU). Patient outcomes were ascertained through a review of the TIER.Net database, physical patient files, registers kept by two non-government organizations that assist with patient tracing, cross-referencing with the HDSS records and supplementary physical tracing. Descriptive statistics were used to compare patient outcomes reported in TIER.Net to their outcome ascertained in the study.

Results: Of 1,074 patients that were eligible for this analysis, TIER.Net classified 533 (49.6%) as LTFU, 80 (7.4%) as deceased, and 186 (17.3%) as transferred out. TIER.Net misclassified 36% of patient outcomes, overestimating LTFU and underestimating mortality and transfers out. TIER.Net missed 40% of deaths and 43% of transfers out. Patients categorized as LTFU in TIER.Net were more likely to be misclassified than patients classified as deceased or transferred out.

Discussion: Misclassification of patient outcomes in TIER.Net has consequences for programme forecasting, monitoring and evaluation. Undocumented transfers accounted for the majority of misclassification, suggesting that the transfer process between clinics should be improved for more accurate reporting of patient outcomes. Processes that lead to correct classification of patient status including patient tracing should be strengthened. Clinics could cross-check all available data sources before classifying patients as LTFU. Programme evaluators and modelers could consider using correction factors to improve estimates of outcomes from TIER.Net.

Keywords: HIV, retention in care, bias, South Africa, health information systems

## INTRODUCTION

At the end of 2017, it was estimated that there were 36.9 million people living with HIV (PLHIV) worldwide, with 70% of the disease burden situated in sub-Saharan Africa (1, 2). The World Health Organization's (WHO) revised HIV treatment guidelines in 2015 call for immediate provision of lifelong antiretroviral therapy (ART) to all people testing positive for HIV. By the end of 2017, 60% of PLHIV in sub-Saharan Africa were on ART (1, 2). Whilst ART initiation rates have been increasing over time, in order to reduce HIV transmission rates and achieve 90- 90-90 AIDS elimination goals, there is a need for accelerated increases in treatment adherence and retention in care (3–5). South Africa has the largest population of PLHIV worldwide, with an estimated 18.8% of the adult population aged 15–49 years old living with HIV, representing 7.2 million people (1, 2). By the end of 2018, an estimated 68% of PLHIV in South Africa were on ART (1, 2).

The rapid growth in access to ART has accentuated the need for an affordable and accurate way to monitor and evaluate treatment programmes (6–8), including documenting the number of people alive and on treatment, and programme impact on mortality. In the past, the progress of patients on ART was mainly monitored through patient cohorts (9) and tallying numbers of services rendered to inform resource allocation (8). However, evaluation of HIV programmes has proved challenging due to multiple data constraints. These include concerns about data reliability (8), and continued use of paper registers which often lack unique identifiers, suffer from incompleteness (10), and are cumbersome to use with increasing patient numbers and length of patient follow-up (11, 12). Another major concern is "silent transfers" whereby patients change clinics informally and without accompanying documentation, a phenomenon which has become more prevalent with the expansion of ART programmes (13, 14). As a result, there is concern that many high-burden countries are ill-equipped to report on the outcomes of patients in care and on treatment (6, 7, 15–17).

In order to address these concerns, many countries are scaling-up the use of electronic patient registers (11). However, challenges persist including insufficient linkages between clinics (10), insufficient training of staff who are responsible for entering this information (10), and staff shortages (18–21), resulting in some staff responsible for data management being stretched across multiple roles (22). This sometimes leads to poor workflow, and staff resistance which results in poor change management. Privacy and security issues (22) are also a major concern.

In 2010, South Africa adopted TIER.Net, a three-tiered monitoring approach involving paper registers (TIER 1 recommended for facilities with <500 patients), an offline electronic register (TIER 2—recommended for facilities with 500–2,000 patients) and networked electronic medical records (TIER 3—recommended for facilities with more than 2,000 patients) (11). This allowed for different tiers to be implemented in each facility based on the context and resources available at the time of implementation and typically involved a phased evolution, beginning with preparation for TIER.Net, installation and training, back capturing, live capturing and finally a live site able to produce monthly and quarterly reports with staff on-site to manage it. In 2014, an estimated 3,000 out of 4,000 public sector clinics in South Africa were using TIER.Net (23, 24) in one of the three phases of implementation. As of 2017, TIER 3 was still in its pilot phase (25).

ART patient outcomes have evolved since the start of national HIV treatment programmes. In several cohort studies of ART programmes in sub-Saharan Africa, there have been reports of higher rates of LTFU among patients who initiated ART in later years compared to earlier years (26–28). This may be explained by patients increasingly initiating treatment while less severely ill (29), as well as a negative consequence of patient numbers increasing such as facility workloads (30), raising concerns about the sustainability of these programmes. Some systematic reviews have shown that the percentage of patients LTFU who have died has decreased in later years as eligibility criteria have evolved to include less immunologically compromised patients, and as the proportion of patients LTFU has increased (13, 31). Furthermore, scale-up and decentralization of these programmes means ART can be offered at clinics closer to patients' homes, which may serve as an incentive to self-transfer in order to continue treatment at more convenient locations (13, 32).

Unpublished TIER.Net analyses from 2018 showed LTFU rates to range from 11 to 15% in the first three months and from 27 to 34% in the first year of ART (Y.Pillay, HIV Think-tank update, March 19, 2019). The high percentages of LTFU present many issues. Firstly, if these patients have really stopped ART then they have a higher mortality risk (33–35), and are more likely to transmit HIV (36–39). Given that patients that are LTFU have poorer outcomes, LTFU can also through misclassification bias event rates such as mortality downwards (40), leading to biased performance indicators for ART programmes. Accurate mortality rates are also important as they are used as parameters for projections such as in the UNAIDS spectrum package (41, 42).

We compare patient outcomes recorded in TIER.Net to the outcomes ascertained through a record review and tracing study for patients deemed lost to follow-up in eight public sector health facilities in rural north-eastern South Africa. We aimed to assess misreporting in TIER.Net and potential biases in the national programme statistics reported from the TIER.Net database.

## METHODS

## Setting

The Agincourt Health and Demographic Surveillance System (Agincourt HDSS) is located in rural north-eastern South Africa, in Mpumalanga province which has the second highest prevalence of HIV at 14.1% (43). HIV prevalence among people 15 years of age or older in the HDSS was estimated at 19.4% in 2010 (44). The Agincourt HDSS comprises of 31 villages covering an area of 475 square kilometers with an estimated population of 115,000 people (45, 46).

There are five primary health facilities and three secondary community health centers located within the Agincourt HDSS. Every HIV-positive patient has a clinical file that is opened when they first register at an ART clinic and updated at each clinic visit. Following the clinic visit, visit-level information from the patient file is entered into the national electronic database, TIER.Net. All health facilities routinely trace patients that are late for a scheduled clinic appointment. This tracing is done in conjunction with two non-profit organizations, Right to Care (RtC) and Home-Based Carers (HBC). Clinic staff must contact all patients first by phone and if this does not yield a satisfactory outcome, a home visit is organized. Patients are classified as lost to follow-up (LTFU) if they have not returned 90 days after their scheduled visit.

### Demographic Surveillance

Data collection aims to capture all demographic events for the Agincourt HDSS population. Fertility, mortality and migration data are based on a comprehensive household registration system that has been in operation since 1992. Following the baseline demographic surveillance survey in 1992 and three update rounds until 1998, the site has conducted annual surveys since 1999 (45, 47–49). Trained fieldworkers visit each household and interview the most knowledgeable adult available. During the visit, individual-level information on all household members is checked and updated and any events that have occurred since the last census round are recorded. Starting in 2017, data have been collected utilizing an electronic data collection system using tablets (50).

## Point-of-Contact Interactive Record Linkage (PIRL)

A key element of the data infrastructure for this study consists of HIV patient visit logs collected by a study fieldworker in the health facilities that provide ART in the area. This work started in April 2014 at seven government facilities and was extended in 2016 to include one additional health facility. In addition to logging patient visits, these records are linked to the Agincourt HDSS using a procedure that we have previously described as Point-of-Contact Interactive Record Linkage (PIRL) (51, 52). In brief, a fieldworker conducts a short uptake interview with patients in the waiting area of the clinic. Patients who consent are asked to declare a few personal identifiers that are used to search a local copy of the Agincourt HDSS database using a probabilistic algorithm. Matches are confirmed in interaction with the patients, and the names of other household members are used as a key attribute to adjudicate between possible matches.

## Record Review and Tracing Study

Through the PIRL database, we identified patients who were more than 90 days late for a scheduled clinic appointment from HIV services on August 15, 2017 (the date of data extraction) at any of the eight health facilities located in the Agincourt HDSS area. These patients were recruited into a cohort and followed up to ascertain their treatment and vital status i.e., whether they were still alive.

All PLHIV aged 18 years or over, who had ever declared residency in Agincourt HDSS, and had enrolled in the HIV treatment and care programme since 2014 (after the Agincourt HDSS record linkage was established at the health facilities) were eligible for inclusion in the study.

Trained and supervised fieldworkers conducted a thorough record review, comparing the list of patients LTFU against (a) TIER.Net (b) patient clinic files, and (c) logbooks kept by the RtC and the HBC. The PIRL database was also reviewed for duplicate patients (different clinical records linking to the same individual in the Agincourt HDSS database, which was taken as evidence of silent transfers), and residency and vital status were extracted from the Agincourt HDSS database. This was done on a case-by-case basis.

HBC conducted a visit to the households of all patients for whom a definitive outcome (defined as death, transfer out, stopped ART, migrated, re-engaged in care, and alive with ART status unknown) could not be established, or for whom routine patient tracing was not done. Finally, all patients who remained LTFU after the HBC visit, were searched for in TIER.Net databases at clinics in close proximity to their home residence to capture any further silent transfers.

We also reviewed the records for a stratified random sample of 162 patient records who were not LTFU as of August 15, 2017, in order to assess whether TIER.Net misclassified any patients that were still in care. This sample was chosen to include 18 patients from every clinic (six men, six non-pregnant women, and six women who initiated ART while pregnant) with the exception of one clinic which had recently merged with another, and from which we sampled 18 patients who had enrolled whilst in each of the clinics prior to the merger.

## Definitions

Definitions of terms used in this article are provided in **Table 1**.

## Statistical Analyses

For patients included in the record review and tracing study, we calculated counts and proportions for socio-demographic and baseline clinical characteristics, TIER.Net treatment status, the final outcome, and cross-tabulated TIER.Net status and the final outcome.

To assess the degree and direction of mis-reporting of patient outcomes in TIER.Net, we graphically present TIER.Net treatment status and the final outcome proportions by some selected patient characteristics. A Pearson's chi-square test was

### TABLE 1 | Definitions of terms used.

used to compare whether TIER.Net treatment status and the final outcome varied by all the categorical variables. We also present a cross-tabulation of patient outcomes from the two sources.

A binary outcome variable was created to indicate whether TIER.Net had misclassified a patient's treatment status, with a second outcome created to identify whether the patient was recorded as LTFU in TIER.Net. All cases where an electronic record could not be found were removed from further analysis.

To explore factors associated with misclassification in TIER.Net, we ran bi-variate analyses with patient-level treatment characteristics, demographic characteristics and facility-level characteristics. All variables with p < 0.1 were included in the multivariable logistic regression model. A parsimonious model was achieved using Wald tests. This same procedure was followed in order to understand what factors were associated with being reported LTFU in TIER.Net. All analyses were conducted in Stata 15 (53) and all data visualization was done using R (54).

### Ethics

Ethical approval was obtained from the London School of Hygiene and Tropical Medicine, the University of Witwatersrand and the Mpumalanga Department of health.

## RESULTS

## Database Population Characteristics

Over the study period, 4,089 patients were added to the PIRL database and met the inclusion criteria. Of these 4,089, 1,325 (32.4%) met the LTFU criteria and were eligible for inclusion into the record review and tracing study. Of these 1,325 patients, 166 (12.5%) did not have an ART initiation date and were assumed to be pre-ART. Further investigation of these 166 patients found 46 (27.7%) had initiated ART after record linkage, 59 (35.5%) were genuine pre-ART patients and 61 (36.7%) had initiated ART before record linkage began. These 61 patients were excluded from further analyses. Of the remaining 1,264 patients, 190 (15.0%) were found to have data errors (mostly due to missing clinic visits in the PIRL database) and were excluded from further analysis (**Supplementary Figure 1**).

Of 1,074 remaining patients, 280 (26.1%) initiated ART for prevention of mother-to-child transmission (PMTCT), 737 (68.6%) met the ART initiation criteria for non-pregnant adults, and 57 (5.3%) had not initiated ART yet (pre-ART).

Thirteen (8.0%) of the 162 patients still in care were excluded from the analysis because they had not declared residency in the HDSS. The remaining 149 from the random sample of patients still in care were also assessed to see if misclassification also occurred for those who remained engaged in care (**Table 2**).

## TIER.Net Treatment Status

Of the 1,074 patients who remained eligible for this analysis, 533 (49.6%) were categorized as LTFU, 222 (20.7%) as still in care, 186 (17.3%) as transferred out, 80 (7.5%) as deceased, and 53 (4.9%) could not be found in the TIER.Net database (**Table 2**).

There was a statistically significant difference (all p < 0.001) in the TIER.Net treatment status by sex, age, ART initiation status and reason, year of ART initiation, baseline CD4 count, time on ART, clinic visit schedule, health facility, and time since a missed appointment (see **Supplementary Figure 2**). Women who initiated ART for PMTCT were less likely to be categorized as deceased and more likely to be LTFU. All 149 patients sampled as still in care were also reported as still in care in TIER.Net.

## Outcomes After Record Review and Tracing Study

Of the 1,074 patients who remained eligible for this analysis, 326 (30.3%) were found to have transferred to another clinic, 234 (21.8%) to have re-engaged in care, 132 (12.3%) were deceased,

### TABLE 2 | Database population characteristics.

Etoori et al. Misreporting HIV Patient Treatment Outcomes

### TABLE 2 | Continued

*Not in sample: All patients eligible for the study but not LTFU in the PIRL database and not included in the still in care sample; Sample: All patients included in the study (149 still in care, 1,264 LTFU* = *1 074 really LTFU* + *190 data errors); Data error: Patients included as LTFU but found to be still in care and* <*90 days late for their last appointment; For ART start year data from 2017 reflects number of ART initiations up to mid-August when data extraction occurred.*

TABLE 3 | A cross-tabulation of TIER.Net treatment status and the final outcome.

117 (10.9%) were alive with ART status unknown, 81 (7.5%) were alive but not on treatment, 53 (4.9%) had migrated to another place of residence, and 131 (12.2%) were still LTFU (**Table 3**). These outcomes differed (all p < 0.001) by sex, age, ART initiation status and reason, baseline CD4 count, time on ART, clinic visit schedule, health facility, whether the patient record was successfully linked to an Agincourt HDSS record,

and time since a missed appointment (some selected variables illustrated in **Figure 1**).

## Differences Between TIER.Net Treatment Status and Final Outcomes

Records of deceased or transferred out patients documented in TIER.Net aligned with patients' final outcome (i.e., no inaccuracies found for these two statuses). However, TIER.Net misclassified 52 (39.4%) of 132 deaths. Of these 52, 38 (73.1%) were classified as LTFU, 6 (11.5%) as still in care, and 8 (15.4%) were not found in the system at all.

TIER.Net also misclassified 53 patients as still in care. Of these, 10 (18.9%) were found to be LTFU, 16 (30.2%) to have transferred, 12 (22.6%) as alive with unknown ART status, 8 (15.1%) alive but not on treatment, 6 (11.3%) to have died, and 1 (1.9%) to have migrated to another place of residence. TIER.Net correctly captured 186 (57.1%) of 326 transfers.

Of 533 patients classified as LTFU by TIER.Net, 116 (21.8%) were found to have transferred to another clinic, 70 (13.1%) to be alive but not on treatment, 47 (8.8%) to have migrated to another place of residence, 38 (7.1%) to have died, and 56 (10.5%) to have re-engaged in care (38 of whom were resolved by new visit data in the PIRL database and so it is possible that their TIER.Net status could have also changed back to still in care) (**Figure 2** and **Supplementary Figure 3**).

As patients classified as LTFU in TIER.Net were more likely to be misclassified we report on the factors associated with being classified as LTFU in TIER.Net (**Supplementary Data** and **Supplementary Table 1**).

## Factors Associated With Misclassification

In the multivariable model (**Table 4**), men (OR: 1.47, p = 0.021) had higher odds of misclassification when compared to women who initiated ART for non-PMTCT reasons (CD4, WHO stage, tuberculosis coinfection). Higher baseline CD4 (CD4 100- 199 OR: 1.95. p = 0.002, CD4 ≥ 500 OR: 1.81, p = 0.014) was also associated with higher odds of misclassification when compared to patients who initiated treatment with CD4 < 100. Health facility also remained statistically significant suggesting that facility level variability plays a role in misclassification. Patients who were linked to an Agincourt HDSS record in the PIRL database (OR: 2.09, p < 0.001) were more likely to be misclassified. Finally, patients who were between 1 and 2 years late (OR: 1.62, p = 0.001) were more likely to be misclassified. Older age (30–44 years OR: 0.73, p = 0.046, 45–59 years OR: 0.63, p = 0.046) was associated with lower odds of misclassification and patients on longer refill schedules (>3 months OR: 0.31, p = 0.009) were less likely to be misclassified.

## DISCUSSION

In this paper, we described the discrepancies between the treatment, vital and residency status of HIV patients enrolled in care between April 2014 and August 2017 in a rural South African setting as recorded in the national treatment database (TIER.Net) and their treatment outcome following a comprehensive record review and tracing study.

We found that TIER.Net misclassified 36% of the patient outcomes. ART initiation reason, baseline CD4, health facility attended, PIRL linkage, time since the last appointment, age, and ART refill schedule were all found to be significantly associated with misclassification. TIER.Net underestimated mortality and overestimated the number of patients who were LTFU. Seventynine percent of patients classified as LTFU in TIER.Net had a final outcome ascertained, mirroring findings from a systematic review of low and middle income country ART programmes which found that tracing generated higher estimates of mortality and lower estimates of LTFU (55). Our findings show that LTFU is still an important problem in ART programmes in this setting, even with routine patient tracing in place. TIER.Net also missed 43% of transfers with these silent transfers being the biggest contributor to misclassification among those documented as LTFU. We also found that 21.8% of patients had re-engaged in care, a phenomenon that was previously not well understood, but which is now increasingly recognized as becoming a common feature of ART programmes (55). Using our findings to

revise LTFU figures to reflect re-engagement may help improve programme evaluation and forecasting.

In our study, we found that 40% of deaths were missed by TIER.Net, indicating that mortality of ART patients would be underestimated if relying on this data source. Given the role that national statistics play in HIV/AIDS projections (41, 42, 56), our findings suggest a need for correction factors for the estimates of the effect of ART on mortality. Although South Africa has a good vital registration system in place (57), these data are not currently linked to clinic-based information. However, with the move to registering patient national ID numbers, clinics should consider matching patients that are LTFU to the national death registry and other available databases such as the national health laboratory services database to ascertain vital status as this has proved useful in other studies in South Africa (58–60). Clinics in the Agincourt HDSS study area and other HDSS sites could also consider using vital status data from annual demographic surveillance to ascertain vital status for all patients.

The number of patient transfers to another clinic that were missed in TIER.Net suggests that communication between clinics is sub-optimal and that the current system for transferring patients between clinics can be improved. With studies reporting patient fear and concern about provider reactions if they return to care after a treatment interruption (61, 62) it is possible that some patients considered it less stressful to self-transfer or restart treatment at another nearby clinic, rather than returning to the facility where they had initiated treatment. These silent transfers could lead to double counting of patients currently on treatment, the second of the 90-90-90 targets, potentially suggesting that the programme is performing better than it is. Furthermore, given that the national treatment programme relies on data from TIER.Net to plan and procure ART based on active patient numbers, misclassification in the database, and more specifically double-counting due to silent transfers may lead to inaccurate drug forecasts and misestimation of medicines and other commodities at the national level. This bias will only increase as the South African ART programme expands with more patients potentially moving into new clinics closer to their homes and more people initiate ART with the move to test and treat. Future work will consider how application of correction factors from this research would change programme statistics and drug forecasts.

It is also important to consider the risk that silent transfers pose with regards to drug resistance, as this misclassifies treatment experienced patients as treatment naïve and may lead to patients being offered regimens that have lost their optimal therapeutic benefit. This is particularly concerning because resistance testing is not commonly used in these settings, and can potentially lead to increases in levels of transmitted drug resistance (63). Better referral systems, patient education, regular information exchange between clinics, and provider training (64), could improve recording of transfers and clinic staff attitudes toward less adherent patients. The WHO also recommends enforcement of unique identifiers as paramount to improve patient safety, improve the efficient use of programme resources by reducing duplications, and to improve programme monitoring and evaluation (65). With national IDs becoming mandatory at clinic registration, information exchange could prove useful in identifying silent transfers. This should also become less of an issue when TIER.Net is upgraded to a fully networked database.

We found several factors to be associated with misclassification of outcome in TIER.Net, with older age and longer ART refill schedules found to be protective factors. Older patients were less likely to be classified as LTFU in TIER.Net which probably explains why they were subsequently

less likely to be misclassified. Given that longer ART refill schedules are synonymous with previous good adherence (66), these patients accounted for 11% of the patients LTFU and were also more likely to have re-engaged in care, a category that contributed very little to misclassification. Patients whose clinic record was successfully linked to an Agincourt HDSS record in the PIRL database were more likely to be misclassified. They were also more likely to have been resolved which could explain this association. Health facilities were also positively associated with misclassification, with the facilities with the highest proportion of patients classified as LTFU in TIER.Net being more likely to misclassify patients. Two of these clinics also had issues with routine tracing, with one clinic not undertaking any physical tracing at all, emphasizing an additional benefit of routine tracing. Finally, patients who had been LTFU for a longer duration were more likely to be classified as LTFU in TIER.Net, more likely to have transferred to another clinic, and less likely to have re-engaged in care which probably explains their higher likelihood of misclassification.

This analysis has several limitations. Firstly, TIER.Net was only consulted at a specific point in time. The cross-sectional nature of TIER.Net outcomes means that some may have changed, but we would have no way to ascertain this. We checked TIER.Net 12 months after the initial record review for all patients whose outcome after record review and tracing was still LTFU and 85% of the outcomes had not changed. However, for patients whose final outcome was resolved through new visit data in the PIRL database, it is likely that their TIER.Net outcome also changed. It is also possible that some of the patients categorized as LTFU in TIER.Net are due to the rigidity of the system as TIER.Net only allows for four possible outcomes; still in care, transferred out, LTFU and deceased (11, 67, 68). It is possible that for some patients, their outcomes were ascertained, but the rigidity meant that they could not be recorded in the database and may call for the inclusion of other possible outcomes in the database. The exclusion of patients for whom an electronic record could not be found from the multivariable analyses might bias our findings. However, given the relatively small number we expect that this bias is fairly small. Finally, we did not adjudicate causes of death, so it is possible that patients died from causes other than those related to HIV/AIDS. A strength of this study is that we attempted to trace all patients that were LTFU and not a sample. Therefore, the findings might be more generalizable to other settings. The multiple methods, data sources and levels of follow-up used to trace patients are also a strength.

In conclusion, although TIER.Net misclassified 36% of patient outcomes, this reflects the various challenges with the processes and upstream factors that lead to this misclassification and calls for their improvement rather than the utility of the database itself, as patients classified as LTFU were most likely to be misclassified. Clinics should consider training staff about ascertaining patient outcomes, putting more emphasis into patient tracing and using other data sources such as the national death register to improve ascertainment of patient treatment outcomes. For policy and planning purposes, programme evaluators should consider using correction factors to improve the accuracy of estimates from TIER.Net.

## DATA AVAILABILITY STATEMENT

The data used in these analyses are not yet publicly available as they currently being utilized for the first author's Ph.D. research. They will be deposited with the Agincourt HDSS data manager and made available on request at the end of his Ph.D. in 2021. Data from the PIRL database are also available by making a data request to the Agincourt HDSS data manager.

## ETHICS STATEMENT

Ethical approval was obtained from the London School of Hygiene and Tropical Medicine, the University of Witwatersrand and the Mpumalanga Department of health.

## AUTHOR CONTRIBUTIONS

The study was conceived by DE, AW, and GR. Fieldwork was planned and executed by DE, FG-O, and CK. Data collection was supervised by DE. Analyses were conducted by DE with input from all authors. The manuscript was drafted by DE with input

## REFERENCES

from BR, JR and all the other authors. All authors contributed to the interpretation of the findings and read and approved the final manuscript.

## FUNDING

This study was made possible with support from the Economic and Social Research Council (ES/JS00021/1), the Bill and Melinda Gates Foundation for the MeSH Consortium (OPP1120138), the Bill and Melinda Gates Foundation ALPHA grant (OPP1164897), and the MRC SHAPE UTT grant (MR/P014313/1).

## ACKNOWLEDGMENTS

The authors would like to thank all the participants in the study.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpubh. 2020.00100/full#supplementary-material

Sub-saharan Africa. J Acquir Immune Defic Syndr. (2017) 75(Suppl 2):S115– 22. doi: 10.1097/QAI.0000000000001343

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer PA declared a shared affiliation, with no collaboration, with the authors to the handling editor at the time of review.

Copyright © 2020 Etoori, Wringe, Kabudula, Renju, Rice, Gomez-Olive and Reniers. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.