On the Addams family of discrete frailty distributions for modeling multivariate case I interval-censored data

Abstract

1 Introduction

Multivariate time-to-event data are commonly encountered in the life sciences. An example is the time to occurrence of multiple non-lethal events within the same individual. In this setting, the individual can be thought of as forming a cluster in which event times are likely to be correlated. An alternative point of view is that there is heterogeneity across individuals (or clusters) due to characteristics that may be difficult or impossible to measure. Random effect (RE) models for time-to-event data, also known as frailty models, offer a conceptually appealing approach to quantify these associations within clusters and to model unobserved heterogeneity across individuals (or clusters) (Hougaard 2000; Duchateau and Janssen 2008; Wienke 2010).

While the majority of the existing literature assumes continuous frailty distributions, such as the gamma (⁠$G$⁠), inverse Gaussian, or log-normal, for some applications discrete frailty distributions may be more appropriate. One such application is the study of infectious diseases transmitted via close contact, where unobserved heterogeneity with respect to becoming infected with pathogens of interest may be represented and analyzed by latent risk groups.

This paper discusses the interpretation of discrete frailty models with parameter dependent support. We propose interpreting the discrete frailties as ordered latent risk categories. We extend this analysis by allowing the parameters of the discrete frailty distribution to depend on covariates. The ordered nature of discrete frailty facilitates comparing latent risk categories within or across different strata (as defined by covariate values) using hazard ratios, along with probabilities of belonging to specific risk categories. This approach facilitates the analysis of non-homogeneous populations or non-homogeneous covariate effects, which might lead to a deeper understanding of the problem under study.

In this context we discuss the “Addams” family (⁠$AF$⁠) of discrete frailty distributions introduced by Farrington et al. (2012). This family is homogeneous in the sense that it consists of discrete distributions, except for the continuous $G$ distribution, which serves as a reference case of constant association within clusters over time. The discrete distributions within the $AF$ are scaled variants of negative binomial, shifted negative binomial, binomial, and Poisson distributions. For this family, we investigate the hazard ratios conditional on risk group membership within or across different strata. Furthermore, we introduce estimation routines for the $AF$ in the context of case I interval-censored data (Sun 2007). Our optimization framework allows for stratification of frailty distributions along covariates. We apply the $AF$ to data from a serological survey of human papillomaviruses (Mollema et al. 2009; Scherpenisse et al. 2012).

The structure of the paper is as follows: Section 2 presents the time-invariant shared frailty model and reviews the literature on discrete frailty models, comparing the $AF$ to other discrete frailty distributions. We further discuss the interpretation of discrete frailty models for which the support is parameter dependent. Section 3 examines the conditional time-to-event model resulting from the $AF$ of distributions. Section 4 outlines the estimation framework and optimization algorithm. In Section 5, the $AF$ is applied to multivariate serological survey data. Finally, Section 6 offers concluding remarks.

2 Discrete frailty models

Let i refer to a cluster, $i=1,…,n$⁠, and let n be the number of clusters observed. Each cluster contains J units, and $Ti(j)>0$ denotes the time-to-event random variable (RV) for the $jth$ unit, $j=1,…,J$⁠, within the $ith$ cluster. All indices are unique, such that $i≠i′$ and $j≠j′$⁠. The random vectors $Ti=[Ti(1),…,Ti(J)]T$ and $Ti′$ are independent given the covariates $xi$ and $xi′$⁠. Within a cluster, $Ti(j)$ and $Ti(j′)$ are independent given the random and unobservable time-invariant cluster-specific frailty $Zi$ and covariates. The non-negative RV $Zi$ has density or probability mass function $g(z|x˜i)≡gi(z)$⁠, where the index i is relevant only if the distribution or its parameters depend on (a subset of) covariates denoted by $x˜i$⁠. Technically, $x˜i$ could include unit-specific covariates. However, as $Zi$ is shared within the cluster i, $x˜i$ will typically also be shared within the cluster i, without containing covariates that differ across the units. For unit-specific covariates, correlated frailties (Hens et al. 2009) may be more appropriate, where the correlation parameter between the frailties of a cluster may depend on cluster-invariant covariates, while unit-specific frailty parameters may depend on unit-specific covariates. Given covariates, $Zi$ and $Zi′$ are independent. We do not distinguish in language between the random $Zi$ and a realization $zi$ and refer to both as frailty or RE.

In much of the existing literature, $Zi$ is typically considered a continuous random variable, with common choices being the log-normal or $G$ distributions. Nonetheless, discrete frailty distributions have also been explored for both univariate and multivariate data. A prominent choice is the K-point distribution (Pickles and Crouchley 1994; Begun et al. 2000; Choi and Huang 2012; Bijwaard 2014; Choi et al. 2014; Palloni and Beltrán-Sánchez 2017; Troncoso-Ponce 2018; Gasperoni et al. 2020), ie a frailty distribution with (ordered) support parameters ${z(1),…,z(K)}∈R≥0K$ and corresponding probability parameters $prm=gi(z(m))>0$ for $m=1,…,K$⁠, and $∑m=1Kprm=1$ (Wienke 2010). Other choices are the binomial (⁠$B$⁠) (Ata and Özel 2013), negative binomial (⁠$NB$⁠) (Caroni et al. 2010; Ata and Özel 2013), geometric (Caroni et al. 2010; Choi and Huang 2012; Choi et al. 2014; Cancho et al. 2021), Poisson (⁠$P$⁠) (Caroni et al. 2010; Choi and Huang 2012; Ata and Özel 2013; Choi et al. 2014; Cancho et al. 2020b), and the hyper-Poisson distributions (de Souza et al. 2017; Mohseni et al. 2020). The framework of the zero-inflated and zero-modified power series (ZMPS) distributions, where the probability of $Z=0$ is modified by an additional parameter as compared to the discrete reference distribution, has been investigated by Cancho et al. (2018, 2020a) and Molina et al. (2021), respectively. In particular, the zero-inflated geometric, $P$ and logarithmic distributions and the zero-modified geometric and $P$ distributions are investigated. The Addams family (⁠$AF$⁠), as conceptualized by Farrington et al. (2012), includes (shifted and) scaled $NB$⁠, as well as scaled $B$ and $P$ distributions, and the $G$ distribution as a continuous special case.

Modeling a whole family of distributions, such as the $AF$⁠, is preferable to modeling a given (discrete) distribution, such as the $P$⁠, since the latter strategy typically severely limits the patterns of heterogeneity as measured by the relative frailty variance (⁠$RFV(t)=Var(Z|T>t)E(Z|T>t)2$⁠,$t=[t(1),…,t(J)]T$⁠) to either monotonically increasing or monotonically decreasing trajectories over time; for examples see Farrington et al. (2012) and Bardo and Unkel (2023). It can be shown that the long-term trajectory of the RFV (or association, as measured by the cross-ratio function [$CRF(t)=1+RFV(t)$]) is determined with positive probability by the smallest value of the frailty distribution (⁠$z(1)$⁠). Specifically, if $z(1)=0$⁠, $RFV$ and $CRF$ approach infinity as time approaches infinity. Conversely, if $z(1)>0$⁠, the RFV approaches zero and the CRF approaches one. Therefore, the choice of a discrete frailty distribution has an enormous impact on the model, even if $gi(0)>0$ or $z(1)>0$ is very small, and even more so if the model’s trajectory of RFV is monotone. The $AF$ achieves greater flexibility in the trajectory of heterogeneity and association by incorporating discrete frailty distributions for which $z(1)=0$ and frailty distributions for which $z(1)>0$⁠, and is thus able to induce increasing and decreasing trajectories of the RFV (CRF). Thus, in the case of $AF$⁠, the trajectory of the RFV can be informed by the data in the fitting process. This is a rare property among frailty distributions, and for the discrete distributions mentioned above it is only possible for the K-point distribution and the ZMPS. However, in these cases the decision between a decreasing or increasing long-term trajectory lies on the edge of the parameter space (Bardo and Unkel 2023), which is not the case for the $AF$⁠. As optimizing on the edge of the parameter space is usually difficult, this is an advantage when using the $AF$⁠. Note, however, that the K-point distribution and the ZMPS are able to induce non-monotone trajectories of the RFV (CRF), which is not the case for the $AF$⁠. We will discuss the $AF$ and its special cases in Section 3.

From an interpretive point of view, the $AF$ and K-point distributions offer a unique perspective because the support is also subject to estimation. The support of other discrete frailty distributions is usually the natural number including zero (⁠$N≥0$⁠). Therefore, the interpretation of discrete frailty models often focuses on the cure rate, ie those who are not susceptible to the event of interest [see e.g. de Souza et al. (2017); Cancho et al. (2018, 2020a, 2020b); Mohseni et al. (2020); Cancho et al. (2021); Molina et al. (2021)]. Subject-related interpretations are also common. Caroni et al. (2010) suggest interpreting discrete frailties as the unobservable number of flaws in a unit or exposure to damage on an unknown number of occasions. Similar interpretations can be found in Ata and Özel (2013) for time-to-event models on earthquake data. Both studies consider discrete frailty distributions with support on $N≥0$⁠. However, due to the latent nature of the frailty, such concrete interpretations are difficult because the hazard ratio (HR) of, say, two events versus one event would be fixed by the model structure at $HR=2$ if the support is fixed at $N≥0$⁠. However, if the true HR is less than two, the probability mass of frailty could be shifted to the right relative to the distribution of the number of hits. Consequently, a more abstract interpretation such as the “effective” number of harms would be more appropriate. The K-point frailty distributions are often interpreted as representing sub-populations, such as unobservable carriers of certain disease genotypes [e.g. Pickles and Crouchley (1994); Begun et al. (2000); Wienke (2010); Bijwaard (2014)]. Palloni and Beltrán-Sánchez (2017) interpret a delayed binary frailty as the effect of adverse early life conditions on adult mortality. Using the K-point distribution, where K is also subject to estimation, Gasperoni et al. (2020) interpret the discrete frailties as (an unknown number of) latent sub-populations. They suggest calculating HRs between these sub-populations, which is interpretable as the support being subject to estimation.

In the present paper, we endorse the interpretation of discrete frailties as latent sub-populations and expand upon it. Note that we focus on discrete frailty distributions for which the support $Ω$ is parameter dependent, which is mainly the case for the K-point distribution and, although more restricted, the $AF$⁠, as will be seen in Section 3. Let $Ωi={zi,(1),zi,(2),…}$⁠, with $zi,(1)≥0$ and $zi,(k)<zi,(k+1)$⁠, represent the support of discrete RV $Zi$⁠, where the distribution parameters might depend on $x˜i$⁠. We define that $x˜i$ constitutes a stratum of the population. We also consider $zi,(k)$ as the conditional hazard-determining value for an individual in the $kth$ risk category (RC) within the stratum which is defined by $x˜i$⁠.

For discrete frailty distributions for which the support depends on parameters, the within-stratum hazard ratio, $HRW(k)=λi(j)(t|Zi=zi,(k+1))λi(j)(t|Zi=zi,(k))=zi,(k+1)zi,(k)$⁠, might be analyzed. The $HRW(k)$ compares the hazard of the $kth$ and $(k+1)th$ RC of stratum i. This is in line with Gasperoni et al. (2020), except for the presence of different strata for the distribution of the frailty, ie the $HRW$ might differ for $x˜i$ and $x˜i′$⁠. Moreover, due to the presence of different strata for the frailty distribution, an across-stratum analysis can be conducted by computing the across-stratum hazard ratio, $HRA(k)=λi(j)(t|Zi=zi,(k))λi′(j)(t|Zi′=zi′,(k))=zi,(k)zi′,(k)$⁠, for $x˜i≠x˜i′$⁠, and equality in the remaining covariates $xi(j)$ and $xi′(j)$⁠. If the support between stratum i and $i′$ is very different, it might be more desirable to analyze $HRA(k;k′)=zi,(k)zi′,(k′)$ for all $k′$ for which $P(Zi′≤zi′,(k′))∈[P(Zi≤zi,(k−1)),P(Zi≤zi,(k)))$ or the closest quantiles of $Zi,Zi′$ if no such $k′$ exists. In order to put the analysis of the $HRW$ and $HRA$ into context, they should always be accompanied by reporting the distribution of the RCs. This allows for a separate but accompanying analysis of the distribution of individual heterogeneity via the distribution of RCs and the magnitude related impact of the RCs on the hazard rates.

The approach of imposing a model on the distribution parameters of the frailty has some similarity to generalized additive models for location, scale, and shape parameters for the population time-to-event distribution, ie for the time-to-event distribution with the frailty marginalized out. However, modeling the determinants of the randomness of hazard rates via covariates might be more intuitive, as the hazard rates are usually the standard approach for modeling time-to-event data. This approach is not new per se. It can be found, for example, in Aalen et al. (2008) and is also quite common in the field of discrete frailty modeling (Choi and Huang 2012; Choi et al. 2014; Cancho et al. 2018, 2020a; Molina et al. 2021), and random slopes could also be interpreted in this way. What is new, however, is the type of analysis with within-stratum and across-stratum hazard ratios. Note that this type of analysis has no counterpart for continuous frailty distributions, as there is no $kth$ RC. Nor does such an analysis make sense for discrete frailty models where the support of the frailty distribution is set by assumption, e.g. to $N≥0$⁠, as this would fix the $HRW$ and $HRA$ by assumption. A tempting alternative in such a case might be to calculate the hazard ratio between, say, the third and the first quartile of the frailty distribution. However, the quartiles of the frailty distribution vary over time due to selection effects, which is particularly important for hazard ratios which inherit the survival condition. Therefore, the hazard ratio is either meaningless if the quartiles at $t=0$ are chosen and maintained throughout the analysis, or the analysis is more complicated if the time-varying quartiles of the frailty distribution are chosen.

Therefore, an analysis via the $HRW$ and $HRA$ is unique to discrete frailty models, where the support of the frailty distribution varies with (a subset of) distribution parameters that may depend on stratum membership. In this case, the data inform the support of the frailty distribution by likelihood criteria making it suitable to represent latent RCs. This allows the analysis of non-homogeneous covariate effects using common measures such as HRs as described above and probabilities of belonging to a particular RC. This might in particular be helpful in communicating the results of heterogeneous (covariate) effects to an audience outside of statistics such as medical doctors.

3 The Addams family of discrete frailty distributions

Let $Λ0(j)(t)$ denote the cumulative baseline hazard rate $∫0tλ0(j)(u)du$⁠. Moreover, $Λi(t)=∑j=1J exp {xi(j)Tβ(j)}$$Λ0(j)(t(j))$⁠. We further ignore that the parameters of the frailty distribution may depend on $x˜i$ in the first part of this section and come back to this issue in the latter part of this section.

Note that $μ$ may be set to one for identification purposes. Hence, the $AF$ is a two parameter distribution. We will retain the parameter $μ$ in our notation, as we will explicitly model it via $μ(x˜)=exp {x˜Tβ(0)}$⁠, where $β(0)∈R$ is an additional vector of parameters. In a setting where the distribution of the frailty is not the focus of the analysis, $β(0)$ is typically interpreted as the log-proportional hazard ratio.

The Laplace transform uniquely determines the distribution of Z, as shown in Table 1 (Farrington et al. 2012). The $AF$ consists of different scaled and possibly shifted discrete distributions. The scaling parameter $ψ$ of the corresponding discrete distributions is equal to $μ|α|$⁠, and the parameter $α$ selects one of the distributions of the $AF$⁠. If $α<0$⁠, the scaled and shifted negative binomial (⁠$ψNB>0$⁠) is chosen, where the support is shifted to the right by the parameter $ν$⁠, resulting in a model without a latent non-susceptible sub-population. If $α>0$⁠, a frailty distribution is chosen that includes a non-susceptible latent sub-population. For $γ>α>0$⁠, the non-shifted scaled negative binomial (⁠$ψNB$⁠) and for $α=γ$ the scaled Poisson (⁠$ψP$⁠) distribution is selected. In the case of $α>γ$⁠, the frailty distribution is scaled binomial (⁠$ψB$⁠), resulting in a model with an upper bound of the frailty that limits the maximum deviation of RC-conditioned survival curves between the susceptible latent sub-populations and a finite number of RCs including a non-susceptible group. Note that for the case $α>γ$⁠, $(α−γ)−1$ has to be an integer in order for $L(s)$ to be a valid Laplace transform (Farrington et al. 2012). The continuous exception in the $AF$ is the $G$ distribution, which results from $α=0$⁠.

The parameter $α$ plays a unique role in the context of discrete shared frailty modeling. For discrete shared frailty models the RFV (CRF) either approaches zero (one) or infinity (Bardo and Unkel 2023). Therefore, it is desirable to have a continuous exception within a family of discrete shared frailty distributions for which the RFV (CRF) does not approach zero (one) or infinity with time approaching infinity. Within the context of the $AF$ this is the $G$ which arises for $α=0$⁠. In that case the RFV (CRF) is constant, a shape that is impossible for a discrete shared frailty model to generate. The nested structure might be utilized to test for a constant RFV (CRF) within the $AF$⁠.

Moreover, the $AF$ can choose between a monotonically decreasing or increasing trajectory of the RFV (CRF) through the sign of $α$⁠. This is unprecedented in the context of discrete shared frailty models. Though the ZMPS distribution is able to create decreasing trajectories of the RFV (CRF), this involves the edge of the parameter space for its deflation/inflation parameter. The opposite is true for the K-point distribution: its RFV (CRF) approaches zero (one) in the long run unless $z(1)$ is equal to zero, which again involves the edge of the parameter space. For the $AF$⁠, the parameter $α$ chooses between a decreasing or increasing RFV (CRF) by determining the support of the discrete distribution without involving the edge of parameter space. If $α>0$⁠, $z(1)=0$ and a cure fraction exists. This induces an increasing trajectory of the RFV (CRF). If $α<0$ instead, $z(1)>0$ and no cure fraction exists. This induces a decreasing trajectory of the RFV (CRF); see Bardo and Unkel (2023) for a discussion of the shape of the RFV (CRF) for discrete frailty models.

The feature of the support being dependent on the distribution parameters offers the possibility of a meaningful interpretation of the $HRW$⁠. Figure 1a shows examples of the $HRW$ for $α>0$ and $α<0$⁠. If $α>0$⁠, $HRW(1)≡∞$ and $HRW(k)=kk−1$ for $k≥2$⁠. Note that there may be an upper bound on the RCs if $ψB$ is chosen. In this case, the upper bound of the frailty as well as the number of RCs chosen through the fitting procedure may be the main component of the analysis, e.g. by comparing the RC-related survival curve of the upper bound versus another RC. If $α<0$⁠, $HRW(k)=ν+kν+k−1$ which approaches $kk−1$ for large k. So there is reasonable flexibility for the first few within-stratum HRs and the focus may be on $HRW(1)$⁠, where the model is more flexible than for later RCs. This is less flexible than the K-point distribution is, provided that K is large enough, which may even show a non-monotone trajectory of the $HRW$⁠. However, modeling the K-point distribution with a large K is difficult given that this involves $2(K−1)$ parameters for a latent distribution, especially if the cluster size is small. On the contrary, within the $AF$ one is remunerated with $|Ω|=∞$ (except for the $ψB$ case) which can be important for extreme observations where events occur very early in the lifespan. However, the parametric constraints of the $AF$ on the trajectory of $HRW$ should always be taken in consideration and challenged by the K-point distribution whenever possible.

Figure 1:

Within- and across-stratum hazard ratios versus RC k for various members of the Addams family. a) Within stratum HR for $α>0$ and $α<0$ versus RC k. b) Across-stratum HR for $αi,αi′>0$ versus RC k. c) Across-stratum HR for $αi,αi′<0$ versus RC k. d) Across-stratum HR for $αi<0$⁠, $αi′>0$ versus RC k.

Open in new tabDownload slide

This analysis can be extended to the across-stratum HR, $HRA$⁠, if there is a model for the distribution of individual heterogeneity within the $AF$⁠. For that purpose, we allow the parameters to depend on a stratum-specifying set of covariates $x˜$⁠, which determine the parameters of individual heterogeneity $α(x˜)=x˜Tζ,γ(x˜)=exp {x˜Tκ}$⁠, with each parameter in the vectors $ζ,κ$ being an element of $R$⁠. For the sake of brevity, we denote $α(x˜i),γ(x˜i),$ and $μ(x˜i)$ by $αi,γi,$ and $μi$⁠, respectively. Figure 1 shows the $HRA$ for varying scenarios of $αi$ and $αi′$⁠. For $αi,αi′>0$ (Fig. 1b), the $HRA=ψiψi′$ for all $k≥2$ and is undefined for $k=1$⁠. However, for $αi,αi′<0$ (Fig. 1c), the $HRA(k)=ψi(νi+k−1)ψi′(νi′+k−1)$⁠, which approaches a constant ratio $HRA(k)=ψiψi′$ for large k. Note that the $HRA$ might be greater or less than one for all k, but can also cross the threshold of one with increasing k. If $x˜$ is for example an experimental treatment indicator (in a univariate context), the $HRA$ represents a heterogeneous treatment effect which might identify sub-groups within the population for which the treatment is harmful. For $αi<0,αi′>0$ (Fig. 1d), $HRA(1)≡∞$⁠, as for stratum $i′$ there is a latent sub-population that is not susceptible to the event of interest, whereas for stratum i all latent sub-populations are susceptible. For $k≥2$ and $HRA(k)=ψi(νi+k−1)ψi′(k−1)$⁠. Another scenario, not explicitly shown in Fig. 1b and d, is that one or both strata (i and $i′$⁠) might have (different) upper bounds of frailty in the $ψB$ case. In such a case, the stratum with the larger upper bound (which could still be $∞$⁠), say $i′$⁠, could be considered more vulnerable, since that stratum has a higher proportion of individuals who are expected to experience the event very early, namely those with a frailty value greater than $ψibi$⁠.

Note that the parameters in the formula of $HRA$ and $HRW$ (and hence the parameters as specified in the legends of Fig. 1) do not uniquely identify the parameters of the frailty distribution, ie for a given trajectory of $HRA,HRW$ there is an infinite set of $(α,γ)$ or $(αi,γi)$ and $(αi′,γi′)$⁠, respectively, that induce the same trajectory but with a different distributions of the RCs which did not need to be specified for Fig. 1. This shows that the analysis of the frailty model has always two branches. The first branch is the analysis of HRs, which indicate the meaning of being in a particular latent RC (in a particular stratum) relative to another latent RC or to another observable stratum in the same latent RC. On the one hand, $HRW$ can help to assess the importance of individual heterogeneity, e.g. if $HRW(k)$ is large, then latent RC membership has a large effect on expected survival. On the other hand, comparing $HRW$ across strata or analysing the $HRA$ gives an account of random covariate effects where, e.g. covariates with a beneficial effect on survival or covariates with partly beneficial, partly detrimental effects might be detected. The second branch of the analysis is the distribution of RCs across strata, which may indicate differences in the distribution of risk-taking behavior and predisposition across strata, e.g. by indicating a stratum with a heavier tail of vulnerable RCs. Taken together, the analysis of HRs and the distribution of RCs can provide thorough analytical explanations in terms of selection and random covariate effects that can help to explain the trajectories of population survival curves (where the RCs are marginalized out), ie explanations for why the survival curves of two strata come closer or even cross over time.

4 Estimation

In this section, all time-dependent quantities are evaluated at the monitored (censored or uncensored) event times of the individuals. We delete the argument from the expression and indicate the corresponding quantity with a subscript, e.g. $Λi(j)=exp {xi(j)Tβ(j)}Λ0(j)(ti(j))$⁠. Furthermore, let$A∈P({1,…,J})$⁠, where $P$ denotes the power set. Then, $Λi(A)=∑j∈AΛi(j)$ and $Λi(−A)=∑j∉AΛi(j)$⁠. Note that we define $Λi(∅)=0$ and $Λi=Λi(1,…,J)$⁠.

Quasi-Newton optimization routines were applied for optimizing the corresponding log-likelihood based on (3). We choose BFGS as the standard method, as implemented in R version 4.2.2 (R Core Team 2022). Standard errors (SE) are obtained via the Hessian of the log-likelihood, which is approximated by Richardson extrapolation as implemented by Gilbert and Varadhan (2019). The delta method was applied where necessary to obtain SE: confidence intervals (CIs) are based on $ln$ or $ln{−ln}$ transformations if the parameter of interest is greater than zero or between zero and one respectively, and are then transformed back to the scale of interest.

We provide algorithms that are able to fit univariate and multivariate frailty models for case I interval-censored data. The frailty distributions can be stratified by a (multi-level) factor. The frailty distribution might either be $AF$ or from the power variance family (both parameters can be estimated). The baseline hazard can be chosen to be piecewise-constant or the parametric generalized gamma distribution (Cox et al. 2007) or one of its special cases, respectively. Covariates can be added in proportional hazards manner. Overdispersion parameters might be added by means of the Dirichlet compound multinomial distribution. Implementations are available on GitHub (https://github.com/time-to-MaBo/Addamsfamily/).

5 Applications

We illustrate the $AF$ in the context of multivariate case I interval-censored data on the human papillomavirus (HPV), obtained from a serological survey in the Netherlands (PIENTER-2); see Mollema et al. (2009) for details on PIENTER-2 and Scherpenisse et al. (2012) for an investigation of the respective HPV dataset. The data were collected in the years 2006 and 2007 and cover people aged 0 to 79. Participants were asked to complete a questionnaire and to provide a blood sample (Mollema et al. 2009). By means of the blood samples, the level of antibodies regarding the high-risk HPV types 16, 18, 31, 33, 45, 52, and 58 were determined in order to detect past infections. Therefore, at the time of observation, it is only known whether the study participants have had an infection in the past or not, but it is never known exactly when the potential infection occurred, resulting in case I interval-censored data, also known as current status data (Sun 2007). Note that at the time of data collection the Dutch national immunization programme did not include a vaccine against HPV.

We analyzed the nationwide sample including oversampled migrants and applied weighting factors to make the sample representative for the Dutch population. We excluded individuals in their first year of life from the analysis, as maternal antibodies could be transmitted to the infant transplacentally or through breastfeeding (Rintala et al. 2005). This left us with a sample size of $n=6384$ individuals aged 2 to 80. The weighted proportion of females in the dataset is 49.9% (unweighted: 54.4%).

The observed time is the individuals’ age at the date of serological monitoring. The event indicator $di(j)=1$ means that individual i is seropositive with respect to pathogen j, $j∈{HPV16,HPV18,HPV31,HPV33,HPV45,HPV52,HPV58}$⁠, and seronegative and still susceptible if $di(j)=0$⁠. Seroprevalence is interpreted as a proxy for past infections. Note, however, that there is a time-lag between the infection and the time of seroconversion, as well as a difference in the number of individuals who were infected with HPV and those who seroconverted: in previous studies, antibodies could not be detected for about 20-50% of females who were carriers of HPV DNA. However, antibody responses are relatively stable over time and hence, the study of the population’s seroprevalence might yield important insights; see Scherpenisse et al. (2012) for a discussion.

We consider the following models for the individual hazard rates, $λi(j)(t)=Ziλ0(sex:j)(t)$⁠,$sex∈{m,f}$⁠, where the target-specific baseline hazard $λ0(sex:j)(t)$ is either sex-stratified (sex-stratified baseline hazard model), or non-stratified $λ0(sex:j)(t)=λ0(j)(t)$ (non-stratified baseline hazard model). The purpose of stratifying baseline hazards by sex is 2-fold. The first is to investigate whether it is justified to estimate baseline hazards jointly for both sexes, and the second is to investigate whether a potential difference in the distribution is better explained by stratified baseline hazards than by different distributions of individual heterogeneity. In any case the target (and potentially sex) specific hazard rate is piecewise constant with a unique parameter within the intervals $[0;5)$⁠, $[5;10)$⁠, $[10;20)$⁠, $[20;30)$⁠, $[30;40)$⁠, $[40;50)$⁠, $[50;65)$⁠, $[65;80)$⁠. The frailty Z is either sex-stratified, ie $Zi∼AF(αsexi,γsexi)$ (sex-stratified RE model), or non-stratified, ie $Zi∼AF(α,γ)$ (non-stratified RE model). Note that in both cases $μm≡1$⁠, $μf∈R>0$ except for the stratified-hazard model where $μf$ is also set to one for the sake of identifiability. The stratified RE model might be better able to reflect differing patterns in individual heterogeneity due to biological and environmental predisposition as well as a different distribution of risk-related behavior across males and females. We combine the stratification status of the baseline hazard with the stratification status of the RE.

An HPV infection can be transmitted via skin-to-skin contact, often—though not exclusively (see, e.g. Syrjänen (2010), Rintala et al. (2005) or Meyers et al. (2014))— via sexual intercourse (Gavillon et al. 2010). Therefore, individual—and typically unobserved—behavior is an important determinant of an individual’s risk of contracting HPV. Frailty models have been used previously to incorporate unobservable individual heterogeneity in the transmission of infectious diseases (see, e.g. Unkel et al. (2014) or Hens et al. (2009)). Moreover, we suspect that there may be distinct jumps in the individual hazard rates due to differences in individual behavior that are relevant for transmission, e.g. comparing individuals who have no sex at all to individuals who have (see, e.g. Richardson et al. (2000), Burchell et al. (2006)), or whether the individuals use condoms or not (see, for example, Lam et al. (2014) or Nielson et al. (2010)). More formally, non-Gaussian and discrete omitted covariates may be the most important drivers of individual heterogeneity, and thus a discrete frailty model may be particularly appropriate here.

We start with a bivariate analysis including HPV16 and HPV18. An extension to higher variate data with data on seven types of HPV follows.

5.1 Bivariate data analysis

We begin by exploiting the nested structure of the models for model selection. The stratification of the RE is statistically significant on conventional levels by means of a likelihood ratio test (LRT) no matter the stratification status of the baseline hazard. The null-hypothesis is $H0$⁠: $Zm∼Zf∼AF(α,γ)$ vs. $H1:Zsex∼AF(αsex,γsex)$⁠. Note that the expectation parameter $μf$ is not included in the hypothesis. In the case of non-stratified baseline hazards the LRT test statistic equals 29.700 on 2 degrees of freedom (p-value $≈$ 0). In the case of sex-stratified baseline hazards the LRT test statistic equals 30.822 on 2 degrees of freedom (p-value $≈$ 0). Better performance of stratified RE models is also suggested by the $ϕ$-plot which can be seen in Fig. 2a. The measure $ϕ$ is an association measure for bivariate current status data introduced by Unkel and Farrington (2012), $ϕ>0$ indicates positive and $ϕ<0$ negative association. Additionally, $ϕ$ tracks $ln{1+RFV(t)}$ with a time-lag. It can be observed that the association between HPV16 and HPV18 is higher for females early in life, but declines more strongly than for males. This is likely to be the reason for the success of stratified RE models here as those models are able to choose a distinct intercept and slope of the RFV across the sexes. We choose the sex-stratified RE model for further analysis.

Figure 2:

Observed association between HPV types in the PIENTER-2 data. a) Observed association between HPV16 and HPV18 in terms of a $ϕ$-plot. Black dots refer to cohort- and sex-specific nonparametric estimates, size proportional to precision. The black solid line is the corresponding LOESS. Other dotted and dashed lines are estimates resulting from corresponding parametric model. b) RFV of bivariate (left) and higher variate (right) dataset of non-stratified hazard, stratified RE model. Note the different scales on the y-axis. The curves reach a plateau at around 30 years of age, hence, the x-axis was cut-off after 40 years. Note that also the seroprevalence curves $(1−P(T(j)>t))$ reach a plateau around this time (not shown).

Open in new tabDownload slide

In terms of AIC, the non-stratified baseline hazard model $λ0(j)(t)$⁠, $j∈{HPV16,HPV18}$ performs better than sex-stratified-baseline hazard model $λ0(sex:j)(t)$ (9647 vs. 9656). Thus, we choose the non-stratified baseline hazard, stratified RE model for further analysis.

A LRT for a constant RFV (CRF), ie $H0:Zsex∼G(γsex)$ vs. $H1:Zsex∼AF(αsex,γsex)$ for males and females, yields a test-statistic of 93.514 on 2 degrees of freedom (p-value $≈0$⁠). Hence, the hypothesis on a constant RFV (CRF) can also be rejected. A constant pattern of the RFV is also not suggested by Fig. 2a, where it can be observed that association is constantly falling up to the age of around 50. Considering all tests and pairwise AIC comparisons we choose the non-stratified baseline hazard, stratified $AF$ model for final analysis.

The RFV parameter estimates for the stratified RE, non-stratified baseline hazard models can be seen in Table 2. The estimated RFV (CRF) is decreasing for both sexes. The estimated RFV parameters indicate higher heterogeneity across clusters or association within a cluster for females early on, as indicated by the intercept of the RFV (⁠$γ^f>γ^m$⁠). However, the descent is more strongly for females (⁠$|α^f|>|α^m|$⁠) and consequently it is estimated that heterogeneity/association is stronger for males from the 4$th$ yr of life onwards; see left-hand panel of Fig. 2b. These results are also supported by the non-parametric estimates of $ϕ$ in Fig. 2a.

The estimated distribution corresponds to a $ψ^sexNB>0(ν^sex,π^sex)$ for males and females. The mean parameter $μ^f=$$0.328 (95%-CI [0.091;1.182])$ and is insignificant as indicated by the $95%$-CI. The resulting distribution parameters can also be found in Table 2. The estimated mean $μ^f$ indicates lower expected frailty (and therefore lower population hazard) for females initially. However, the mean parameter has to be interpreted in the context of its distribution. Let $μsex(t)=E(Z|T>t,sex)$⁠. The limit of the conditional expectation of the frailty is $μsex([∞,∞])=ψsexνsex$⁠. With the estimates from Table 2, $μ^f([∞,∞])=$ 0.01 $>μ^m([∞,∞])=$ 0.006 follows and the initial order of the expectations is reversed. In this example, this leads to the estimated population seroprevalence $P^(T(HPV16)≤age)$ being higher for males early in life but from 12 yrs of life onwards, females start to catch up and finally cross the curve of males at 25 yrs of life (not shown). We will discuss the reason for the switching order of $μ^f(t)$⁠, and $μ^m(t)$ that finally leads to crossing seroprevalence curves by analysing the distribution of the frailties in the paragraphs below.

Table 3 shows an excerpt of the distribution of the RCs. We interpret the distribution of the RCs as the distribution of stratum-relative risk-related behavior and predisposition. The bulk of the population is estimated to be in the lowest RC, though there is more probability to the right of the lowest RC for males. The ratio of the cumulative probabilities between females and males is always above one, also indicating a more heavy tail for males. The heavier tale of the distribution of latent RCs for males is the reason for $μm>μ^f$⁠.