Marginal Methods for Multivariate Failure Times Under Event-Dependent Censoring

Many chronic diseases put individuals at increased risk of several different types of adverse clinical events. Typically these events are combined to define composite events which are then used as the basis of treatment evaluation. A potentially more efficient approach is to conduct separate marginal assessments of the effect of treatment on each component and then to synthesize this information across each type of event. While there is considerable potential for more powerful tests of treatment effect in this setting, it is possible that dependent censoring can cause problems. This happens when the occurrence of one type of event increases the risk of withdrawal from a study and hence alters the probability of observing events of other types. The purpose of this article is to formulate a model which reflects this type of mechanism, to evaluate the effect on the asymptotic and finite sample properties of marginal estimates, and to examine the performance of estimators obtained using flexible inverse probability weighted marginal estimating equations. Data from a motivating study are used for illustration.


Introduction
Many chronic disease processes put individuals at risk for several types of events and it is often of interest to examine the effect of treatment on the risk of occurrence for each type (Dabrowska, 2006;Hougaard, 2000).In settings involving life history processes, multiple events can occur during a particular period of observation and composite endpoints are also routinely used as a basis for treatment assessment (Freemantle, Calvert, Wood, Eastaugh, & Griffin, 2003).In time to event data, a composite endpoint uses the time of the first event as the response, regardless of the type, and is appealing since it permits the use of standard methods for survival analysis (Kalbfleisch & Prentice, 2002;Lawless, 2003).It does not, however, lead to treatment comparisons based on a full characterization of the disease process.For this reason, in clinical trials investigators have increasingly turned to use of multiple endpoints and regulatory agencies are increasingly requiring demonstration of efficacy new interventions based on such analyses (Buzney & Kimball, 2008;Fleming & Lin, 2000;Freemantle & Calvert, 2007;Wei & Glidden, 1997).
There are three common frameworks for the analysis of multivariate failure time data including frailty-based models (Therneau & Grambsch, 2000), copula models (Liang, Self, Bandeen-Roche, & Zeger, 1995;Nelsen, 2006), and marginal methods (Wei, Lin, & Weissfeld, 1989).While frailty models and copula models yield multivariate distributions, they require assumptions regarding the frailty distribution or the copula function respectively.Such fully specified models are appealing if interest lies in estimating the degree of association between two or more event times or for prediction, but when assessing treatment effects in clinical trials it is generally desirable to make minimal assumptions and maintain robustness.The marginal approach of Wei et al. (1989) requires only specification of a Cox regression model for each type of event but no dependence structure is specified for distinct failure times do not have a specific dependence structure.Simultaneous inference regarding the estimates of the marginal regression coefficients is carried out through use of a robust sandwich type variance estimator.This method is easily implemented in most major statistical software packages such as R/S-PLUS and SAS (Therneau & Grambsch, 2000) and is widely used in clinical trials (Lin, 1994).
The marginal approach of Wei et al. (1989) is based on a working independence assumption and the robust covariance matrix and hence has similarities with the generalized estimating equations approach of Liang and Zeger (1986) for dealing with clustered categorical data.A number of methodological advances have been made in the field of multivariate failure time data analysis which are based on a similar framework (Cai, Fan, Jiang, & Zhou, 2007;Cai, Fan, H. Zhou, & Y. Zhou, 2007;Cai & Prentice, 1995, 1997;Cai & Schaubel, 2004;Clegg, Cai, & Sen, 1999;Clegg, Cai, Sen, & Kupper, 2000;Kang & Cai, 2009;Lee, Wei, & Amato, 1992;Liang, Self, & Chang, 1993;Spiekerman & Lin, 1998;Greene & Cai, 2004;Yin & Cai, 2004, 2005).Since the marginal approach of Wei et al. (1989) is based on a partially specified model, however, it is only valid if censoring is completely independent of the failure time process.In studies of life history processes, when individuals are to be followed after the occurrence of events, it is common for censoring to be associated with occurrence of one or more particular types of events, yielding event-dependent censoring.For example, if the occurrence of the first event alerts a physician to the fact that the current treatment "failed" for a patient, it may increase the risk that they will be withdrawn from the study.In general, when marginal regression models are applied to multivariate failure times under such a dependent censoring scheme, biased (martingale) estimating equations are specified and the resulting estimators are inconsistent (Cook, Lawless, & Lee, 2007;Cook, Lawless, Lakhal-Chaieb, & Lee, 2009).The purpose of this article is to study the asymptotic bias that may result from a dependent censoring scheme of this sort under a working independence assumption.We also explore the utility of inverse probability of censoring weighted estimating equations developed by Rotnitzky and Robins (1995) and Robins and Finkelstein (2000).for marginal analysis of multivariate failure time data as a means of correcting for event-dependent censoring.
The remainder of this article is organized as follows.In Section 2 we provide the notation and formulate the models for the multivariate failure times and censoring processes.In Section 3 we investigate the asymptotic and empirical biases arising from a naive marginal approach.The inverse probability of censoring weighted marginal estimating equations are developed in Section 4, where we investigate the large sample properties of resulting estimators of treatment effect and the baseline hazard functions.Sections 5 and 6 contain simulation studies and an illustrative example respectively, and concluding remarks are made in Section 7.

Model Formulation for Multivariate Failure Times
Let T k denote the time of the type k event and {N k (s), 0 < s} denote the corresponding right-continuous counting process, where N k (t) = I(T k ≤ t) indicates that the type k event has occured at or before time t, dN k (t) = 1 if a type k event occurs at time t, and dN k (t) = 0 otherwise.We further let N(t) = (N 1 (t), . . ., N K (t)) and remark that the multivariate counting process {N(s), 0 < s} is specified when interest lies in jointly modeling the occurrence of all event types.Suppose Z(t) is a vector of fixed, exogeneous or endogenous covariates and let {Z(s), 0 < s} denote the covariate process.The full history at t contains information on the number and times of events over [0, t) and covariate data over [0, t] and is denoted The intensity function for type k events is Δt , where ) is the number of the events occurring over the interval [t, t + Δt).The inclusion of a dependence on the history for process in the intensity for type k events accommodates the association between processes.For continuous time processes where at most one event can occur at any time, these intensity functions fully define the multivariate counting processes (Gill, Keiding, & Andersen, 1993).
While this formulation completely specifies the multivariate process, in the context of clinical trials it is undesirable to assess treatment effects conditional on endogeneous variables (Kalbfleisch & Prentice, 2002) and hence intensity functions do not offer an appealing framework for analyses.Instead treatment effects are more naturally expressed in terms of marginal proportional hazards regression models of the form where λ 0k (t) is an unspecified positive function, β k is a regression parameter and Z is a fixed covariate which equals 1 for individuals receiving the experimental treatment and zero for those receiving a control therapy.The marginal hazard ratio reflecting the effect of treatment on type k events is then simply exp(β k ).The cumulative baseline hazard function is Λ 0k (t) = t 0 λ 0k (u)du and the marginal survivor function is ), where θ k contains the regression coefficient β k and the parameters indexing λ 0k (•).When marginal models of this type are specified it is necessary to address the association in the failure times differently than is done for intensity-based analyses.This is conveniently achieved using copula functions (Nelsen, 2006).
A copula function C φ (u 1 , . . ., u K ) in K dimensions defines a multivariate distribution on the unit hypercube [0, 1] K with uniform margins.Parametric copula functions are indexed by a parameter denoted by φ, which characterizes the association between the components of the marginal quantities.Such functions offer a convenient way of constructing multivariate distributions with marginal distributions of a specified form.Specifically, the marginal probability integral transformation of each random variable can be applied to create a K dimensional vector of uniform random variables.These in turn are then viewed as the components of a multivariate uniform random variable with their joint distribution governed by a given copula.Thus the joint survival function F 12 (t 1 , t 2 |z) = P(T 1 ≥ t 1 , T 2 ≥ t 2 |z) can be specified by linking the two marginal survival functions via a copula function where Ω = (θ , φ) with θ = (θ 1 , θ 2 ) .The Clayton copula is widely used in survival analysis and yields a joint survival distribution of the form The degree of association between two failure times is often expressed in terms of Kendall's τ which is given by τ = φ/(φ + 2) (0 ≤ τ ≤ 1) for the Clayton copula where τ = 0 and τ = 1 correspond to the cases of independence and perfect association respectively.

A Model for Event-Dependent Censoring
When multiple clinical events arise investigators often withdraw patients from trials if there is a perception that the randomized treatment is no longer appropriate.If subjects are censored at the time of study withdrawal, this creates a type of event-dependent censoring which leads to inconsistent parameter estimates under partially specified models.Consider a setting in which the intention is to follow individuals over the interval [0, C † ) where C † is a time of administrative censoring.Let C denote a random time of withdrawal where 0 and {N C (s), 0 < s} be the counting process for the random censoring time where dN C (t) = 1 if random withdrawal occurs at time t and dN Z} denote the observed history for the event and censoring processes.
The intensity for the random censoring time C is which accommodates a dependence between the censoring, event times, and possibly the treatment assignment.
When the analysis of the failure times adopted for multivariate failure time data (Wei et al., 1989) is often based on a working independence assumption, it is the dependence on the event times that is particularly problematic.
The dependence on the event history can take many forms, but in what follows we consider a particular model with the censoring intensity where λ c 0 (t) is a baseline intensity for censoring and (α 1 , α 2 ) are regression coefficients which reflect how the risk of withdrawal changes upon the occurrence of type 1 and type 2 events; we write dΛ c 0 (t) = λ c 0 (t)dt.Thus exp(α k ) is the multiplicative factor by which the intensity of censoring increases upon the occurrence of a type k event, k = 1, 2, and if α 1 = α 2 = 0, min(C, C † ) is an independent right-censoring time.

Asymptotic Biases of Marginal Estimators
In this section, we investigate the asymptotic bias caused by event-dependent censoring when marginal estimating equations are specified based on a working independence assumption, as is done for the multivariate approach of Wei et al. (1989).We first present the general framework and then study the one-sample estimates in detail.The data for a sample of n independent individuals are denoted by {( Ni (s), . ., n} where we introduce the subscript i to index individuals.We assume that the marginal distribution of T ik |z i is exponential with λ k (t|z i ) = λ ik (t) = λ k exp(β k z i ) and the joint distribution of (T i1 , T i2 )|Z i is defined through a Clayton copula.The naive marginal estimating equations are where The profile estimate of the cumulative baseline hazard for type k events is then Λ0k (t; The estimate β k is obtained as the solution to and upon substitution of With completely independent censoring the above estimating equations yield consistent estimators of the cumulative baseline hazard function Λ 0k (t) = t 0 λ 0k (u)du, as well as the regression coefficient β k , k = 1, 2. If censoring is governed by an intensity featuring a dependence on the event history, (4) and ( 5) may yield inconsistent estimators (Robins & Finkelstein, 2000;Datta & Satten, 2001).The limiting value of the estimator of the cumulative baseline hazard function under a general censoring scheme is where β * k is the limiting value of β k obtained as the implicit solution to The expectation E(•) in ( 8) and ( 9) is taken with respect to the true process defined here in terms of the marginal distributions, the Clayton copula, and the event-dependent censoring intensity (3).Some details on these derivations are given in Appendix A.
To illustrate the bias of the naive marginal approach in estimation of the cumulative hazard function, we consider a separate analysis of two treatment groups and hence focus on the one-sample problem; in this case we restrict attention to (4) with z i = 1 for the treatment group and z i = 0 for the control group.The Nelsen-Aalen estimator of Λ 0k (t) is The bias of this estimator is investigated empirically by simulating a data set of n = 2000 individuals having an equal risk of events with λ k = 2, k = 1, 2 and C † chosen to give 14% administrative censoring.A Clayton copula model was used to induce an association between the failure times with τ = 0.2 and τ = 0.6.The censoring intensity was based on (3) with α 1 = log 1.3 and α 2 = log 3.5, and λ c 0 (t) = λ c 0 was chosen to give 50% net censoring for type 1 events.Figure 1 displays the true cumulative hazard function and naive estimates based on the Nelson-Aalen estimator when τ = 0.2 and τ = 0.6; the left panel contains the results for type 1 events and the right panel for type 2 events.The plots demonstrate that the naive method yields a conservative estimate of the cumulative hazard function with the magnitude of this bias increasing with time.The stronger the association between the failure times, the greater the empirical bias.It is interesting to note that the bias is greater for the estimated cumulative hazard for type 1 events since the strength of the dependence between type 2 event times and censoring is greatest.

IPCW Weighted Estimating Equations
The estimating Equations ( 4) and ( 5) can be modified to yield consistent estimators by introducing inverse probability of censoring weights (Robins, 1993).If, as in (3), the censoring intensity for individual i at time t depends on the history Hi (t), then let G i (t) = P(C i ≥ t| Hi (t)) which by the product integration (Gill et al., 1993) can be written as G i (t) = u<t [1 − dΛ c (u| Hi (u))].Furthermore let G(t) be survival function of the random right censoring time under the scenario of independent random censoring, in which case the censoring intensity is a hazard function with λ c (t)dt = dΛ c (t); note we drop the subscript 0 here since it is no longer a baseline censoring intensity.Then again by product integration we obtain The marginal estimating functions corresponding to (4) and ( 5) are then defined as respectively.The weights are introduced to ensure that the resulting marginal estimating equations are unbiased and hence that consistent estimators are obtained.We show E(U k1 (t)) = 0 by taking the expectation in stages.First we note that and then and taking the expectation of this with respect to the process history gives zero.By similar arguments it can be shown that the expectation of ( 12) is also zero.
From the derivations above it is clear that G(t) is not necessary to guarantee unbiasedness of the estimating equations.In fact it will have no role in estimation of the baseline hazard function since it cancels in the numerator and denominator of Breslow's estimator given in (11).Robins (1993) showed, however, that inclusion of G(t) in ( 12) yields estimators of β k which are more efficient than those obtained when G(t) = 1.When G(t) is set to 1 we refer to these as regular IPCW weights and when G(t) is defined more generally we refer to these as stabilized weights.
In practise, of course, to use ( 11) and ( 12) the functions G(t) and G i (t) must be consistently estimated.Let Λ c (t) = t 0 dΛ c (u)du where dΛ c (u) = λ c (u)du is the crude censoring hazard under the working independence assumption between the censoring and event processes.In this case Λ c (u) is estimated simply as which gives the usual Kaplan-Meier estimate of the survival function for the censoring distribution.
Correct specification of the model for G i (t) is more crucial since it is what renders the inverse weighted estimating functions unbiased.If one believes the censoring intensity function ( 3) is correct, one can adopt it in the following.Alternatively, we prefer to relax the proportionality assumptions in (3) and consider a more robust stratified model for the censoring intensity with dΛ c (t|H , where l, m = 0, 1.The corresponding nonparametric estimate is then if N i (t − ) = (l, m) .Then again by product integration we have Upon substituting these estimates into (11), we obtain the weighted Breslow estimator of the cumulative baseline hazard function for type k events , where β k is the estimate obtained from the weighted score function for the kth type of event: The limiting distribution of estimated regression coefficients are given by the following theorem.

in distribution to a vector of zero-mean normal random variable with a variance that can be consistently estimated by I
converges in distribution to a vector of zero-mean multivariate normal random variable with a covariance matrix that can be consistently estimated by Σ where the (l, m) element is n Remark 4.2 The variance estimator is the usual robust sandwich estimator that can be directly obtained from R/S-PLUS or SAS using a suitably constructed dataframe in the counting process format.
Remark 4.3 The global estimate of treatment is simply a linear combination of all component-specific estimates of treatment effect β k , k = 1, . . ., K, where β = ( β 1 , . . ., β K ), the weight c( β) = Σ( β) −1 J[ J Σ( β) −1 J] −1 is chosen to estimate the weight matrix to minimize the variance in the class of all linear estimators; Σ( β) is the estimate for the variance-covariance matrix of β given in Remark 4.1, and J = (1, . . ., 1) .An asymptotically equivalent combined estimate can be obtained by fitting a single Cox model by stratifying on event type and constraining the coefficients to be the same for the different types of events (Therneau & Grambsch, 2000).

Empirical Investigation
Simulation studies were conducted to assess the finite-sample performance of the estimators obtained through the IPCW marginal estimating equations.The failure times T 1 and T 2 were generated using a Clayton copula with exponential margins given the treatment assignment.Without loss of generality we set C † = 1.For a given value of the association parameter φ, λ 1 and λ 2 were determined to give a particular stochastic ordering q = P(T 1 < T 2 |z = 0) and rate of administrative censoring p for T = min(T 1 , T 2 ) (i.e.P(T < C † ) = p).
We define the intensity for the random censoring time according to (3) with λ c 0 (t) = λ c 0 .For given (α 1 , α 2 ) , λ c 0 is specified to ensure a prescribed probability of observing the first event in the control arm P(T 1 < C|z = 0) = π is satisfied.We set β 1 = β 2 = log 0.80, τ = 0.4 and varied q = P(T 1 < T 2 |z = 0) over 0.25, 0.50 and 0.75.We set α 1 = log 1.3 and α 2 = log 3.5 and set λ c 0 so that P(T 1 < C|z = 0) = 0.4.The regression coefficients were obtained by solving ( 4) and ( 5) to obtain unweighted estimates, and ( 11) and ( 12) to obtain weighted estimates with G(t) = 1 or more generally.The estimate β c of ( 15) is a pooled estimate of β 1 and β 2 .The summary statistics of the estimated regression coefficients are reported in Table 1 including the empirical bias (Bias), the empirical standard error (ESE), and the average robust standard error (ASE) based on the large sample results.The simulations were conducted with 2000 samples each of n = 500 individuals.There is generally very good agreement between the empirical and average asymptotic standard errors in all settings.While the biases are generally quite small in the unweighted analyses, it is apparent that the impact of dependent censoring is different for the two marginal parameters and the magnitude of bias is influenced by both the stochastic ordering of the events as well as the (α 1 , α 2 ) parameters.That is, the bias is greater for type 1 events since the dependence between the time of the type 2 event and censoring is greatest (this is what induces the dependent censoring from marginal analyses of type 1 events).The empirical biases of the weighted estimators are generally smaller indicating the advantages of inverse weighting.The estimates obtained using the weight 1/G i (t) have considerably larger standard errors than the estimates obtain using the stabilized weight G(t)/G i (t), which are in turn much closer to the standard errors of the unweighted analyses.Thus the weight function G(t)/G i (t) leads to estimator with the best performance in that it provides protection against bias induced by dependent censoring at the price of a relatively small increase in the standard error.

Application
Here we consider data from a trial of 380 breast cancer patients with skeletal metastases (Hortobagyi et al., 1998).
The experimental treatment is a bisphosphonate which is studied for its palliative effect of reducing the incidence of fractures and the need for radiotherapy for the treatment of bone pain.Following randomization patients were followed for up to 24 months for the occurrence of these types of events and interest lies in the marginal analysis of the time to first fracture and the time to first round of radiotherapy, as well as a global assessment of the treatment effect using a Wei-Lin-Weissfeld analysis.To address the issue of the competing risks of death we adapt the marginal analyses to be based on fracture-free survival and radiation-free survival.
Figure 2 gives plots of the cumulative intensity for censoring by fracture status (left panel) and radiotherapy status (right panel).The slope of the cumulative intensity for censoring following the occurrence of the first fracture is considerably steeper than it is for fracture-free individuals, revealing a dependence between fracture status and censoring.The same pattern is seen in the right panel in that the rate of censoring for patients who have had one round of radiotherapy is higher (reflected by the steeper slope) than those who have not required radiotherapy.These plots are suggestive of a need to deal with dependent censoring for marginal analyses.The estimates of the cumulative baseline hazard for the analysis based on fracture-free survival and radiotherapyfree survival are given in Figure 3.There is empirical evidence of a greater effect of dependent censoring in the fracture-free survival analysis in that there is a bigger difference between the unweighted and weighted estimates than is seen for the radiotherapy-free survival analysis.This is compatible with the simulation results in that the larger difference between the cumulative intensities for censoring by radiotherapy status in Figure 2 (in comparison to the estimates of the censoring intensity by fracture status), suggests the dependence between censoring and radiotherapy is greater.This in turn will have a greater effect on the estimates related to the fracture-free survival endpoint.While the effects are not large, there is a suggestion that the unweighted analysis yields a conservative estimate of the event rates since the naive estimate is lower than either of the weighted estimates.Table 2 contains the results of the marginal and global regression analyses.Here we see the estimates of the treatment effect from the use of stabilized weights are slightly larger than those obtained from an unweighted analysis.The relative risk reduction for fracture-free survival was 22.6% for the unweighted analysis compared to a 24.9% relative risk reduction from analysis using stabilized weights.Moreover, in contrast to the unweighted analysis, the results based on the stabilized weights give statistically significant evidence of a treatment benefit for fracture-free survival (p= 0.0465).Very similar estimates are seen for the radiotherapy-free survival endpoint for unweighted and weighted analyses using stabilized weights.Finally, the use of stabilized weights incurs a relatively small penalty in terms of efficiency as the standard errors are very close to those of the unweighted analyses.For the global analysis, there is a 32.3% relative risk reduction from the unweighted analysis and a 33.2% reduction from the weighted analysis using stabilized weights.Note: Weighted estimates are obtained by inverse probability weighted estimating equations with weights 1/G i (t) and stabilized weights G(t)/G i (t).

Discussion
Multivariate failure time data are frequently encountered in clinical trials and observational studies.Frailty models are a popular choice for the analysis of multivariate failure time data, but they do not yield estimates of treatment effect which have a simple marginal interpretation.Such models can be formulated using copula functions, but it is undesirable to base inferences on a particular parametric model and for this reason the marginal methods proposed by Wei et al. (1989) are preferred.Inference regarding regression coefficients in this framework are carried out by use of a robust sandwich-type variance estimate easily computed in SAS or R/S-PLUS (Therneau & Grambsch, 2000).
With multivariate failure time data, studies are usually designed to follow individuals for the occurrence of all types of events up until some administrative censoring time.In this setting however, event occurrence may cause investigators to withdraw patients from a trial if it is thought that following the protocol is no longer in the patients' best interests.This can occur, for example, when one or more clinical endpoints are observed.The simplicity of the marginal analysis of Wei et al. (1989) arises from the working independence assumption.This enables the use of standard software for point estimation, but the validity of this hinges on the censoring being independent of the event processes.When this is not satisfied inconsistent estimates are obtained for all marginal parameters including the cumulative baseline hazard functions and the regression coefficients from the marginal Cox models.Use of inverse probability of censoring weights are known to address this problem (Robins, 1993) and we have shown that this strategy can be put to good use in the context of multivariate failure time data.In this study, we proposed a marginal IPCW approach to analyze multivariate failure times with event-dependent censoring and demonstrated the effectiveness of the proposed approach using simulation studies.It offers a viable approach to correcting for event-dependent censoring when data from each type of event are to be used for analyses.An alternative, less efficient approach to this problem would be to artificially censor individuals after the occurrence of the first event.This creates a competing risk problem for which one could use standard cause-specific Cox regression models.Alternatively models based on the cumulative incidence function could be employed as developed by Fine and Gray (1999).

A. The Limiting Value of Unweighted Estimators
The numerator of ( 8), E( Ȳik (t)dN ik (t)), is calculated by noting This gives In the same way, the denominator E( Ȳik (t) exp{βz i }) can be obtained as The probabilities P( Ȳik (t) = 1|z i = z) can be obtained analytically under the marginal models, Clayton copula (2) and event-dependent censoring mechanism (3).Hence the limiting value of the estimator of the baseline cumulative hazard function dΛ * 0k (t) can be obtained.The limiting value of estimators of the regression coefficients are obtained following these calculations from (9).
Proof.We first establish the asymptotic normality of √ n( β − β 0 ) using the true weight function w i (t) = G(t)/G i (t), and then prove that the effect of using a nonparametric estimation of w i (t) can be ignored in the variance estimation.Consider the weighted Cox log partial likelihood that leads to the partial score function ( 14), where R(β, u) = nS (0) (β, u).By Lemma A2 of Hjort and Pollard (1993), we can expand log R(β, u) around the true value β 0 log R(u, where V(u) = n i=1 w i (u)Y i (u) exp{Z i (u)β 0 }{Z i (u) − Zi (u)} ⊗2 /R(β 0 , u).The reminder term ν(u, x) is bounded by which can be further written as where ) and goes to zero as n → ∞.Hence, (18) can be approximated by U n t − 1 2 t I n t, which can be maximized at t = I −1 n U n .Note that its concavity in t follows from the convexity of log R(β, u) in β.Suppose β is a solution to the estimating equation ( 14) that maximizes the log partial likelihood (16), then √ n( β − β 0 ) maximizes (18).By the assumption of ( 5) and the extension of the Basic Corollary of Hjort and Pollard (1993), and the asymptomatic normality of β can be established if the asymptomatic normality of U n is established.
We now follow the arguments in Lin, Wei, Yang, and Ying (2000) to establish the asymptotic normality of U n .Let where M(t) = i≤n M i (t) and MZ (t) = i≤n t 0 Z i (u)dM i (u).For fixed time t, M(t) and MZ (t) are sum of iid zero-mean random variables.By the multivariate central limit theorem, (n − 1 2 M(t), n − 1 2 MZ (t)) converges in finite

Figure 2 .
Figure 2. Plot of cumulative intensity function for censoring (13) by fracture status (left panel) and radiotherapy status (right panel) for placebo arm patients

Table 1 .
Empirical results from simulation studies examining the frequency properties of estimators of the marginal regression coefficients (Type 1 and Type 2) and global estimates (Global) under dependent censoring

Table 2 .
Estimates obtained by fitting separate marginal Cox models and using the global Wei-Lin-Weissfeld analysis (Global) in the analysis of data from the trial of breast cancer patients with skeletal metastases; unweighted and weighted analyses