Parametric frailty models for clustered data with arbitrary censoring: application to effect of male circumcision on HPV clearance
© Kong et al; licensee BioMed Central Ltd. 2010
Received: 4 February 2010
Accepted: 6 May 2010
Published: 6 May 2010
In epidemiological studies, subjects are often followed for a period during which study outcomes are measured at selected time points, such as by diagnostic testing performed on biological samples collected at each visit. Although test results may indicate the presence or absence of a disease or condition, they cannot provide information on when exactly it occurred. Such study designs generate arbitrarily censored time-to-event data, which can include left, interval and right censoring. Adding to this complexity, the data may be clustered such that observations within the same cluster are not independent, such as time to recovery of an infectious disease of family or community members. This data structure is observed when evaluating circumcision's effect on clearance of penile high risk human papillomavirus (HR-HPV) infections using data collected from the male circumcision(MC) trial conducted in Rakai, Uganda, where the multiple infections within individual and HPV testings performed at trial follow-up visits gave rise to the clustered data with arbitrary censoring.
We describe the use of parametric proportional hazards frailty models and accelerated failure time frailty models to examine the relationship between explanatory variables and the survival outcomes that are subject to arbitrary censoring, while accounting for the correlation within clusters. Standard software such as SAS can be used for parameter estimation.
Circumcision's effect on HPV infection was a secondary end point in the Rakai MC trial, and HPV genotyping was conducted for penile samples of a subset of trial participants collected at enrollment, 6, 12 and 24-month follow up visits. At enrollment, 36.7% intervention arm men (immediate circumcision) and 36.6% control arm men (delayed circumcision at 2 years) were infected with HR-HPV, with the number of infections per man being 1-5. The proposed models were used to examine whether MC facilitated clearance of the prevalent infections. Results show that clearance of multiple infections within each man is highly correlated, and clearance was 60% faster if a man was circumcised.
Parametric frailty models provide viable ways to study the relationship between exposure variables and clustered survival outcome that is subject to arbitrary censoring, as is often observed in HPV epidemiology studies.
In epidemiological studies, subjects are often followed over time and study outcomes are measured at selected time points. The study outcomes may be diagnostic testing results based on biological samples such as tissue, blood, or urine samples. The testing at each time point can detect the presence or absence of a condition (for example, infection of an infectious disease), but it does not provide the exact information of when the infection or condition occurred. The best knowledge about the actual event time is that it occurred during the interval where discordant test results are observed at the start and end of the interval, yielding so-called interval censoring of time-to-event (survival) data . If the event is observed at the first follow-up, we only know that the event occurred before the first scheduled testing time, and this generates left censored data. On the other hand, if a subject drops out of the study or remains event free at the end of the study, the time to event could only be after the last observed testing time, which corresponds to the well-known right censoring. When it is of interest to evaluate the effect of treatment on time to event occurrence, many analysts use either the end point  or the mid-point of the interval where discordant results were observed as the actual event times. The former way introduces "time aggregation" bias when estimating the hazard rate, while the latter mid-point method reduces time-aggregation bias under certain conditions . For unbiased estimation of the treatment or exposure effect, approaches that directly model the likelihood of the arbitrarily censored data (including left, interval and right censoring) can be used . Other than the presence of arbitrarily censored data, what can further complicate the analysis is the presence of clustered data where the study subjects are correlated within each cluster, such as patients visiting the same clinic in a multi-site study, members within the same family or community, or repeated testing on the same subject. The Cox regression model has been extended for clustered time-to-event data to model the marginal distributions without full specification of the correlation structure between the clustered observations [5, 6], though it has only been established for data with non-informative right censoring. For clustered data with interval censoring,  introduces the use of conditional proportional hazards model (i.e. semiparametric frailty model), and briefly discusses the advantages of parametric frailty models. In this paper we describe the use of parametric frailty models to assess treatment or exposure effect on time to a well-defined event for clustered survival data with arbitrary censoring, which includes left, right, as well as interval censoring. The model estimation can be carried out using existing software, such as SAS (SAS Institute, Inc., Cary, NC) PROC NLMIXED. The method is illustrated through the application on data collected from a randomized clinical trial of male circumcision (MC) conducted in Rakai, Uganda, and the current study purpose is to evaluate the effect of MC on clearance of penile high risk human papillomavirus (HR-HPV) in HIV-negative men.
The following section provides a brief introduction about the HR-HPV data to exemplify the problem of interest. Statistical notations and the proposed model are then presented through the context of the example. Parameter estimation can be realized using SAS and the SAS code for this example is provided. The analytical result on the HR-HPV data is subsequently presented. We conclude with discussion about the general use of the proposed model for clustered survival data with arbitrary censoring.
Example: Effect of Male Circumcision on HR-HPV clearance
Enrollment prevalence of HR-HPV infections and observed clearance proportions in the intervention arm (I) and the control arm (C).
n = 330 (%)
n = 314 (%)
Rate Ratio (I/C)
No. of men infected with Any HR-HPV
Total No. of infections No. of infections cleared
Clustered data structure arises in this dataset, as each participant (a cluster) had testing results for the fourteen HR-HPV genotypes. A person with multiple infections may contribute multiple events (i.e. clearances) that are likely to be highly correlated. The exact clearance time point is unavailable, but it is known to be either before the first subsequent visit (left censored), or lie within an interval between two visits (interval censored), or after the last follow-up visit (right-censored).
Let T denote the random variable for the time to event. Assume there are n clusters in the study, and the ith cluster (i = 1, 2, ⋯, n) has j = 1, 2, ⋯, n i observations. In the HR-HPV example, T is the time to clearance of an HR-HPV infection, and each participant is a cluster. For a person with single genotype HR-HPV infection, n i = 1; while if a person has multiple HR-HPV genotype infections, n i > 1, and the clearance of these multiple infections is not independent. For the jth HR-HPV genotype infection (j = 1, 2, ⋯, n i ) on the ith participant, t ij is the actual clearance time which is not exactly observed. Let (a ij , b ij ] denote the interval where the clearance occurs, i.e. a ij <t ij ≤ b ij . If a prevalent infection of genotype j for person i is observed to be cleared at first follow-up, then a ij = 0, and t ij ≤ b ij = 6 months, which corresponds to a left censored observation; if it is never cleared during the trial period, then a ij = 24 months (the last follow-up time), and a ij <t ij < ∞ (the upper bound b ij can be viewed as ∞), which is a right censored observation.
Suppose there are p explanatory variables, which can include the variable indicating treatment arm and other covariates that may be of interest, and let x ij denote the vector of explanatory variables for the jth infection in person i. Without loss of generality, we only consider one explanatory variable for treatment arm, and x ij = 1 denotes intervention and x ij = 0 denotes control. Note that treatment arm for the fourteen HR-HPV genotypes is identical for each participant, although in general applications, x ij can be different for different observations belonging to the same cluster.
Parametric Proportional Hazards Frailty Model
Without the ξ i term, expression 1 is just the ordinary Cox proportional hazards model, where h ij (t; x) is the "hazard" for clearance of the jth genotype infection in person i, h 0(t) is the baseline hazard for a person with all explanatory variables being zero. The ξ i term in expression 1 represents a random effect realized on the ith person. It is assumed to follow a prior distribution, such as a normal distribution, or equivalently exp(ξ i ) ~ a log-normal distribution. The use of the normal random effect ξ i on t (through the hazard function) is one way of introducing correlation within the ith cluster, and is similar to that in linear mixed effects or generalized linear mixed effects models. Because the same realization value of ξ i (from its prior distribution) is shared by observations on the multiple HR-HPV genotype infections within the ith person (thus its' subscript does not depend on j which indexes genotype), it is therefore taken into account of the dependence between the clearance times of these multiple HR-HPV infections. One advantage of using normally distributed random effect is that more complicated correlation structures between observations, such as multi-level correlations, can be handled naturally by extending the use of a univariate normal random effect to multivariate normal random effects [11, 12].
where the random frailty effect is assumed to follow a normal distribution with zero mean (i.e. exp (ξ i ) ~log-normal distribution). Model 2 implies that conditional on the random frailty ξ i , the clearance of any HR-HPV infection follows a Weibull distribution with shape parameter γ and scale parameter , where circumcision's effect on clearance is quantified by its' coefficient parameter β. ξ i is the random effect shared by the multiple infections within the same person. Given a value of ξ i , the hazard ratio of HR-HPV clearance between the intervention and control arm is exp (β). Therefore, if β estimated from the data is significantly larger than 0, then exp(β) > 1, indicating the instantaneous clearance rate is larger if a participant was circumcised than not circumcised.
Parametric Accelerated Failure Time Frailty Model
Therefore, for a specific participant, the clearance process for an uncircumcised man is e β times of the clearance process if the man was circumcised, implying circumcision accelerates HR-HPV clearance with a factor of exp( β).
where ξ i ~Normal (0, σ 2), S 0(·) is the baseline survivor function of a parametric survival distribution, such as Weibull distribution, log-normal distribution, generalized gamma distribution, log-logistic distribution, generalized F distribution , and inverse Gaussian distribution . Compared to parametric proportional hazards models where only few distributions have the proportional hazard property, this family of models comprehends a broader class of parametric survival distributions and allows for more flexibility on the shape of the conditional hazard function. It can be shown that AFT models can be expressed in log-linear model form , where it is easy to see that the predictor variables act additively on the logarithm of the survival time T, or equivalently multiplicatively on T itself. Conditional on a given value of the frailty effect ξ i in the AFT model 5, exp (β) has the interpretation of the ratio of the median (or any percentile) survival times between intervention and control arms. Moreover, Additional file 1 shows that the AFT frailty model also provides a marginal interpretation of the treatment effect, where -β is the average log ratio of the clearance times between intervention and control arms, i.e. .
It is important to notice that, however, except for the Weibull distribution (including the exponential distribution), other aforementioned distributions for AFT frailty models do not have the proportional hazards property and thus cannot be modeled as a proportional hazards frailty model given in 1, exp (β) consequently does not have the interpretation of conditional hazard ratio.
where a ij = exp(β 0 + x ij β + ξ i ), and ξ i ~Normal (0, σ 2).
Estimation of model parameters
where S ij (t|ξ i ) is the conditional survivor function for the jth observation in person i, and given ξ i , the conditional survivor functions corresponding to the multiple infections within person i are independent. is the density function for the normal random effect ξ i . The full likelihood is then obtained by integrating over all the possible values of the random effect.
The maximum likelihood estimate (MLE), denoted as and , can be attained by maximizing the likelihood in expression 7 using iterative procedures, and the variance-covariance matrix is estimated as the inverse Hessian matrix. To test the significance of the explanatory variable, i.e. H 0 : β k = 0, the likelihood ratio (LR) statistic (difference in -2Log(L) between the hierarchical models containing and not containing the variable) can be used by comparing it to a χ 2 distribution with 1 degree of freedom . Alternatively, Wald type inference can be drawn with its MLE and variance. For the random effect parameter σ 2, the null hypothesis H 0: σ 2 = 0 corresponds to the situation where all of the observations within a cluster are independent, and thus may be of interest of testing. However, since σ 2 ≥ 0, = 0 is the boundary of the parameter space of σ 2.  showed that the asymptotic distribution of the LR statistic (comparing between the models containing and not containing the random effect) is a 50:50 mixture of the and distribution. Thus a simple rule of computing P-value for testing H 0: σ 2 = 0 is that if the LR statistic is 0, then P = 1; otherwise, the P-value is half of the P-value obtained from comparing the LR statistic to distribution .
Since we assume a normal prior on the random frailty effect, the likelihood function normally does not have a close analytical form, although SAS PROC NLMIXED can numerically compute the integrals and maximize the approximated likelihood. The procedure provides commonly estimated statistics such as the MLE, Wald-type confidence intervals, and -2log-likelihood. In order to perform likelihood ratio inference for a variable, the hierarchical models with and without the variable have to be estimated respectively. The SAS code for analyzing the HR-HPV clearance data using model 2 is listed in Additional file 2 as an illustration, and relevant computation details are also discussed. The format of the input dataset and some useful options needed when calling the procedure are also described.
Results of Study on Male Circumcision Effect on HR-HPV Clearance
Parameter estimates from the Weibull and Log-logistic frailty model with normal random effect for studying circumcision's effect on HR-HPV clearance.
LR P-value (2-sided)
LR P-value (2-sided)
There are several limitations of this application, however. One is that it was found that at the follow-up visits, a significant higher proportion of samples collected on circumcised men could not be amplified rendering more missing HPV results in intervention arm. In this illustration, it was assumed that the infection persisted during the interval if the testing result was missing for the follow-up visit at the end of the interval. This may underestimate circumcision's effect. Techniques for dealing with missing data, such as multiple imputation [18, 19], can be utilized for this specific dataset. Another limitation for this study is that the follow-up intervals of 6, 12 and 24 months were too long to capture the clearance of HPV infection. It has been known that the median duration of genital HPV infection in woman is 8 months, and persistent detectable infection rate is approximately 30% after 1 year and 9% after 2 years . In a recent observational study on men , it is reported that the median time to clearance was 5.9 months (95% CI, 5.7-6.1 months), and 75% infections had cleared by 24 months after initial detection, implying the clearance rate is low when t is large. With the effective immune response to HPV, most infections had cleared by 6 months as observed in the trial participants, meaning that there are no data describing the early phase of the clearance process. The limited data determine that there is "little to choose between alternative distributional models"  and different models may yield similar results in the range where data are observed. It is suggested to adopt the model "most convenient for the purpose in hand " . However, any extrapolation on the functional form of the clearance process from the estimated model should be conducted with caution. This trial was primarily designed to study male circumcision's effect on preventing HIV acquisition. If it is of interest to examine HPV clearance in a future study, the study design should allow for a more frequent testing interval in order to capture the whole process. The presented parametric frailty models implicitly assume the same conditional baseline hazard functions for the clearance of different genotype infections within an individual. The different genotypes all belong to the papillomavirus family, and the mechanism of immune response is the same when fighting against the different type of HPV infections. Therefore it is reasonable to apply the same form of conditional baseline hazard functions for the clearance of different genotype infections.
Conclusion and Discussion
Arbitrarily censored survival data is not uncommon in epidemiological studies, and the censoring nature should be considered during analysis to reduce estimation bias, or when the disease onset and diagnosis are two steps that need to be differentiated . Moreover, clustered data may arise and the correlation within each cluster should also be accounted for. In the current study we particularly describe the use of parametric frailty models to explore treatment effect's on survival when data are clustered and subject to arbitrary censoring. The two main classes of models are parametric proportional hazards frailty model and accelerated failure time frailty model. Most commonly used survival distributions can be used in these models, providing abundant choices of parametric forms to appropriately model the data of interest. For example, Weibull distribution and Gamma distribution can be used for survival problems where the hazard monotonically changes with time, and log-logistic and log-normal distribution can be used when the hazard is a unimodal function of time. The main advantage of adopting a parametric form is for computational ease. On the other hand, with the presence of arbitrary censoring and clustering, it is difficult to perform model diagnostics on the assumption of the parametric form. A clear understanding of the scientific nature of the problem to be addressed is essential for choosing an appropriate parametric distribution in analysis. When using normal random effect, the presented models can be estimated using SAS PROC NLIMIXED, and the code for analyzing the example HPV dataset is provided in the Additional file 2. Models with random effects following other distributions may be estimated using PROC NLIMIXED by transforming the normal random effect using appropriate probability transformation function provided by SAS [21, 22]. Alternatively, for gamma frailty model or log-t proportional hazards frailty model for data with arbitrary censoring, the "frailty()" function provided in the R package "survival" can be used.
Genital HPV infection has high prevalence in both men and women, and the high risk types of HPV are well known to be associated with anogenital cancers, especially cervical cancer . Current diagnostic tools allow for simultaneous detection of multiple HPV genotypes, though the actual infection or clearance time is unknown. Therefore clustered data with arbitrary censoring are normally generated from such studies. The presented modeling approach can be used to study factors associated with HPV clearance (or persistence), or to compare the clearance process between different genotypes to examine type-specific persistence. However, as pointed earlier, the design for such studies need to allow for appropriate short testing intervals to capture the entire process.
We owe thanks to the male circumcision trial participants and the Rakai field work teams, without whom we would not have the example dataset. The Rakai trial was supported by a grant (UO1 AI11171-01-02) from the National Institutes of Allergy and Infectious Disease (NIAID), Division of AIDS, National Institutes of Health (NIH), and in part by the Division of Intramural Research, NIAID, NIH.
- Sun J: The Statistical Analysis of Interval-censored Failure Time Data. 2006, SpringerGoogle Scholar
- Koshiol J, Schroeder J, Jamieson D, Marshall S, Duerr A, Heilig C, Shah K, Klein R, Cu-Uvin S, Schuman P, Celentano D, Smith J: Time to clearance of human papillomavirus infection by type and human immunodeficiency virus serostatus. International Journal of Cancer. 2006, 119 (7): 1623-1629. 10.1002/ijc.22015.View ArticlePubMedGoogle Scholar
- Petersen T: Time-Aggregation Bias in Continuous-Time Hazard-Rate Models. Sociological Methodology. 1991, 21: 263-290. 10.2307/270938.View ArticleGoogle Scholar
- Collett D: Modelling survival data in medical research. Chapman & Hall Ltd: London, New YorkGoogle Scholar
- Lin D: Cox regression analysis of multivariate failure time data: the marginal approach. Statistics in medicine. 1994, 13: 2233-2247. 10.1002/sim.4780132105.View ArticlePubMedGoogle Scholar
- Wei L, Lin D, Weissfeld L: Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. Journal of the American Statistical Association. 1989, 84 (408): 1065-1073. 10.2307/2290084.View ArticleGoogle Scholar
- Hougaard P: Statistical Models and Methods for Biomedical and Technical Systems-Semiparametric Regression Models for Interval-Censored Survival Data, With and Without Frailty Effects. 2008, Boston: Birkháuser Boston, 307-317.Google Scholar
- Tobian A, Serwadda D, Quinn T, Kigozi G, Gravitt P, Laeyendecker O, Charvat B, Ssempijja V, Riedesel B, Oliver A, Nowak R, Moulton L, Chen M, Reynolds S, Wawer M, Gray R: Male Circumcision for the Prevention of HSV-2 and HPV Infections and Syphilis. New England Journal of Medicine. 2009, 360 (13): 1298-1309. 10.1056/NEJMoa0802556.View ArticlePubMedPubMed CentralGoogle Scholar
- Gray R, Kigozi G, Serwadda D, Makumbi F, Watya S, Nalugoda F, Kiwanuka N, Moulton L, Chaudhary M, Chen M, Sewankambo N, Wabwire-Mangen F, Bacon M, Williams C, Opendi P, Reynolds S, Laeyendecker O, Quinn T, Wawer M: Male circumcision for HIV prevention in men in Rakai, Uganda: a randomised trial. Lancet. 2007, 369 (9562): 657-666. 10.1016/S0140-6736(07)60313-4.View ArticlePubMedGoogle Scholar
- Gray R, Serwadda D, Kong X, Makumbi F, Kigozi G, Gravitt P, Watya S, Nalugoda F, Ssempijja V, Tobian A, Kiwanuka N, Moulton L, Sewankambo N, Reynolds S, Quinn T, Iga B, Laeyendecker O, Wawer M: Circumcision of HIV-infected men: Effects on High Risk Human Papillomavirus Infections in a Randomized Trial in Rakai, Uganda. Journal of Infectious Diseases.Google Scholar
- Vaida F, Xu R: Proportional hazards model with random effects. Statistics in Medicine. 2000, 19: 3309-3324. 10.1002/1097-0258(20001230)19:24<3309::AID-SIM825>3.0.CO;2-9.View ArticlePubMedGoogle Scholar
- Ripatti S, Palmgren J: Estimation of multivariate frailty models using penalized partial likelihood. Biometrics. 2000, 56: 1016-1022. 10.1111/j.0006-341X.2000.01016.x.View ArticlePubMedGoogle Scholar
- Kalbfleisch J, Prentice R: The statistical analysis of failure time data. 2002, Hoboken, NJ: John Wiley & SonsView ArticleGoogle Scholar
- Agresti A: Categorical Data Analysis. 2002, Wiley-Interscience, full_text.View ArticleGoogle Scholar
- Self S, Liang K: Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of American Statistical Association. 1987, 82 (398): 605-610. 10.2307/2289471.View ArticleGoogle Scholar
- SAS/STAT 9.2 User's Guide. 2008, Cary, NC: SAS Institute Inc, 2234-2235.Google Scholar
- Morse S, Ballard R, Holmes K, Moreland A: Atlas of sexually transmitted diseases and AIDS. 2003, Elsvier Science, ThirdGoogle Scholar
- Rubin DB: Multiple imputation for nonresponse in surveys. 2004, Wiley-IEEEGoogle Scholar
- Allison PD: Missing data. 2001, SAGEGoogle Scholar
- Giuliano A, Lu B, Nielson C, Flores R, Papenfuss M, Lee J, Abrahamsen M, Harris R: Age-specific prevalence, incidence, and duration of human papillomavirus infections in a cohort of 290 US men. Journal of Infectious Diseases. 2008, 6: 827-835. 10.1086/591095.View ArticleGoogle Scholar
- Lambert P, Collett D, Kimber A, Johnson R: Parametric accelerated failure time models with random effect and an application to kidney transplant survival. Statistics in Medicine. 2004, 23: 3177-3192. 10.1002/sim.1876.View ArticlePubMedGoogle Scholar
- Liu L, Yu Z: A likelihood reformulation method in nonnormal randomeffects models. Statistics in medicine. 2008, 27: 3105-3124. 10.1002/sim.3153.View ArticlePubMedGoogle Scholar
- Pinheiro J, Bates D: Approximations to the Log-Likelihood Function in the Nonlinear Mixed-Effects Model. Journal of Computational and Graphical Statistics. 1995, 4: 12-35. 10.2307/1390625.Google Scholar
- Liu L, Huang X: The use of Gaussian quadrature for estimation in frailty proportional hazard models. Statistics in medicine. 2008, 27: 2665-2683. 10.1002/sim.3077.View ArticlePubMedGoogle Scholar
- Lesaffre E, Spiessens B: On the Effect of the Number of Quadrature Points in a Logistic Random-Effects Model: An Example. Journal of the Royal Statistical Society. Series C (Applied Statistics). 2002, 50 (3): 325-335. 10.1111/1467-9876.00237.View ArticleGoogle Scholar
- Henschel V, Engel J, Hózel D, Mansmann U: A semiparametric Bayesian proportional hazards model for interval censored data with frailty effects. BMC Medical Research Methodology. 2009, 9: 10.1186/1471-2288-9-9.Google Scholar
- The pre-publication history for this paper can be accessed here:http://www.biomedcentral.com/1471-2288/10/40/prepub
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.