Piecewise Cause-Specific Association Analyses of Multivariate Untied or Tied Competing Risks Data

2.1 Multivariate competing risks data

In this paper, we consider clustered competing risks data. Suppose that in each cluster there are M+1 failure times (T0,T1,…,TM) and their corresponding cause indicators (∈0,∈1,…,∈M) from a mother and her M children, where M is random and is assumed to bear no information on the association. We only consider two causes here ∈j=1,2,j=0,…,M, as we are primarily interested in association in cause 1 events and can always group multiple competing events together. In addition, there may also be independent and non-informative censoring times (C0,C1,…,CM). Note that the data were collected retrospectively, hence it is reasonable to assume censoring times are independent of event times. This is in contrast with a perspective study where there may be some dependence among censoring times from mothers and children due to their inherent age differences. Hence one observes Yj=min(Tj,Cj), and ηj=∈jI{Tj≤Cj}, for j=0,1,…,M, where I is the indicator function. This notation incorporates two types of multivariate competing risks data. One is the sibship data, where pairwise exchangeability can be assumed among siblings, i.e. all siblings have the same marginals and the pairwise joint densities are symmetric. The other consists of multiple mother and child pairs where we do not impose the exchangeability assumption between a mother and her child, but assume exchangeability among all children. The observed data containing n i.i.d. clusters are denoted as {(Yi0,…,Yi,mi,ηi0,…,ηi,mi,mi), i=1,…,n}. The following two association measures are considered for each of the two types of multivariate competing risks data.

2.2 Bivariate cause-specific hazard ratio for clustered data and its plug-in estimator

Cheng and Fine [7] proposed a cause-specific association measure based on bivariate cause-specific hazard functions and demonstrated its equivalence to the conditional cross hazard ratio [4]. We now extend their association measure to the clustered competing risks data. To simplify the notation we focus on association analysis of cause 1 events, though the methods proposed here are applicable to other causes. Analogous to Cheng and Fine [7], for any child–mother pair (Tj,T0),1≤j≤M, we define a bivariate cause-specific hazard function for double cause 1 events from both subjects under the exchangeability assumption for all children

and define bivariate hazard functions of a single cause 1 event from the child or mother with the other member being at risk as the following:

Then the association measure for cause 1 events based on these bivariate hazard functions is

It quantifies the elevated risk that both a child and his/her mother failed at times (s,t), as compared to the product of the risks that only one of them failed at s or t, conditionally on both the child and the mother being at risk at (s,t).

A simple plug-in estimator of ζCM(s,t;1,1) can be constructed analogously to Cheng and Fine [7] based on all child–mother pairs. We first need to extend the event processes to the more complicated family structure with multiple children in a family. Let N11∗(s,t)=∑jMI(Yj≤s,ηj=1,Y0≤t,η0=1), N10*(s,t)=∑jMI(Yj≤s,ηj=1,Y0≥t), N01*(s,t)=∑jMI(Yj≥s,Y0≤t,η0=1), H*(s,t)=∑jMI(Yj≥s,Y0≥t), where M is the number of children in a family which may vary from family to family. Cheng et al. [10] divided the time scale as follows: 0=τ1,0<τ1,1<τ1,2<⋯<τ1,n1=τ1,0=τ2,0<τ2,1<τ2,2<⋯<τ2,n2=τ2, where τ1 and τ2 are the largest observed event times. These grids are predetermined based on practical interests as well as the feasibility of analyzing the data (i.e. the number of at risk or the number of the events should not be too small in any time grid). Assuming constant association measure ζCM(Ωqr) over the region Ωqr≡(τ1,q,τ1,q+1)×(τ2,r,τ2,r+1) with 0≤q≤n1−1,0≤r≤n2−1, we have

where E is the expectation operator and w(s,t) is some bounded deterministic function. The above equation is derived by recasting the hazard functions in eq. (1) in terms of the expectations of N11∗,N10∗,N01∗, and H∗, similar to the one given in Cheng and Fine [7]. That is, EN11∗(s,t)=E(M)P(Y1≤s,η1=1,Y0≤t,η0=1) and the expectations of the other three stochastic processes can also be expressed as the products of E(M) and the probabilities of a single event or at risk. Therefore, ζCM(Ωqr)=EN11∗(ds,dt)EH∗(s,t)EN10∗(ds,t)EN01∗(s,dt) and the above equation follows automatically.

Define the empirical process PnN11∗(s,t)=1n∑i=1n∑j=1miI(Yij≤s,Yi0≤t,ηij=1,ηi0=1), and similarly for PnH∗(s,t), PnN10∗(s,t), and PnN01∗(s,t). Since the number of children M is random, the above formula naturally handles different cluster sizes. Plugging the empirical processes into the above equation yields a plug-in estimator

where wˆ is a weight function satisfying ∫(s,t)∈Ωqrwˆ(s,t)dsdt=1, ||wˆ−w||∞→0, and n(wˆ−w)=OP(1). The selection of the weight function wˆ deserves some attention. If the association is constant within the region Ωqr, different weight functions usually would yield comparable estimates of ζCM. In this case, the uniform weight wˆ(s,t)=1(τ1,q+1−τ1,q)(τ2,r+1−τ2,r) is a convenient choice. However, if the association is time-varying within Ωqr, the weight will play an important role in summarizing the time-varying association. When wˆ(s,t)=1PnN10∗(ds,t)PnN01∗(s,dt)/∫Ωqr1PnN10∗(ds,t)PnN01∗(s,dt), we are averaging the point-wise association across Ωqr with equal weight. One may consider other weighting schemes such as putting more weights on early time points at which the number of at risk is large.

The consistency and asymptotic normality of ζˆCM(Ωqr) can be established using empirical processes techniques [15, 16]. Given EM2<∞ we can show that all the above indicator processes are P-Donsker. Since ζˆCM is a compact differential map of empirical processes indexed by Donsker classes, its consistency and asymptotic normality follow. In addition, there is an asymptotically linear representation n1/2(ζˆCM−ζCM)=n−1/2∑iIi+oP(1), where Ii are influence functions, upon which the variance of ζˆCM can be estimated explicitly. Detailed proofs are given in the section “Asymptotic properties of ζˆCM” in Appendix.

The association measure ζ(s,t;1,1) for bivariate competing risks data can be similarly adapted to sibship data denoted as ζCC(s,t;1,1). To estimate ζCC, we define the empirical processes PnN11#(s,t)=1n∑i=1n∑1≤j,j′≤miIYij≤s,ηij=1,Yij′≤t,ηij′=1 and similarly PnH#(s,t), PnN10#(s,t), and PnN01#(s,t) based on all within-cluster pairs. The asymptotic property of the plug-in estimator for ζCC can be established using similar arguments to those in the section “Asymptotic properties of ζˆCM” in Appendix except assuming EM4<∞ as a sufficient condition, which is reasonable as the family size M is random but bounded.

2.3 Cause-specific cross hazard ratio for clustered data

For mother and all children data, under the exchangeability assumption for the children, the child–mother cause-specific density was defined as

for k,l=1,2 and j=1,…,M. Assuming g is absolutely continuous, a child–mother pair (Tj,T0) has a joint overall survival function

Cheng et al. [10] defined a child–mother cause-specific association measure denoted by θCM(s,t;k,l) as follows:

which shares the same interpretation as the conditional cross hazard ratio proposed by Bandeen-Roche and Liang [4]. That is, the risk of a child experiencing a cause k event by time s would be accelerated for θCM>1 or decelerated for θCM<1 had his mother developed a cause l event by time t as compared to the case that his mother had no event by time t. This extended conditional cross hazard ratio was estimated in a piecewise manner by a U-statistic [10]. We now briefly describe the notation and methods given in Cheng et al. [10]. Similar notation will be used when we develop our modified U-statistic in Section 2.5 that specifically deals with tied events.

Let Yi,a,Yi,b,1≤a<b≤mi be the failure times of two distinct members of the ith cluster, and Yj,c,Yj,d,1≤c<d≤mj be the failure times from the jth cluster. θCM(s,t;k,l) was assumed to be a constant over each grid Ωqr, 0≤q≤n1−1,0≤r≤n2−1. Cheng et al. [10] defined Y(iajc)=min(Yia,Yjc),η(iajc)=ηiaI(Yia<Yjc)+ηjcI(Yia>Yjc) and similarly Y(ibjd) and η(ibjd) and defined the concordant indicator

and the discordant indicator

The extended child–mother association measure was estimated by

which was the ratio of the number of concordanst pairs to that of discordant pairs, among those pairs where concordance status was determinable, and the times and causes of failure matched those in θCM.

For any 1≤j<j′≤M in the sibship data, Cheng et al. [10] defined the bivariate cause-specific densities fjj′(s,t;k,l)=lim(Δs,Δt)↓0P(s≤Tj≤s+Δs,∈j=k,t≤Tj′≤t+Δt,∈j′=l)/(ΔsΔt), for k,l=1,2 and the overall survival function Sjj′(s,t)=∫s∞∫t∞∑k=12∑l=12fjj′(u,v;k,l)dudv assuming that fjj′ is absolutely continuous. The subscripts j and j′ were dropped under the exchangeability assumption among siblings, such that f(s,t;k,l)=f(t,s;l,k), for any pair of causes k,l at any time point (s,t). The sibship cause-specific association measure assuming exchangeability was defined as

The estimation of the sibship association θCC is more complicated than that of the child–mother association θCM, as under pairwise exchangeability it is meaningful to define concordance and discordance indicators based on two different pairings (Yia,Yjc) with (Yib,Yjd) as well as (Yia,Yjd) with (Yib,Yjc) [17]. Hence θCC is estimated by

The estimators (4) and (6) are U-statistics whose asymptotic properties can be obtained with U-statistic theory [10].

2.4 Maximum-pseudo (composite)-likelihood estimator

One can show that ζCC is equivalent to θCC defined in eq. (5) and that ζCM is equivalent to θCM defined in eq. (3), following the same argument as shown in Cheng and Fine [7]. Even though the data are multivariate, the association measures are still defined for paired event times, and the argument in Cheng and Fine [7] can be applied directly. Hence the two sets of estimators, the U-statistics and the plug-in estimators, are estimating the same quantities. We now discuss another estimating procedure for θ s adapted from the Clayton’s [13] pseudo-likelihood function for bivariate survival data. The likelihood was also called a composite likelihood as discussed in Lindsay [18]. Ning and Bandeen-Roche [19] recently adapted the pseudo-likelihood to bivariate competing risks data under a regression setting. We now further extend the likelihood to multivariate competing risks data using a stratified method.

We first divide mother–multiple children data into several strata, each containing pairs of mothers and their dth children. Here d ranges from 1 to mmax, where mmax≡maxi=1nmi is the maximum number of children among all families. Without loss of generality, we consider the first stratum of mothers and the corresponding eldest children. Define the cause-specific set A(s,t;k,l) of size a(s,t;k,l) pairs satisfying (1) both members are not censored by time (s,t), and a(s,t;k,l)≥1; (2) there is at least one type l event that occurred at time t from one of the mothers and no event just prior to time s from her eldest child; and (3) there is at least one type k event that happened at time s from one of the eldest children and no event just prior to time t from the corresponding mother. Let Δ(s,t;k,l)=1 if the types k and l events were from the same pair and 0 otherwise. It can be shown that P{Δ(s,t;k,l)=1|A(s,t;k,l)}=θCS(s,t;k,l)a(s,t;k,l)−1+θCS(s,t;k,l), and P{Δ(s,t;k,l)=0||A(s,t;k,l)}=a(s,t;k,l)−1a(s,t;k,l)−1+θCS(s,t;k,l).

We now construct the pseudo-likelihood function based on the observed mother–multiple children data {(Yi0,…,Yi,mi,ηi0,…,ηi,mi,mi),i=1,…,n} from the Cache County Study. To incorporate the multiple children structure, we adopt a stratified method. Let Dd denote the set containing all the families that have a dth child, where d=1,…,mmax. Then mother–child pairs belonging to Dd form a bivariate competing risks data set {(Yid,ηid,Yi0,ηi0),i∈Dd} for which the pseudo-likelihood can be constructed. As discussed before, we assume a constant θCM(Ωqr) over the region Ωqr. To find all (s,t) pairs in Ωqr such that a(s,t;k,l)≥1 and Δ(s,t;k,l)=1, we count each child–mother pair in Dd as a concordant pair, if two event times from a pair are contained in Ωqr and their causes are (k,l), denoted by the indicator function ϕi,d0qr=I{(Yid,Yi0)∈Ωqr,(ηid,ηi0)=(k,l),}, i∈Dd. The contribution of such pair to the log-likelihood function is ϕi,d0qrlog(θCMa(Yid,Yi0,k,l)−1+θCM). Next, to find (s,t) such that Δ(s,t;k,l)=0, i.e. the events k and l are not from the same pair, we check any two different pairs (Yid,ηid,Yi0,ηi0) and (Yjd,ηjd,Yj0,ηj0) in Dd. If the times of the two pairs are discordant, s will be the minimum of children’s failure times Y(idjd) and t will be the minimum of mothers’ failure times Y(i0j0). The pair of the minimum failure times (Y(idjd),Y(i0j0)) needs to be in Ωqr and the cause indicators of the minimums (η(idjd),η(i0j0)) need to match (k,l). We call two such pairs “discordant pairs” and define the discordant indicator ψij,dd00qr=I{(Yid−Yjd)(Yi0−Yj0)<0, (Y(idjd),Y(i0j0))∈Ωqr,(η(idjd),η(i0j0))=(k,l)}, where i<j and i,j∈Dd. Their contribution to the log-likelihood function is ψij,dd00qrlog(a(Y(idjd),Y(i0j0),k,l)−1a(Y(idjd),Y(i0j0),k,l)−1+θCM). Combining all the bivariate competing risks data from different strata Dd together, d=1,…,mmax, we have the following pseudo-log-likelihood function:

The maximum-pseudo-likelihood estimator θ˜CML can be obtained by using some standard minimization procedure, e.g. function “nlminb” in R, for the negative log likelihood. As pointed out by Oakes [14], the variance of θ˜CML cannot be obtained from the second derivative of the pseudo-log-likelihood function as ψij,dd00qr may be dependent when θCM(Ωqr)≠1. Instead, we use some results for Z-estimators in empirical processes [15] to derive the asymptotic property of θ˜CML based on the estimation equation Ψn(θCM)=θCM (2n)L˙n(θCM)=0, where the superscript dot denotes the derivative. The inference will be based on the nonparametric bootstrap method.

The above framework can be further extended to clustered sibship data for which exchangeability can be assumed among all siblings. The adaptation of the maximum-pseudo-likelihood estimation to the sibship data is more complicated as there are two meaningful ways to form a pair, (Yia,Yib) and (Yib,Yia), for any two children a,b from the ith family, under the exchangeability assumption. In addition, there may be too many different combinations of a and b from multiple children families. To simplify the problem, we continue to use our stratified method and denote the set Dab, 1≤a<b≤mmax, which consists of all families having ath and bth children. To adapt the exchangeability assumption, we consider the bivariate competing risks set consisting of all child–child pairs {(Yia,ηia,Yib,ηib),i∈Dab} together with its switched set containing all pairs {(Yib,ηib,Yia,ηia),i∈Dab}. The cause-specific set A(s,t,k,l) is naturally extended to this bivariate set and its switched set, with the size of a(s,t,k,l). As before let Δ(s,t;k,l)=1 if cause k event at time s and cause l event at time t were from the same pair and Δ(s,t;k,l)=0 otherwise. We can still show that P{Δ(s,t;k,l)=1|A(s,t;k,l)}=1−P{Δ(s,t;k,l)=0|A}= θCC(s,t;k,l)a(s,t;k,l)−1+θCC(s,t;k,l).

We now construct the pseudo-log-likelihood function based on the stratum Dab. As before we assume constant sibship association θCC within time region Ωqr. To find all (s,t) pairs in A(s,t;k,l)∩Ωqr such that Δ(s,t;k,l)=1, we count each concordant sibship pair in the bivariate set and its reversed set and define concordant indicators ϕi,abqr=I{(Yia,Yib)∈Ωqr,(ηia,ηib)=(k,l)} and ϕi,baqr=I{(Yib,Yia)∈Ωqr,(ηib,ηia)=(k,l)}, i∈Dab. It is more complicated to find all time points (s,t)∈Ωqr such that Δ(s,t;k,l)=0, because for any two pairs, i<j, from Dab, we may have four different pairings: (Yia,Yib) with (Yja,Yjb), (Yia,Yib) with (Yjb,Yja), (Yib,Yia) with (Yja,Yjb), and (Yib,Yia) with (Yjb,Yja). For the first pairing we define the discordant indicator ψij,aabbqr=I{(Yia−Yja)(Yib−Yjb)<0, (Y(iaja),Y(ibjb))∈Ωqr,(η(iaja),η(ibjb))=(k,l)}. Similarly, we define discordant indicators ψij,abbaqr, ψij,baabqr, and ψij,bbaaqr. Combining all the strata Dab together we have the following pseudo-log-likelihood function:

The maximum-likelihood estimator θ˜CCL for the constant association measure θCC can be similarly obtained by using the function “nlminb” in R. The consistency and asymptotic properties of θ˜CCL can be proved analogously as for θ˜CML in the section “Asymptotic properties of θ˜CML” in Appendix.

For each extended association measures θCC(ζCC) and θCM(ζCM), we have discussed three sets of estimators. The plug-in estimators ζˆCC and ζˆCM and the U-type estimators discussed in Cheng et al. [10] (θˆCMU,θˆCCU) are computationally straightforward to obtain. The extended maximum-pseudo-likelihood estimators proposed in this paper (θ˜CML,θ˜CCL) do not have explicit variance estimators and hence are more computationally intensive. Asymptotic normality can be established for all the three estimators. Their relative efficiencies are examined by numerical studies in Section 3.

2.5 Modified U-statistic for tied data

In practice, the data are often rounded to the nearest integers, resulting in tied events. The adaptation of the plug-in estimator to the tied data is straightforward. For the U-statistic, we followed Cheng et al. [10] and excluded those pairs where one or both event times were tied. For the pseudo-likelihood function (7), the tied observations do not affect concordant pairs much, since a mother–child pair with double events is still counted as a concordant pair after rounding. However, rounding tends to underestimate the number of discordant pairs. Suppose we have two discordant pairs (Yi0,ηi0,Yi1,ηi1) and (Yj0,ηj0,Yj1,ηj1) from two distinct families i and j, where Yi0≤Yj0,ηi0=k and Yi1≥Yj1,ηj1=l. When ties exist, e.g. Yi1=Yj1, and if ηi1=l, then the pairing of these two pairs would not be considered as discordant any longer.

We will show through simulation studies in Section 3.4 that the plug-in estimator tends to underestimate while the other two tend to overestimate the true positive association with the presence of tied data. For the pseudo-likelihood estimator, the ties decrease the number of discordant pairs, and hence increase the estimated strength of association. For the U-statistic, rounding causes the removal of roughly similar amounts of concordant and discordant pairs. When the association is positive, i.e. the numerator is larger than the denominator, subtracting the same positive value from both the numerator and the denominator makes the ratio larger. Therefore, the U-statistic overestimates the positive association when ties are present. The similar arguments can be used to show that the U-statistic results in an estimate that is smaller than the true value when the association is negative (0<θ(ζ)<1) and ties are present. On the other hand, the pseudo-likelihood estimator still overestimates θ, since the number of discordant pairs is getting smaller. However, it is not clear theoretically how rounding errors would affect the plug-in estimator. Our simulation study suggests that the plug-in estimator tends to bias toward independence for both positive and negative associations.

To address potential loss of accuracy in the three estimators for tied data, we propose a modified U-statistic which essentially adds back those removed concordant and discordant pairs due to tied events. When two mother–child pairs have tied observations – for example, when both mothers had the same event time – we consider their corresponding cause indicators. If mother 1 died and mother 2 had dementia at the same age, we assume that mother 1 would have dementia later than mother 2, had mother 1 not died. If child event times were not tied, we could determine concordant status. For those tied pairs that can be equally concordant or discordant, we add half a pair to both concordant and discordant pairs. It becomes more complicated for sibship data under the exchangeability assumption. The section “Formula for corrected U” in Appendix provides details on how to add back all missing concordant and discordant pairs for sibship data.

The idea may be used analogously for the pseudo-likelihood estimator. In the section “Asymptotic properties of θ˜CML” in Appendix, we show that the pseudo-likelihood estimator is a Z-estimator based on the estimating equation which is a U-statistic of order 2 with the kernel function h(Xi,Xj;θ) given in eq. (11). One may let both ϕij,dd00qr and ψij,dd00qr in the numerator of eq. (11) be 0.5 for tied events following the same scheme as in the section “Formula for corrected U” in Appendix. However, it is not clear how to fix the plug-in estimator directly when there are rounding errors. Some random perturbation method may be used by adding simulated small error terms to tied observations. However, multiple imputation is needed which is beyond the scope of this paper.

2.6 Grid selection

The selection of grids plays an important role in capturing the local association structure. If we use too few grid points, the resulting estimates may be too smooth, missing crucial local associations. On the other hand, if there are too many grid points, the number of events would be too few and there would be too much variability in the estimators. In practice, if there is no particular scientific definition of time regions in place, we recommend that we start with some crude divisions and then further refine the grid points. We will carry out formal hypothesis tests on two neighboring associations being the same, based on which we decide whether we should keep the finer grid points. That is, we would like to test the following hypotheses:

where Ωqr≡(τ1,q,τ1,q+1)×(τ2,r,τ2,r+1) with 0≤q≤n1−1,0≤r≤n2−1. For any pair of two adjacent regions, e.g. Ωqr and Ωq+1,r, we construct a Wald type of statistic based on the difference in the two local associations. That is,

where var{ζˆCM(Ωqr)} and var{ζˆCM(Ωq+1,r)} are the variance estimators of ζˆCM(Ωqr) and ζˆCM(Ωq+1,r) computed based on the influence functions given in eq. (10). When the number of events in each grid is not too few, we can compare this test statistic with its asymptotic chi-squared distribution with one degree of freedom. For the two neighboring regions Ωqr and Ωq,r+1, similar tests can be constructed. Due to the equivalence of ζ and θ, the above hypotheses can be equally expressed in terms of θ and be tested based on the U-estimator or the pseudo-likelihood estimator, though for the latter estimator the bootstrap method is used to obtain the variance. The same idea can be applied analogously to the sibship association analysis.

3.1 Constant association

We conducted simulation studies to evaluate small-sample performance of the plug-in estimator ζˆ, the maximum-pseudo-likelihood estimator θ˜L, and the U-statistic θˆU. We first simulated mother–multiple children data by extending Bandeen-Roche and Liang’s [4] frailty model. They modeled cause-specific between-subject dependency by an overall frailty A and proportions of cause allocations B which was between 0 and 1 and independent of A. As the Archimedean copula is associative, we naturally extended their frailty model for the bivariate competing risks to multivariate settings in our simulations. Three different strengths of association were considered. For each scenario, we generated 500 replicates of 500 families with varying family sizes. More details of our simulation are given below.

For each family, we first drew the familywise propensity A from a Gamma (1,1). The number of children for a family M was determined from a multinomial distribution ranging from 1 to 5 with probabilities 0.09, 0.18, 0.24, 0.24, and 0.25 which were estimated from the Cache County Data. Next we generated M+1 failure times from Exp(A) for the mother and her M children. Since the Cache County Study data were collected retrospectively, some previous investigation [20] showed that the marginal cumulative incidence function (CIF) for a mother is similar to that for her eldest child. For the sake of simplicity, we used the same distribution for mother and children in this simulation. Cause indicators were generated from a Bernoulli distribution with probability B for each subject, where B was drawn from a Beta(R1,1−R1). Three values of R1, 0.2,0.5, and 0.8, were considered with corresponding theoretical values of θCM(s,t;1,1) being 6.0,3.0, and 2.25 for all s and t. Finally, independent censoring times were generated from Uniform(5,10) for each subject in this family, imposing about 12% censoring.

We first adopted the same frailty model as those used in Bandeen-Roche and Liang [4] and Cheng and Fine [7]. The model assumes the same strength of association between causes and results in constant association over time. Hence the estimation is carried out over the entire observation region with uniform weight. We applied the three methods to each of the 500 simulated data sets. The mean of estimates, mean of standard error (for the plug-in and U-estimators), or bootstrap standard error (for the pseudo-likelihood estimator) using family-based resampling with replacement with 500 replications, empirical standard error, and coverage rate of 95% confidence intervals are summarized in the upper left panel of Table 1. With moderate sample size, light censoring, and 50–80% of cause 1 events, all the three estimators perform well. The estimates are very close to the true values. The model-based or bootstrap standard errors agree well with the empirical standard errors. The coverage rates are close to the nominal level 0.95. There is some bias in the pseudo-likelihood estimate and U-statistic when R1=0.2, since there were fewer cause 1 events. Similarly, we simulated 500 sibship data sets using the same simulation strategy. Hence the child–child association θCC takes on the same values as the mother–child association θCM. Each data set contained 250 families. The results are summarized in upper right panel of Table 1. Since there are 250 families instead of 500 in each sibship data, there is some discrepancy between the estimates and the true value of association when the number of cause 1 events is small (R1=0.2).