Biometrika. 2010 March; 97(1): 133–145.
PMCID: PMC3633199

# A semiparametric random effects model for multivariate competing risks data

Thomas H. Scheike
Department of Biostatistics, University of Copenhagen, Øster Farimagsgade 5, Copenhagen DK-1014, Denmark, ; kd.uk.tatsoib@st
Yanqing Sun
Department of Mathematics and Statistics, The University of North Carolina at Charlotte, 9201 University City Boulevard, Charlotte, North Carolina 28223, U.S.A., ude.ccnu@nusay
Mei-Jie Zhang
Division of Biostatistics, Medical College of Wisconsin, 8701 Watertown Plank Road, Milwaukee, Wisconsin 53226, U.S.A., ude.wcm@eijiem

## Abstract

We propose a semiparametric random effects model for multivariate competing risks data when the failures of a particular type are of interest. Under this model, the marginal cumulative incidence functions follow a generalized semiparametric additive model. The associations between the cause-specific failure times can be studied through dependence parameters of copula functions that are allowed to depend on cluster-level covariates. A cross-odds ratio-type measure is proposed to describe the associations between cause-specific failure times, and its relationship to the dependence parameters is explored. We develop a two-stage estimation procedure where the marginal models are estimated in the first stage and the dependence parameters are estimated in the second stage. The large sample properties of the proposed estimators are derived. The proposed procedures are applied to Danish twin data to model the cumulative incidence for the age of natural menopause and to investigate the association in the onset of natural menopause between monozygotic and dizygotic twins.

Some key words: Binomial modelling, Copula function, Cross-odds ratio, Cumulative incidence function, Danish twin data, Estimating equation, Inverse-censoring probability weighting, Two-stage estimation

## 1. Introduction

This work is motivated by a dataset from the Danish twin registry (Skytthe et al., 2003) on the timing of natural menopause among female Danish twins, where natural menopause and surgical menopause are two competing events and time to menopause may be correlated within each pair of twins.

Much of the early work for analyzing competing risks data has been focused on estimation and modelling the cause-specific hazard function (Kalbfleisch & Prentice, 2002). The cumulative incidence function by contrast measures the absolute cause-specific risk, that is, the probability of having experienced failure from a particular cause over time. The effect of a covariate on the cause-specific hazard function and cumulative incidence function can be very different. The two measures provide different perspectives for cause-specific failure times; see Fine & Gray (1999), Katsahian et al. (2006), Latouche et al. (2007) and Beyersmann & Schumacher (2007).

Gray (1988) and Fine & Gray (1999) proposed a proportional hazards model for the subdistribution hazard to model the cumulative incidence function directly. A general semiparametric model for the cumulative incidence function is studied by Scheike et al. (2008), where some covariates have time-varying effects and others have constant effects. Scheike & Zhang (2008) proposed some goodness-of-fit procedures for model checking for these semiparametric models. Here we extend this work to deal with clusters by providing robust standard errors, and thereby also extend the recent work by Chen et al. (2008) that considered the k-sample case.

Multivariate failure time data analysis has received much attention in the literature. Dabrowska (1988) proposed a nonparametric estimator for the multivariate survival function. Gill (1993) proposed the concept of the iterated odds ratio to model dependence between failure times. The marginal model approach of Wei et al. (1989) fits univariate models to the marginal distributions without modelling the within-cluster dependence. The efficiency of the marginal approach can be improved by introducing weights that account for failure time dependencies (Cai & Prentice, 1995). A frailty model approach, on the other hand, attempts to model this dependence by postulating a cluster-specific frailty to account for each cluster’s particular susceptibility to failures; see Clayton (1978), Oakes (1982, 1989), Clayton & Cuzick (1985) and Glidden (2000), among others.

There has been little work on developing statistical methods for correlated cause-specific failure times in the competing risks setting. Recently, Bandeen-Roche & Ning (2008) and Bandeen-Roche & Liang (2002) introduced a nonparametric cause-specific cross hazard ratio association measure for multivariate cause-specific failure times; see also Cheng & Fine (2008). Cheng et al. (2007) studied nonparametric estimation of bivariate cause-specific hazard functions and bivariate cumulative incidence functions. They also proposed two nonparametric association measures in terms of the bivariate cause-specific hazard functions and bivariate cumulative incidence functions. These methods do not deal with covariates and require specification of joint cause-specific intensities for all cause combinations. It may be difficult to apply these measures to cause-specific failure time regressions.

To analyze clustered competing risks data, Katsahian et al. (2006) proposed a frailty model for the subdistribution hazard by multiplying lognormal random cluster effects to the proportional subdistribution hazard model. This model allows subject-specific regression effects. One drawback of this approach is that the frailty parameter can be identified solely from the marginal models and that the marginal models depend on the distribution of the cluster-specific effects.

We suggest an alternative approach to avoid this problem by working on a new semiparametric random effects model. Under this model, the marginal cumulative incidence functions follow a generalized semiparametric additive model that is independent of the frailty parameters and frailties account only for correlations between competing risks failure times within clusters. The frailty variance will thus reflect solely the amount of variation due to clusters. The concept of cross-odds ratio is introduced to measure the degree of association between two cause-specific failure times within the same cluster.

## 2. Estimation procedures with multivariate competing risks data

### 2.1. Model specification

Let K be the number of clusters and nk the number of subjects within the kth cluster. The failure time and cause of failure for the ith subject in the kth cluster are denoted by Tki and ki, respectively. Assume that cause ki = 1 is the primary cause of interest and ki = 2 for other competing causes. Let Xki = (1, Xki,1, . . . , Xki,p)T and Zki = (Zki,1, . . . , Zki,q)T be the associated covariates. Let Cki be the censoring time. The observed competing risks data can be represented by {(ki, ∊̃ki, Xki, Zki)} for i = 1, . . . , nk and k = 1, . . . , K, where ki = min(Tki, Cki), Δki = (Tki Cki) is an indicator for uncensored failure time and ∊̃ki = ki Δki. The clusters are assumed to be independent but the competing risks failure times within clusters are generally correlated. Let [0, τ] be the time period from which data are collected. We assume that (Tki, ki) are independent of Cki given the covariates (Xki, Zki).

Let θk be independently distributed random variables. Assume that the subjects within clusters are independent conditional on θk and (Xki, Zki). Let $P1*(t|xki,zki,θk)=pr(Tki≤t,∊ki=1|θk,xki,zki)$ be the cluster-specific conditional cumulative incidence function given θk and (Xki, Zki) = (xki, zki). We consider the following random effects model:

$P1*(t|xki,zki,θk)=1-exp(-θkΨνk-1[exp{-η(t)Txki-(γTzki)t}]),$
(1)

where η(t) is a (p + 1)-dimensional vector of regression functions, γ is a q-dimensional vector of parameters, Ψνk(t) = Eθk{exp(−θkt) | xki, zki} is the Laplace transform of the random effects θk conditional on (xki, zki) and νk are the dependence parameters. For example, here we take Ψνk(t) as the Laplace transform of a gamma distribution with mean 1 and variance νk, where νk can depend on some cluster-level covariates Qk that may be a subset of (xki, zki). This generalization allows for different degrees of association among different types of clusters. For example, the associations between the times to natural menopause for monozygotic twins may be different from that of dizygotic twins. Here, we consider the following simple regression model:

$νk=αTQk,$
(2)

where α is an r-dimensional vector of parameters and Qk is a regression design vector.

Under (1), the marginal cumulative incidence function, P1(t | xki, zki) = pr (Tki t, ki = 1 | xki, zki), follows the generalized semiparametric additive model:

$-log{1-P1(t|xki,zki)}=η(t)Txki+(γTzki)t.$
(3)

The joint survival function of the cause-specific failure times has an Archimedean copula function representation,

$E[{1−I(Tki≤t,∊ki=1)}{1−I(Tkj≤s,∊kj=1)}|xki,zki,xkj,zkj] =Ψνk[Ψνk−1{1−P1(t|xki,zki)}+Ψνk−1{1−P1(s|xkj,zkj)}] =Ck{1−P1(t|xki,zki),1−P1(s|xkj,zkj)},$

where $Ck(u,v)=Ψνk{Ψνk−1(u)+Ψνk−1(v)}$. Various frailty distributions for θk can be used to describe correlations among members of each cluster. When the conditional distribution of θk given (xki, zki) follows the gamma distribution with mean 1 and variance of νk, Ck(u, v) = (uνk + vνk − 1)−1/νk is the Clayton copula; see Clayton (1978). Other Archimedean copula families include those of Frank (1979) and Gumbel (1960), which are generated by the Laplace transforms of the log series distribution and positive stable distribution; see also Hougaard (2000).

### 2.2. Estimation of marginal models

The random effects model (1) implies the semiparametric additive model (3) for the marginal cumulative incidence functions. This suggests that model (1) can be estimated through a two-stage procedure, where the parameters of the marginal models are estimated in the first stage and the dependence parameter is estimated in the second stage. Ignoring the dependence between the failure times, the parameters of the marginal models can be estimated by modifying the procedure of Scheike et al. (2008) for univariate competing risks data. This resembles the marginal approach of Wei et al. (1989).

Let Gc(t | Xki, Zki) = pr(Cki > t | Xki, Zki) and Nki (t) = (Tki t, ki = 1) be the counting process associated with cause 1. Under conditional independence between Cki and (Tki, ki) given the covariates, we have

$E{ΔkiNki(t)Gc(Tki|Xki,Zki)|Xki,Zki}=P1(t|Xki,Zki).$

We consider estimating equations based on ΔkiNki (t) / Gc(Tki | Xki, Zki). Although in practice the censoring distribution Gc(t | xki, zki) is unknown, it can be estimated by the Kaplan–Meier estimator or by using a regression model to improve the efficiency. For simplicity, we use the Kaplan–Meier estimator Ĝc(t). In the regression case, Gc(t | xki, zki) will be estimated based on the data {(ki, Δki, Xki, Zki), i = 1, . . . , nk, k = 1, . . . , K}.

Let $P1(k)(t,η,γ)$ be the nk × 1 vector of P1(t | xki, zki) for i = 1, . . . , nk, R(k)(t) be the nk × 1 vector of adjusted responses ΔkiNki (t)/Ĝc(Tki), $Dη(k)(t,η,γ)$ and $Dγ(k)(t,η,γ)$ be the nk × (p + 1) and nk × q matrices with the ith rows equal to Dη,ki (t, η, γ) = P1(t | xki, zki)/η(t) and Dγ,ki (t, η, γ) = P1(t | xki, zki)/γ, respectively. Let W(k)(t) = diag{wki(t)} where wki (t) are the weights.

Following Scheike et al. (2008), the marginal regression functions η(t) and regression parameters γ can be estimated based on the following estimating equations:

$Uη(t,η,γ)=DηT(t,η,γ)W(t){R(t)−P1(t,η,γ)}=0,$
(4)

$Uγ(τ,η,γ)=∫0τDγT(t,η,γ)W(t){R(t)−P1(t,η,γ)}dt=0,$
(5)

where W(t) = diag{W(k)(t)}, R(t) = [{R(1)(t)}T, . . . , {R(K)(t)}T]T, and

$P1(t,η,γ)={{P1(1)(t,η,γ)}T,…,{P1(K)(t,η,γ)}T}T,Dη(t,η,γ)=[{Dη(1)(t,η,γ)}T,…,{Dη(K)(t,η,γ)}T]T,Dγ(t,η,γ)=[{Dγ(1)(t,η,γ)}T,…,{Dγ(K)(t,η,γ)}T]T.$

The estimating equations (4) and (5) can be solved by an iterative procedure similar to that given in Scheike et al. (2008); see the Appendix. The procedure can be implemented using standard software for generalized linear models. However, the standard errors obtained from such software are conservative. Correct estimation of the standard errors is described in the following.

Suppose that there exists a constant N0 such that nk < N0 for all k. Under regularity conditions similar to Scheike et al. (2008) we show in the Appendix that as K → ∞, $K12(γ^−γ)$ and $K12{η^(t)−η(t)}$ are asymptotically equivalent to the following decompositions

$K12Iγ−1∑k=1K∑i=1nkWγ,ki(τ), K12{Iη(t)}−1∑k=1K∑i=1nkWη,ki(t),$

respectively, for which the formula for Iγ, Iη(t), Wγ,ki (τ) and Wη,ki (t) are given in the Appendix. The sums $∑i=1nkWγ,ki(τ)$ and $∑i=1nkWη,ki(t)$ are independent and identically distributed with mean zero for k = 1, . . . , K. Thus $K12(γ^−γ)$ converges in distribution to a mean zero normal random vector by a simple application of the central limit theorem. By an application of Theorem 19.5 of van der Vaart (1998), $K12{η^(t)−η(t)}$ converges weakly to a mean zero Gaussian process on t [0, τ].

Let Îγ, Ŵγ,ki(τ), Îη(t) and Ŵη,ki (t) be the estimators of their population counterparts with the plug-in estimates of all needed quantities. The asymptotic variance of $K12(γ^−γ)$ and $K12{η^(t)−η(t)}$ can be estimated by

$Σ^γ=KI^γ−1[∑k=1K{∑i=1nkW^γ,ki(τ)}⊗2]I^γ−1,Σ^η(t)=K{I^η(t)}−1[∑k=1K{∑i=1nkW^η,ki(t)}⊗2]{I^η(t)}−1,$

where a2 = aaT for a column vector a. These estimators can be used to construct confidence intervals for γ and η(t).

The Gaussian multiplier resampling method can be used to construct simultaneous confidence intervals for η(t) over a fixed time interval and make inferences. By application of empirical process theory (van der Vaart, 1998), $K12{η^(t)−η(t)}$ converges weakly to the same mean zero Gaussian limit processes as $K12{I^η(t)}−1∑k=1K∑i=1nkGkW^η,ki(t)$, given the observed data, where G1, . . . , GK are independent standard normal random variables. We refer to Lemma 1 of Sun & Wu (2005) for a result that applies directly to this setting. Hence a (1 − α) × 100% asymptotic confidence band for (t) can be constructed by repeated resampling of G1, . . . , GK.

### 2.3. Estimation of dependence parameters

Let

$Vkij(t)={Gc(Tki)}−1ΔkiNki(t){Gc(Tkj)}−1ΔkjNkj(t),$
(6)

where i, j are two subjects within cluster k. It follows that, for ij,

$vkij(t,α,η,γ)=E{Vkij(t)|Xki,Zki,Xkj,Zkj} =1−exp{−Λki(t)}−exp{−Λkj(t)}+Ψνk(Ψνk−1[exp{−Λki(t)}]+Ψνk−1[exp{−Λkj(t)}]),$
(7)

where νk = αTQk under (2) and Λki (t) = η(t)TXki + (γTZki) t.

Let Hk denote the index set for the kth cluster. Replacing Gc(t) by Ĝc(t) in (6) to obtain kij(t), we consider the following estimating equation for α:

$Uα(α,η^,γ^,G^c)=∫0τ∑k=1K∑i,j∈Hk,i≠jDα,kij(t,α,η^,γ^){V^kij(t)−vkij(t,α,η^,γ^)}dt=0,$
(8)

where Dα,kij(t, α, , ) = vkij (t, α, , )/α, and (t), are the estimators of η(t), γ from the marginal model. Let be the estimator of α solving (8). Then the cluster-specific dependence parameter is given by k = TQk. By Taylor approximation,

$K12(α^−α)=−{K−1∂Uα(α,η,γ,Gc)∂α}−1K−12Uα(α,η^,γ^,G^c)+op(1).$

We show in the Appendix that $K12(α^−α)$ has the following approximation:

$K12(α^−α)=K−12{K−1ℐα(α,η,γ)}−1{∑k=1KWα,k(α,η,γ,Gc)}+op(1),$

for which the formulae for Wα,k(α, η, γ, Gc) and α(α, η, γ) are given in the Appendix. It follows that $K12(α^−α)$ is asymptotically normal with mean zero and its variance can be estimated by

$Σ^α=Kℐ^α−1{∑k=1K(W^α,k)⊗2}ℐ^α−1,$

where Îα and Ŵα,k are the empirical counterparts of α (α, η, γ) and Wα,k(α, η, γ, Gc), respectively.

### 2.4. Cross-odds ratio

There seems to be no simple direct extension of the association measures for multivariate failure times to the competing risks data. The popular Kendall’s tau-like interpretation of dependence parameter is no longer applicable to multivariate cause-specific failure times since cumulative incidence functions are in general not distribution functions.

We introduce a cross-odds ratio to describe the dependence between cause-specific failure times. Consider the case of bivariate failure times (Tki, Tkj) with causes (ki, kj). We define the following association measure for the failure times due to cause 1 between individual i and j in the same cluster k:

$πi|j(t)=pr(Tki≤t,∊ki=1|Tkj≤t,∊kj=1)/pr{(Tki≤t,∊ki=1)c|Tkj≤t,∊kj=1}pr(Tki≤t,∊ki=1)/pr{(Tki≤t,∊ki=1)c},$
(9)

where (A)c means the event complementary to A. This measure of association resembles the odds ratio for the event probabilities and it may vary with time. If the events (Tki t, ki = 1) and (Tkj t, kj = 1) are independent, then πi|j (t) = 1. If they are positively associated with pr(Tki t, ki = 1 | Tkj t, kj = 1) > pr(Tki t, ki = 1), then πi|j (t) > 1; otherwise πi|j (t) < 1. If pr(Tki t, ki = 1 | Tkj t, kj = 1) = 1, thus perfectly positive correlation, then πi|j (t) = ∞. If pr(Tki t, ki = 1 | Tkj t, kj = 1) = 0, then πi|j (t) = 0. Further, this measure is symmetric about the individuals i and j provided that the two have the same marginal subprobabilities, i.e., pr(Tki t, ki = 1) = pr(Tkj t, kj = 1).

In the presence of covariates, the conditional cross-odds ratio given the covariates, πi|j (t | xki, zki, xkj, zkj), can be similarly defined as in (9) by replacing the probabilities with the conditional probabilities given the covariates. An important advantage of our definition of the cross-odds ratio is that it is readily available once model (1) is estimated.

When the random effects θk follow a gamma distribution with mean 1 and variance νk, Ψνk(x) = (1 +νkx)−1/νk, for x > −1/νk and its inverse function is $Ψνk−1(y)=(y−νk−1)/νk$ for y > 0. Once the marginal cumulative incidence model and the frailty variance are estimated the cross-odds ratio πi|j is simple to estimate. The cross-odds ratio is closely related to the dependence parameter νk of the Clayton Archimedean copula family. Indeed, when all the covariates are cluster level covariates, i.e. (xki, zki) = (xkj, zkj), πi|j (t | xki, zki, xkj, zkj) = 2 for νk = 1, and πi|j (t | xki, zki, xkj, zkj) ≈ νk + 1 for moderate values of νk in the range of [0, 2]. The association modelled by the Clayton copula is always positive. A larger value of νk indicates stronger positive association in terms of times to cause-specific failure times within clusters. Other copula functions may be used to model negative as well as positive associations.

## 3. Danish twin data

We identified 503 pairs of female twins born between 1931 and 1952 through the Danish Twin Registry. There are 269 monozygotic twin pairs and 234 dizygotic twin pairs identified by Skytthe et al. (2003). Before they reach natural menopause, some women may experience surgical menopause as a result of removal of the uterus, cervix or ovaries and others may be treated with hormones that affect the timing of menopause. Surgical menopause and hormonal treatment are therefore competing risks for natural menopause. Of 1006 twins, 416 twins experienced natural menopause, 57 twins had surgical related menopause, 183 twins underwent hormonal treatments before natural menopause, 246 twins had censored event times. There were 104 twins with unknown failure types or censoring status and they were excluded from the analysis. These summaries do not reflect that the probabilities of different events changed considerably over the time cohorts. We partition the data into four cohorts of almost equal size according to birthdate. Cohort 1 from 1931–1936 has 120 twin pairs, cohort 2 from 1937–1942 has 123 twin pairs, cohort 3 from 1943–1947 has 121 twin pairs, and cohort 4 from 1948–1952 has 139 twin pairs. In cohort 1, there was zero censoring, 154 natural menopause, 22 surgical menopause, 40 hormonal treatments, and 24 cases with unknown failure types or censoring status. These numbers were (3, 137, 13, 64, 29) for cohort 2, (65, 89, 11, 54, 23) for cohort 3, and (178, 36, 11, 25, 28) for cohort 4. In the analysis, we allow the censoring distribution to depend on cohorts.

We first estimate the marginal cumulative incidence of natural menopause excluding the cohort structure. Let xki = {1, Iki (dizygotic twin)}T, where Iki (dizygotic twin) is 1 if the i th twin from the kth family is dizygotic for i = 1, 2 and k = 1, . . . , 503. Under model (1) without covariate zki, the marginal cumulative incidence is − log{1 − P1(t | xki)} = η(t)Txki, which leads to the nonparametric estimates for the distributions of the ages to natural menopause for the monozygotic and dizygotic twins. The estimates, 95% confidence intervals and 95% confidence bands based on the Gaussian multiplier resampling technique are shown in Fig. 1(a) for the monozygotic twins. The plots for the dizygotic twins are similar and omitted to save space.

Estimated probability of natural menopause in Danish twins data. The estimated marginal cumulative incidence functions for monozygotic twins are shown in panel (a) and the upper curves in panel (b) are for cohort 1 and the lower curves are for cohort ...

To examine whether the cumulative incidence of natural menopause changes over the different cohorts, we estimate the cumulative incidence function for natural menopause by adjusting for time cohort effects. Let Iki (cohort 2) be the cohort indicator of whether the twins in the kth family are born between 1937 and 1942. In addition to xki defined above, we introduce cohort covariates zki to the random effects model (1), where zki = {Iki (cohort 2), Iki (cohort 3), Iki (cohort 4)}T. The marginal cumulative incidence functions then follow the semiparametric additive model − log{1 − P1(t | xki, zki)} = η(t)Txki + (γTzki) t. This model assumes separate nonparametric baseline functions for monozygotic and dizygotic twins and constant time cohort effects. The goodness-of-fit procedures of Scheike & Zhang (2008) indicate that the cohorts can well be described as having constant effects. Figure 1 (b) shows the estimated cumulative incidence functions of natural menopause for cohort 1 and 4 for monozygotic twins. There is a clear trend over time in that the later cohort experience their natural menopause later. For cohort 4 we only estimate the cumulative incidence curve for the ages between 40 and 55 due to censoring. The estimated cohort effects are (−1.03, −3.88, −5.72) × 10−2 with the standard errors of (0.65, 0.58, 0.53) × 10−2, respectively.

We now examine whether there is an association in the occurrence of natural menopause between twins of the same family under the random effects model (1). To allow for different degrees of association for different types of twins, we let νk = α0 + α1xk under (2), where xk is 0 for monozygotic twins and 1 for dizygotic twins. Let ρ1 = α0 and ρ2 = α0 + α1. Using the estimation procedure presented in § 2.3, we get 1 = 1.04 with the standard error of 0.42, and 2 = 0.39 with the standard error of 0.29. Hence, a significant association exists in the monozygotic twins but not for the dizygotic twins.

The joint probability of both twins experiencing natural menopause by age t depends on the degree of association between the twins. This probability is vkij(t) given in (7) under the random effects model (1). For monozygotic twins in cohort 1, this joint probability is plotted in Fig. 1 (c) against age. The random effect θk in model (1) represents family-to-family variability and induces dependence between twins. Let θk,0.95 and θk,0.05 be the 95th percentile and the 5th percentile of the gamma distribution with mean 1 and variance νk, respectively. The joint probabilities of both twins experiencing natural menopause by age t in a twin pair with θk = θk,0.95 and θk = θk,0.05 can be expressed as $P1*(t|xki,zki,θk,0.95)P1*(t|xkj,zkj,θk,0.95)$ and $P1*(t|xki,zki,θk,0.05)P1*(t|xkj,zkj,θk,0.05)$, respectively. For monozygotic twins in cohort 1, these joint probabilities are given in Fig. 1 (c). The plot of the joint probability of both twins experiencing natural menopause by age t, P1(t | xki, zki)P1(t | xkj, zkj), under the marginal model (3) and the presumed independence between twins is also given in Fig. 1 (c). Similar plots for dizygotic twins are given in Fig. 1 (d). Other quantities such as the conditional probabilities of the second twin experiencing natural menopause given that the first twin has already experienced her natural menopause can also be calculated and plotted.

Finally, we calculate the cross-odds ratio to show how the onset of natural menopause for one twin affects the probability of the onset for the other twin. The cross-odds ratios for monozygotic twins are around 2. This indicates that, if one twin has experienced natural menopause, the odds of the second twin experiencing natural menopause by the same age essentially doubles the unconditional odds of natural menopause for the second twin. For the dizygotic twins, we find a cross-odds ratio that lies around 1.4.

## 4. Simulation study

We conducted a small simulation study for model (1) to evaluate small sample performance of the proposed method. We take K = 100 and 200 pairs and let the random effects follow a Gamma distribution with mean 1 and variance ν. We consider a situation with one covariate that is uniformly distributed on [0, 1]. The baseline function is set as η(t) = 0.5t and the regression coefficient γ = 0.5. The censoring times are generated from the uniform distribution on [0, 2]. The probability of observing the cause of interest is 0.42, which is comparable to the Danish twin dataset where natural menopause occurred in about 40% of the twins.

Table 1 summarizes the simulation results for and based on 1000 realizations for γ = 0.5 and ν equal to 2, 1 and 0.5. The estimator performs satisfactorily at all levels of the random effect variance ν and for both sample sizes. The estimator is essentially unbiased and its variance is well estimated by the robust variance estimator. The biases of the estimator for ν are reasonably small, and the estimator for the variance of perform well. The simulations also find the baseline function well approximated by its estimator, which are not reported here.

Summaries of the simulation results based on 1000 realizations for K twin pairs

## 5. Discussion

The random effects model (1) is estimated through a two-stage approach where the marginal parameters are estimated separately from the dependence parameters. Additional efficiency can be gained by jointly estimating both the marginal parameters and the dependence parameters and by using all the joint observations of the inverse probability of censoring weighted responses for the dependence parameters. This issue deserves further attention and investigation.

The asymptotic distribution of i|j (t | xki, zki, xkj, zkj) can be derived using the delta method. Confidence intervals and hypothesis testing procedures about the cross-odds ratio can be derived accordingly. The cross-odds ratio may also be extended in different ways to show how the association persists across the time-scale.

## Acknowledgments

The research was partially supported by grants from the National Cancer Institute to Scheike and Zhang, and by grants from the National Science Foundation and the National Institutes of Health to Sun. Thomas Scheike was also supported by a grant from the Danish Research Council. We appreciate the valuable comments of the reviewers and the editor, which have improved the content and presentation of this paper.

## Appendix

#### An iterative procedure for the estimators and (t)

The estimating equations (4) and (5) can be solved simultaneously in a way similar to Scheike et al. (2008). By the Taylor expansion of P1(t, η, γ) around the current value of the estimates ((m)(t), (m)), we have

$P1(t,η,γ)=P1(t,η^(m),γ^(m))+Dγ(t,η^(m),γ^(m)){γ−γ^(m)}+Dη(t,η^(m),γ^(m)){η(t)−η^(m)(t)}+op(K−12).$

Let $Dη(m)(t)=Dη(t,η^(m),γ^(m))$, $Dγ(m)(t)=Dγ(t,η^(m),γ^(m))$ and $P1(m)(t)=P1(t,η^(m),γ^(m))$. We have the following approximated estimating equations:

$Dη(m)(t)TW(t)[R(t)−P1(m)(t)−Dη(m)(t){η(t)−η^(m)(t)}−Dγ(m)(t){γ−γ^(m)}]=0,$
(A1)

$∫0τDγ(m)(t)TW(t)[R(t)−P1(m)(t)−Dη(m)(t){η(t)−η^(m)(t)}−Dγ(m)(t){γ−γ^(m)}]dt=0.$
(A2)

Solving the equations (A1) and (A2) for γ and η(t), the (m + 1)th iterative estimators are given by $γ^(m+1)=γ^(m)+{Iγ(m)}−1Bγ(m)$ and $η^(m+1)(t)=η^(m)(t)+{Iη(m)(t)}−1{Dη(m)(t)}TW(t){R(t)−P1(m)(t)−Dγ(m)(t){Iγ(m)}−1Bγ(m)}$, where

$Iγ(m)=∫0τDγ(m)(t)TW(t)H(m)(t)Dγ(m)(t)dt, Bγ(m)=∫0τDγ(m)(t)TW(t)H(m)(t){R(t)−P1(m)(t)}dt,H(m)(t)=I−Dη(m)(t){Iη(m)(t)}−1Dη(m)(t)TW(t),Iη(m)(t)=Dη(m)(t)TW(t)Dη(m)(t).$
(A3)

#### On the asymptotic distributions of $K12(γ^−γ)$ and $K12{η^(t)−η(t)}$

Let Dη(t), Dγ(t) and P1(t, η, γ) be the corresponding terms of $Dη(m)(t)$, $Dγ(m)(t)$ and $P1(m)(t)$ evaluated at the true values of {η(t), γ}. Similarly, Iγ, Iη(t), Bγ and H(t) are the corresponding terms for (A3). By applying the Taylor approximation,

$K12(γ^−γ)=K12Iγ−1Bγ+op(1).$

Let $Q(t)=DγT(t)W(t)Dη(t){Iη(t)}−1$, $ζγ,ki(t)={Dγ,ki(t)−Q(t)Dη,ki(t)}wki(t)$, $B˜γ=∑k,i∫0τζγ,ki(t)dt$, {Gc(Tki)−1kiNki (t) − P1(t, xki, zki)}, and $Δγ=∑k,i∫0τ{Dγ,ki(t)−Q(t)Dη,ki(t)}wki(t)ΔkiNki(t){G^c(Tki)Gc(Tki)}−1{Gc(Tki)−G^c(Tki)}dt$. Then Bγ = γ + Δγ. Note that

$ΔkiNki(t)G^c(Tki)−Gc(Tki)G^c(Tki)Gc(Tki)=−ΔkiNki(t)Gc(Tki)∑l=1K∑j=1nl∫0τℐ(s≤Tki)dMljC(s)Y(s)+op(K−12),$

where Y(s) = ∑k,iYki (s) = ∑k,i(ki s) and $MljC(s)=ℐ(T˜lj≤s,Δlj=0)−∫0sYlj(u)d{−logGc(u)}$. We have $Δγ=∑k,i∫0τψγ,ki(t)dt+op(K12)$, where $ψγ,ki(t)=∫0τ{Y(s)}−1qγ(s,t)dMkiC(s)$ and qγ(s, t) = ∑l,j {Gc(Tlj)}−1{Dγ,lj (t) − Q(t)Dη,lj}wlj (tljNlj(t)(s Tlj t). Thus,

$K12(γ^−γ)=K12Iγ−1∑k,iWγ,ki(τ)+op(1),$
(A4)

where $Wγ,ki(τ)=∑k,i∫0τ{ζγ,ki(t)+ψγ,ki(t)}dt$. Let Îγ and Ŵγ,ki (τ) be estimators of Iγ and Wγ,ki (τ) obtained by inserting the estimators (t), and $M^kiC(t)$. The asymptotic variance of $K12(γ^−γ)$ can be consistently estimated by $Σ^γ=KI^γ−1[∑k=1K{∑i=1nkW^γ,ki(τ)}⊗2]I^γ−1$.

Similarly, by using a Taylor expansion, we have

$K12{η^(t)−η(t)}=K12{Iη(t)}−1∑k,iWη,ki(t)+op(1),$
(A5)

where

$Wη,ki(t)={ζη,ki(t)+ψη,ki(t)−DηT(t)W(t)Dγ(t)Iγ−1Wγ,ki(τ)},ζη,ki(t)=Dη,ki(t)wki(t){ΔkiNki(t)Gc(Tki)−P1(t,xki,zki)},ψη,i(t)=∫0τ{∑l,jDη,lj(t)wlj(t)ΔljNlj(t)Gc(Tlj)ℐ(s≤Tlj≤t)}dMkiC(s)y(s).$

Both Iη(t) and Wη,ki (t) can be consistently estimated by the plug-in estimators, denoted as Îη(t) and Ŵη,ki (t), and the variance of $K12{η^(t)−η(t)}$ can be estimated by $Σ^η(t)=K{I^η(t)}−1[∑k=1K{∑i=1nkW^η,ki(t)}⊗2]{I^η(t)}−1$.

#### On the asymptotic distribution of $K12(α^−α)$

To establish the asymptotic distribution for $K12(α^−α)$, we consider the decomposition

$Uα(α,η^,γ^,G^c)=∫0τ∑k=1K∑i,j∈Hk,i≠jDα,kij(t,α,η^,γ^){Vkij(t)−vkij(t,α,η,γ)}dt+∫0τ∑k=1K∑i,j∈Hk,i≠jDα,kij(t,α,η^,γ^){V^kij(t)−Vkij(t)}dt−∫0τ∑k=1K∑i,j∈Hk,i≠jDα,kij(t,νk,η^,γ^){vkij(t,α,η^,γ^)−vkij(t,α,η,γ)}dt.$
(A6)

Since ${G^c(Tki)}−1{G^c(Tki)−Gc(Tki)}=−∑l=1K∑m=1nl∫0τℐ(s≤Tki){Y(s)}−1dMlmC(s)+op(K−12)$, the second term of Uα(α, , , Ĝc) in (A6) equals

$2∫0τ∑k=1K∑i,j∈Hk,i≠jDα,kij(t,α,η,γ)Vkij(t)∑l=1K∑m=1nl∫0τℐ(s≤Tki)dMlmC(s)Y(s)dt+op(K12) =∑l=1K∑m=1nl∫0τ∫0τqα(s,t)dMlmC(s)dt+op(K12),$
(A7)

where $qα(s,t)=2{∑k=1K∑i,j∈Hk,i≠jDα,kij(t,α,η,γ)Vkij(t)ℐ(s≤Tki)}{Y(s)}−1$.

By a Taylor expansion at the true values of {η(t), γ}, the third term of (A6) equals

$K∫0τqx(t,α,η,γ){η^(t)−η(t)}dt+K∫0τqz(t,α,η,γ)(γ^−γ)dt+op(K12),$
(A8)

where

$qx(t,α,η,γ)=K−1∑k=1K∑i,j∈Hk,i≠jDα,kij(t,α,η,γ){∂vkij∂Λki(t)xkiT+∂vkij∂Λkj(t)xkjT},qz(t,α,η,γ)=K−1∑k=1K∑i,j∈Hk,i≠jDα,kij(t,α,η,γ){∂vkij∂Λki(t)zkiT+∂vkij∂Λkj(t)zkjT}t.$

By a Taylor approximation,

$K12(α^−α)=−{K−1∂Uα(α,η,γ,Gc)∂α}−1K−1/2Uα(α,η^,γ^,G^c)+op(1).$

Let $ℐα(α,η,γ)=∑k=1K∑i,j∈Hk,i≠j∫0τ{∂vkij(t,α,η,γ)/∂α}⊗2dt$. It follows from equations (A4), (A5), (A6), (A7) and (A8) that

$K12(α^−α)=K−1/2{K−1ℐα(α,η,γ)}−1∑k=1KWα,k(α,η,γ,Gc)+op(1),$

where

$Wα,k(α,η,γ,Gc)=∑i,j∈Hk,i≠j∫0τDα,kij(t,νk,η,γ){Vkij(t)−vkij(t,α,η,γ)}dt+∑i=1nk∫0τ∫0τqα(s,t)dMikC(s)dt−∑i=1nk∫0τqx(t,α,η,γ)nk{Iη(t)}−1Wη,ki(t)dt−∑i=1nk∫0τqz(t,α,η,γ)nkIγ−1Wγ,ki(τ)dt,$

and Iη(t), Wη,ki (t), Iγ and Wγ,ki (τ) are given following (A3), (A4) and (A5). By the central limit theorem $K12(α^−α)$ is asymptotically normal with mean zero and a variance that can be estimated by $Σ^α=Kℐ^α−1∑k=1K{W^α,k}⊗2ℐ^α−1$, where α and Ŵα,k are the empirical counterparts of α (α, η, γ) and Wα,k (α, η, γ, Gc), respectively, obtained by inserting the estimators of α, η(t), γ and by replacing $M^kiC(t)$ for $MkiC(t)$.

## References

• Bandeen-Roche K, Liang K-Y. Modelling multivariate failure time associations in the presence of a competing risk. Biometrika. 2002;89:299–314.
• Bandeen-Roche K, Ning J. Nonparametric estimation of bivariate failure time associations in the presence of a competing risk. Biometrika. 2008;95:221–32. [PubMed]
• Beyersmann J, Schumacher M. Letter to the editor: misspecified regression model for the subdistribution hazard of a competing risk. Statist Med. 2007;26:1649–52. [PubMed]
• Cai J, Prentice RL. Estimating equations for hazard ratio parameters based on correlated failure time data. Biometrika. 1995;82:151–64.
• Chen BE, Kramer JL, Greene MH, Rosenberg PS. Competing risks analysis of correlated failure time data. Biometrics. 2008;64:172–79. [PubMed]
• Cheng Y, Fine J. Nonparametric estimation of cause-specific cross hazard ratio with bivariate competing risks data. Biometrika. 2008;95:233–40.
• Cheng Y, Fine JP, Kosorok MR. Nonparametric association analysis of bivariate competing-risks data. J Am Statist Assoc. 2007;102:1407–15.
• Clayton DG. A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika. 1978;65:141–51.
• Clayton DG, Cuzick J. Multivariate generalisations of the proportional hazards model. J. R. Statist. Soc A. 1985;148:82–117.
• Dabrowska DM. Kaplan–Meier estimate on the plane. Ann Statist. 1988;16:1475–89.
• Fine JP, Gray RJ. A proportional hazards model for the subdistribution of a competing risk. J Am Statist Assoc. 1999;94:496–509.
• Frank MJ. On the simultaneous associativity of f(x, y) and x + y − f(x, y) Aequationes Mathematicae. 1979;19:194–226.
• Gill RD. Multivariate survival analysis. Theory Prob Appl. 1993;37:18–31.
• Glidden D. A two-stage estimator of the dependence parameter for the Clayton-Oakes model. Lifetime Data Anal. 2000;6:141–56. [PubMed]
• Gray RJ. A class of k-sample tests for comparing the cumulative incidence of a competing risk. Ann Statist. 1988;16:1141–54.
• Gumbel EJ. Bivariate exponential distributions. J Am Statist Assoc. 1960;55:698–707.
• Hougaard P. Analysis of Multivariate Survival Data. New York: Springer; 2000.
• Kalbfleisch JD, Prentice RL. The Statistical Analysis of Failure Time Data. 2nd ed. New York: Wiley; 2002.
• Katsahian S, Resche-Rigon M, Chevret S, Porcher R. Analysing multicentre competing risks data with a mixed proportional hazards model for the subdistribution. Statist Med. 2006;25:4267–78. [PubMed]
• Latouche A, Boisson V, Chevret S, Porcher R. Misspecified regression model for the subdistribution hazard of a competing risk. Statist Med. 2007;26:965–74. [PubMed]
• Oakes D. A model for association in bivariate survival data. J. R. Statist. Soc B. 1982;44:414–22.
• Oakes D. Bivariate survival models induced by frailties. J Am Statist Assoc. 1989;84:487–93.
• Scheike TH, Zhang M-J. Flexible competing risks regression modelling and goodness-of-fit. Lifetime Data Anal. 2008;14:464–83. [PubMed]
• Scheike TH, Zhang M-J, Gerds T. Predicting cumulative incidence probability by direct binomial regression. Biometrika. 2008;95:205–20.
• Skytthe A, Pedersen N, Kaprio J, Stazi M, Hjelmborg J, Iachine I, Vaupel J, Christensen K. Longevity studies in genomeutwin. Twin Res. 2003;6:448–54. [PubMed]
• Sun Y, Wu H. Semiparametric time-varying coefficients regression model for longitudinal data. Scand J Statist. 2005;32:21–47.
• van der Vaart AW. Asymptotic Statistics. Cambridge: Cambridge University Press; 1998.
• Wei LJ, Lin DY, Weissfeld L. Regression analysis of multivariate incomplete failure time data by modelling marginal distributions. J Am Statist Assoc. 1989;84:1065–73.

Articles from Biometrika are provided here courtesy of Oxford University Press

 PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers.