Example: Effect of Male Circumcision on HR-HPV clearance

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).

Notations

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

Let h(t) denote the hazard function representing the "hazard" (i.e. the instantaneous rate) of clearance at time t. To examine the treatment effect on clearance, we can model the hazard to be a function of the explanatory variables, while at the same time including a random effect to account for the correlation between the multiple HR-HPV genotype infections on such infected participants. For example, as in [7], for clearance of the jth infection on person i, we can use the conditional proportional hazards model (or semiparametric frailty model):

(1)

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 ₀(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].

Direct estimation of the coefficient parameters from model 1 by maximizing the likelihood may not be available with standard software such as SAS. To reduce the computation complexity, we further assume a parametric form on the clearance time for participants in the control arm, for example, conditional on the random effect ξfor a person in the control arm, let time to clearance of any HR-HPV infection have a Weibull distribution. That is h ₀(t) = γα ₀ t ^γ-1, where γ(γ > 0) and α ₀ are the shape and scale parameters respectively for the Weibull distribution. The hazard function for Weibull distribution is a monotone function of t [13] (Figure 1a). The Gompertz distribution introduced for describing human mortality [4] is another parameteric distribution that has the proportional hazards property, and can be considered in the proportional hazards frailty model 1 .

Plugging h ₀(t) into expression 1, for the jth HR-HPV infection in person i with explanatory variables x _ij , j = 1, 2, ⋯, n _i , i = 1, 2, ⋯, n, the Weibull frailty model for clearance is:

(2)

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

The Weibull frailty model given in 2 can be equivalently expressed in terms of the survivor function of the Weibull distribution as:

(3)

(4)

where S ₀(t) is the baseline survivor function of the conditional Weibull distribution, that is, S ₀(t) = exp(-α ₀ t ^γ ), where α ₀ = exp(β ₀). Thus we have .

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( β).

Conditional on the random effect, expression 4 corresponding to Weibull distribution in fact belongs to the family of parametric accelerated failure time (AFT) models [4]. For clustered survival data, the general AFT model with normal random frailty effect can be written as:

(5)

where ξ _i ~Normal (0, σ ²), S ₀(·) 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 [13], and inverse Gaussian distribution [4]. 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 [4], 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.

For the HR-HPV data, the proportional hazard frailty model 2 with Weibull distribution assumption implicitly imposes a constant instant clearance rate ratio (conditional on the random effect) over time between intervention arm and control arm (Figure 1b). The clearance rate ratio may however change over time. From Table 1, the majority of infections had cleared by year 2 in either arm, thus the clearance rate ratio should be close to 1; whereas the rate ratio by month 6 may be different from 1. The log-logistic distribution was also fit for the survivor function S(t|ξ _i ) in the AFT frailty model 5. The hazard of log-logistic distribution is a unimodal function of t when its shape parameter is larger than 1 (Figure 1a), which may better capture the natural history of HPV infection. It also allows the clearance rate ratio between arms to change over time(Figure 1b). The log-logistic AFT random frailty model is:

(6)

where a _ij = exp(β ₀ + x _ij β + ξ _i ), and ξ _i ~Normal (0, σ ²).

Estimation of model parameters

With parametric assumptions in model 1 or model 5, parameters (including fixed effects β and the variance of the random effect σ ²) can be estimated using the maximum likelihood method. Recall that for the jth infection in the ith participant, the clearance time is observed to be in the interval of (a _ij , b _ij ], where a _ij = 0 for left censored observations and b _ij = ∞ for right censored observations. For the n _i genotype HR-HPV infections in person i, their clearances can either be left, right, or interval censored. Without loss of generality, let l _i denote the number of left censored observations and r _i denote the number of right censored observations, and hence n _i -l _i -r _i is the number of interval censored observations on person i. Thus the full likelihood on all participants can be written as:

(7)

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 ₀ : β _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 χ ² distribution with 1 degree of freedom [14]. Alternatively, Wald type inference can be drawn with its MLE and variance. For the random effect parameter σ ², the null hypothesis H ₀: σ ² = 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 σ ² ≥ 0, = 0 is the boundary of the parameter space of σ ². [15] 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 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 [16].

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.