Analysis of Complex Correlated Interval-Censored HIV Data from Population Based Survey ()
1. Introduction
Interval-censored data arise in research settings where the exact time an event occurs is not directly observed. In epidemiological studies that are cross-sectional in design, the time of the occurrence of an event is not observed but only an indication that it has either occurred or not. Knowing exact time of the occurrence of an event is crucial when interest is in estimating incidence of diseases such as HIV. HIV prevalence data is the most commonly available surveillance data. However, HIV incidence provides a more useful measure of epidemiological trends and sensitive indicator for evaluating impact of interventions. Incidence data are rarely available because their collection requires difficult, time-consuming and expensive follow-up of large cohorts. Thus, alternative laboratory based assays that can distinguish recent from long-term HIV infections from changes in antibody characteristics after seroconversion are becoming popular [1] . One such method is the BED capture enzyme immunoassay technique which can detect if HIV positive blood is infected within the last six months or longer [2] and thus leading to interval-censored data for those that are infected within the last six months, left censored data for those that are infected longer than six months and right censored data for those that are HIV negative at the time of the study. Estimation techniques for interval-censored data often stem from Cox proportional hazards model [3] . The model has been generalized to handle interval-censored data [4] .
In multi-stage cluster sampling, individuals belonging to the same cluster such as an EA share the same unobserved cluster-specific frailty effect and thus making them positively correlated. For example, in the HIV prevalence survey conducted in South Africa as described in Reference [5] , individuals from the same EA are considered to be more similar than those that are from different EA with respect to their observed risk behaviour, risk of HIV infection and other unobserved factors. These data are the subject of our analysis in Section 2. Frailties are considered mutually as independent random variables specified by some parametric distribution. Frailties are assumed independent given the parameters of their marginal distribution. The topic of frailty models has received considerable attention in demography [6] and statistics [7] in recent years. Gamma distribution for frailties, commonly assumed in frailty models, is conjugate to the right-censored data likelihood and thus simplifies computations. Several authors in both frequentist [8] [9] and the Bayesian viewpoints [10] - [12] have used gamma frailty.
The inclusion of frailty effects in the interval-censored data likelihoods proposed in references [4] [13] results in complex and intractable likelihood functions that often require numerical integration methods. The gamma frailty distribution is no longer conjugate. In the univariate context, the common fix-up approach is to assume that the event occurred at the beginning or end of each examination time and then use standard estimation techniques. However, this approach has been shown to lead to misleading results [14] [15] . The alternative approach is to ignore correlation in the data and utilize standard univariate techniques for interval-censored data [4] [13] . But, naive standard errors obtained through this approach can lead to invalid statistical inference [16] [17] . Other approaches that consider both infection time and frailties as unobserved random variables thus facilitating the use of expectation- maximization algorithm have been proposed [12] [18] .
In cross-sectional prevalence survey, the numbers of recently infected, and HIV negative can be used to estimate HIV incidence [19] . However, these estimation methods only provide point estimates of incidence with no functional forms for modeling the data and inclusion of key determinants of incidence. In this paper, we present an approach for joint modeling of left, right and interval-censored data when event times are correlated. The model is implemented using a population based household HIV prevalence survey that used BED techniques to classify infection as either recent or latent to estimate HIV incidence by different reporting domains [19] .
2. The Data
The data is from a household based population survey conducted in 2005 in South Africa. The design of the survey was based on a multistage stratified, clustered sampling scheme. The survey is conducted every after three years since 2002 to measure, partly, the impact of HIV. The survey targeted all persons that are 2 years old and above living in households. The survey excluded individuals living in educational institutions, old-age homes, hospitals and uniformed services barracks but included those living in hostels. The survey applied a multi-stage stratified sampling approach based on a master sample of 1 000 EA used by Statistics South Africa for the 2001 census. An EA consists of about 250 households. In this survey, 15 households were systematically selected within each EA. Three persons in each household were potentially eligible to be selected for the survey. And only one person was selected from each of the following age groups: 2 to 14 years old, 15 to 24 years old and 25 years old and above.
A total of 23,275 individuals aged 2 years old and older from 10,584 households were interviewed in the 2005 survey and 15,851 respondents agreed to be tested for HIV. Linked anonymous HIV testing was performed using dried blood spot (DBS) specimens. Socio-demographic and behavioural data were collected using a detailed questionnaire. Among those that were tested, 1449 tested HIV positive and among those that were HIV positive, 181 were diagnosed as recent infections using the BED [1] . This meant that 1268 (8%) were considered left censored, 14,402 (90.9%) right censored and 181 (1.1%) interval-censored. Point estimates of incidence from these data have been reported elsewhere [19] . The average number of people who were tested in an EA were 22 (standard deviation = 9.3). Covariates in our analysis include sex, race, age, locality type and marital status. Summary statistics for these variables are presented in Table 1.
3. Model and Parameter Estimation
The data to be analysed in this paper is clustered within an EA. Let 
denote the event time for the jth individual in cluster i. Due to interval-censored nature of the data, the exact event time 
is unobserved. Instead, we only observe clinical endpoints 
encompassing interval-censored event time 
for those that are considered recently infected (infected within the last six months), 
and 
for those that are left and right censored respectively. Define censoring indicator 
and 
if left censored; 
and 
if interval-censored and 
and 
if right censored. This is similar to defining 
if right censored. The ith cluster specific frailty is denoted by
. The frailties are assumed to operate multiplicatively on the hazard function. Furthermore, Weibull baseline hazards 
are assumed. The corresponding integrated baseline hazards are![]()
.
The model is framed as a shared frailty model defined as
![]()
where ![]()
is the conditional hazard function of the jth subject from the ith EA. Conditional on unobserved
![]()
Table 1. Summary statistics for variables used in the hazard model.
frailty![]()
, observations from cluster i are independent. The contribution from cluster i can thus be expressed as
![]()
where ![]()
is the assumed density function for unobserved frailties. The conditional survival distribution
is ![]()
where ![]()
is the integrated hazards corresponding to
![]()
, ![]()
and ![]()
is the vector of covariates and β represents the corresponding covariate effect. The likelihood ![]()
is an extension of a classical case of general interval censoring described by [13] with
![]()
where ![]()
is the joint density function of ![]()
In this model, the density function for unobserved frailties is assumed to follow a Normal ![]()
distribution. The likelihood has no closed form solution and thus numerical integration methods will be used before maximising the marginal likelihood function. Since ![]()
is assumed to be normally distributed then Gaussian quadratures will be used for numerical integration. The full data likelihood is then approximated by
![]()
Maximum likelihood estimates are calculated using Newton-Raphson optimization for the corresponding (approximate) full data log-likelihood function.
A similar model with parameterization consistent with the accelerated failure time model formulation used in PROC LIFEREG in SAS has been used before [20] . Here, we have formulated the likelihood in the proportional hazards settings by adding a frailty term to the (linear) log-hazard component of the model. PROC NLMIXED in SAS can be used to find maximum likelihood estimates for this model. NLMIXED uses adaptive Gaussian quadrature to approximate the cluster-specific likelihood contribution. This is an attractive feature of this approach that it can be readily implemented using existing commercial statistical computing software.
4. Results
The methods of estimating point estimate of incidence described by [1] uses an adjustment formulae that corrects for both sensitivity and specificity specified by the BED analysis. Using this method, the incidence of HIV among persons aged 2 years old and older is estimated at 1.4% per 100 person-years [95% CI: 1.0 - 1.18] with 571,000 new HIV infections estimated for 2005 in South Africa as first reported in [19] . Incidence of HIV appears to vary by age, sex, race, geographical type of the area and by other socio-behavioural factors. Being female or aged between 15 to 24 years old or being an African is associated with high incidence of HIV. Those residing in informal settlements have the highest incidence of HIV compared with other residential areas. Similar results are observed for those that are never married.
Table 2 presents univariate and multivariate results of the model presented in Section 3. In univariate analysis nearly all of the covariates are statistically significant, and the estimated effects of the covariates are in close agreement with findings from previous research. The hazards of being infected with HIV is significantly lower for males than females [HR = 0.72, 95% CI: 0.64 - 0.81]. The hazards of HIV infection are [HR = 1.68, 95% CI: 1.34 - 2.11] times higher for youth aged 15 to 24 years old compared with those aged 2 to 14 years old and similar to those aged 25 years old and above [HR = 1.57, 95% CI: 1.22 - 2.01]. Africans are associated with significantly higher risk of HIV infection than any other race [HR = 7.43, 95% CI: 6.17 - 8.95]. Urban informal areas have the highest risk of HIV infection compared with other residential areas [HR = 2.38, 95% CI: 2.24 - 2.54], rural formal [HR = 1.56, 95% CI: 1.30 - 1.87] and rural informal [HR = 1.94, 95% CI: 1.71 - 2.20] in relation to urban formal. Being married or ever been married has a protective effect against being infected with HIV. Those who have never been married are [HR = 2.91, 95% CI: 2.57 - 3.30] times more likely to be infected with HIV compared to those who have ever been married or are cohabiting.
![]()
Table 2. Parameter estimates of univariate and multivariate analyses.
Table 2 further presents results from a multivariate frailty models. The model with frailty term contains the effects of other covariates not specifically included in the model, which are common to individuals from the same EA. For example, the clustering effects at an EA may result from common community norms, exposure to common underlying risk factors of HIV such as other sexually transmitted infections that are endemic in some geographical areas. The frailty parameter represents the EA effect. The parameter is interpreted as the variance of the random effects distribution. Large values of the variance indicate greater heterogeneity between EAs and stronger association within members from the same EA. The estimated variance of the EA is 0.4947 (SD = 0.0437), highly significant when compared to the one-sided Z-test with critical value of 1.645 at 5 per cent level of significance. The estimate of the variance implies moderate degree of correlation in the risk of HIV infection between members from the same EA even after controlling for observed covariate effects. The estimated correlation is 0.198. Standard errors are slightly magnified in the model with frailty indicating that the fixed effect parameters are now estimated more realistically, but with lower precision. Therefore, treating each individual response as independent gives false impression that there is more information in the data than there really is.
The hazard estimates for locality type are drastically reduced in the multivariate frailty model. The inclusion of EA random effect term accounts for unobserved proneness of individuals within specific EAs which is a locality type in its own right. Therefore, the magnitude of the effect of locality type is reduced in the model that includes EA effect. This is important especially for intervention programs as it indicates the need to focus interventions more in some locality types than others. The lack of striking differences between the model with and without frailty is reassuring since the effect of EA is rarely considered in studies of HIV incidence.
5. Conclusions
Analyses conducted in this paper provide an approach to a complex problem of correlated interval-censored data, where the event time is not observed precisely. The model demonstrates the possibility of multivariate analyses of HIV incident data and also includes frailty effect in this complex model. Inclusion of frailties in the analyses of correlated HIV infection risk at an EA level is of critical programatic value. Evaluation of the impact of HIV intervention programs and their planning relies heavily on incident data. The importance of clustering effect at an enumerator level demonstrates the importance of focusing interventions at key geographical areas. HIV in South Africa is at an epidemic level with considerable geographical variability in HIV distribution with some areas carrying much higher burden at local levels than others [21] . Even in high HIV prevalence areas, there is considerable evidence of localised epidemics [22] . Analyses of data that adjusts for clustering effect provides an opportunity to investigate factors and identify potential geographical areas that contribute disproportionately to the epidemic of HIV. Even after controlling for some important risk factors, the risk of HIV infection varies considerably across EA. The importance of EA frailty term indicates that some geographical areas are at increased risk of HIV infection compared with others. Control measures of HIV should extend further from high-risk individuals to high-risk areas.
In a study conducted at Africa center [23] the hazard of HIV seroconversion was approximately twice higher in people who were currently unmarried but had a partner than among people who were currently married. In our univariate analyses the hazard of HIV seroconversion was approximately three times higher in unmarried. However, in the multivariate model, marital status was excluded due to collinearity with age. Baseline hazard parameters indicated a higher risk of HIV infection at young ages that gradually subsided for older people.
In the estimation phase, Gaussian-Hermite quadratures are used to approximate intractable integrals [24] . The final results are based on six quadrature points. Brillinger and Preisler [25] reported that results did not change much for quadrature points greater than eight. Moreover, even five quadrature points have been considered sufficient [26] [27] . The inherent disadvantage of using large number of quadrature points is the need for strong assumption of normally distributed random effects. For example, ten quadrature points fit a symmetric distribution thus forcing the tails of random effects distribution to be normal.
In conclusion, the results provide further insight to the understanding of the epidemic of HIV. The models fitted are quite complex. However, the ability of fitting these complex models using the standard software is encouraging.
Acknowledgements
Authors acknowledge generous funding received from Nelson Mandela Foundation, Swiss Agency for Development and Cooperation, USA Centers for Disease Control and Prevention and the Human Sciences Research Council. We wish to thank our colleagues and collaborators who participate in the survey. Finally, we wish to sincerely thank all the people of South Africa who willingly open their doors and give researchers some of the most private information and agree to also provide specimen for HIV testing.