The Canadian Journal of Statistics
Vol. 42, No. 3, 2014, Pages 365–383
La revue canadienne de statistique
365
Semiparametric methods for survival analysis
of case-control data subject to dependent
censoring
Douglas E. SCHAUBEL1 *, Hui ZHANG2 , John D. KALBFLEISCH1 and Xu SHU1
1 Department
2 U.S.
of Biostatistics, University of Michigan, Ann Arbor, MI, USA
Food and Drug Administration, Silver Spring, MD, USA
Key words and phrases: Case-control study; Cox regression; dependent censoring; estimating equation;
inverse weighting
MSC 2010: Primary 62N01; secondary 62D05
Abstract: Case-control sampling can be an ef?cient and cost-saving study design, wherein subjects are
selected into the study based on the outcome of interest. It was established long ago that proportional hazards
regression can be applied to case-control data. However, each of the various estimation techniques available
assumes that failure times are independently censored. Since independent censoring is often violated in
observational studies, we propose methods for Cox regression analysis of survival data obtained through
case-control sampling, but subject to dependent censoring. The proposed methods are based on weighted
estimating equations, with separate inverse weights used to account for the case-control sampling and to
correct for dependent censoring. The proposed estimators are shown to be consistent and asymptotically
normal, and consistent estimators of the asymptotic covariance matrices are derived. Finite-sample properties
of the proposed estimators are examined through simulation studies. The methods are illustrated through an
analysis of pre-transplant mortality among end-stage liver disease patients obtained from a national organ
failure registry. The Canadian Journal of Statistics 42: 365–383; 2014 © 2014 Statistical Society of Canada
Re´ sume´ : L’´echantillonnage cas-t´emoins peut constituer un plan d’exp´erience ef?cace et e´ conomique dans
le cadre duquel les sujets sont choisis pour l’´etude en fonction du ph´enom`ene e´ tudi´e. Il est e´ tabli depuis
longtemps que le mod`ele de r´egression a` risques proportionnels peut s’appliquer a` des donn´ees cas-t´emoins.
Cependant, toutes les techniques d’estimation existantes supposent que les temps de d´efaillance sont cen´
sur´es de fac¸on ind´ependante. Etant
donn´e que l’ind´ependance de la censure est souvent bafou´ee dans le
cadre d’´etudes observationnelles, les auteurs proposent des m´ethodes pour la r´egression de Cox de donn´ees
de survie sujettes a` la censure d´ependante obtenues par un e´ chantillonnage cas-t´emoins. Les m´ethodes propos´ees se fondent sur des e´ quations d’estimation pond´er´ees dont les poids s´epar´es et inverses permettent de
tenir compte de l’´echantillonnage cas-t´emoins et de corriger le biais li´e a` la censure d´ependante. Les auteurs montrent que les estimateurs propos´es sont convergents et asymptotiquement normaux. Ils obtiennent
e´ galement des estimateurs convergents pour les matrices de covariance asymptotique. Ils examinent les propri´et´es de ces estimateurs sur des e´ chantillons de taille ?nie par voie de simulation et illustrent les m´ethodes
au moyen d’une analyse de donn´ees sur le taux de mortalit´e pr´etransplantation chez les patients atteints
d’une maladie h´epatique en phase terminale provenant d’un registre national d’organes d´efaillants. La revue
canadienne de statistique 42: 365–383; 2014 © 2014 Soci´et´e statistique du Canada
1. INTRODUCTION
Case-control and case-cohort sampling schemes are designed to over-sample subjects expected
to provide relatively large amounts of information in terms of parameter estimation. In the case* Author to whom correspondence may be addressed.
E-mail: [email protected]
© 2014 Statistical Society of Canada / Soci´et´e statistique du Canada
366
SCHAUBEL ET AL.
Vol. 42, No. 3
control study, subjects observed to experience the event of interest (the cases) are over-sampled,
while only a fraction of subjects not experiencing the event of interest (the controls) are selected. In
the original formulation of case-control study, all cases are selected but, more generally, cases are
over-sampled, such that their sampling fraction would be many times greater than that expected
through random sampling. In the context of proportional hazards regression (Cox, 1972), the casecohort design (Prentice, 1986) was proposed, in which one selects all cases, as well as a random
sample of the study population (the sub-cohort). The case-cohort design is especially well-suited
to Cox regression since the sub-cohort representing the risk set can be used to compute partial
likelihood contributions at each observed event.
A number of methods have been proposed for the regression analysis of case-control and
case-cohort studies under the proportional hazards model. Prentice & Breslow (1978) adapted
the Cox model for use in case-control studies, and studied in detail the correspondence between
Cox regression in the presence of prospective and retrospective sampling. Oakes (1981) provided a partial likelihood interpretation for the likelihood function used by Prentice & Breslow
(1978) and, hence, a theoretical justi?cation of the estimation procedure. Self & Prentice (1988)
provided a rigorous theoretical evaluation of the case-cohort design of Prentice (1986), including
detailed asymptotic derivations and various ef?ciency calculations. Kalb?eisch & Lawless (1988)
developed pseudo-likelihood methods for semiparametric estimation of the illness-death model,
applicable to case-cohort and other retrospective sampling designs. Lin & Ying (1993) considered
general missing data problems in the context of Cox regression, with the case-cohort design cast
as a special case. The asymptotic development by Lin & Ying (1993) leads to a different variance
estimator for the Prentice (1986) case-cohort estimator than that derived by Self & Prentice (1988).
Much of the subsequent research on the case-cohort design was focussed on modi?cations
to the sampling scheme and/or analysis in order to increase ef?ciency. For example, Chen & Lo
(1999) derived a class of estimators for case-cohort sampling; essentially, the methods involve
adding future cases to the risk set (in addition to the sub-cohort) in order to improve ef?ciency
relative to the original method by Prentice (1986). Borgan et al. (2000) proposed strati?ed casecohort designs to improve ef?ciency. Motivated by family-based studies, Moger, Pawitan, &
Borgan (2008) developed case-cohort methods for ?tting frailty models to large clustered data sets.
The methods were inspired by the work of Kalb?eisch & Lawless (1988) and Borgan et al. (2000).
The motivation for the work of Moger, Pawitan, & Borgan (2008) was to reduce computational
burden, rather than dollar cost.
With respect to the case-control study, much of the development in the last 15 years has been
directed at clustered data. For instance, in the case-control family study, one selects a random
sample of cases, random sample of controls, as well as the family members of the cases and
controls. Li, Yang, & Schwartz (1998) developed parametric methods for ?tting copula models to
clustered survival data obtained through case-control family sampling. Shih & Chatterjee (2002)
extended this work to the semiparametric setting through quasi-likelihood procedures. Hsu &
Gor?ne (2006) developed pseudo-likelihood methods for ?tting frailty model to case-control
family data using an approach similar to that of Shih & Chatterjee (2002). Hsu et al. (2004)
developed semiparametric methods for frailty models based on case-control family sampling,
with the focus being on the conditional regression parameter. Gor?ne, Zucker, & Hsu (2009)
further evaluated frailty models based on case-control family data through a pseudo-full likelihood
approach, and derived rigorous large-sample theory.
There have been several evaluations of cohort sampling in a general sense. Chen & Lo (1999)
described the close theoretical connection between the case-cohort and case-control study within
the framework of Cox regression. Chen (2001) considered “generalized case-cohort sampling”
(nested case-control, case-cohort and case-control, all within the same framework) in the context
of Cox regression. Gray (2009) studied various cohort sampling designs, focusing on Inverse
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
367
Probability of Selection Weighted (IPSW) estimators applied to the Cox model (e.g., Binder,
1992; Barlow, 1994; Borgan et al., 2000; Lin, 2000; Pan & Schaubel, 2008; Nan, Kalb?eisch, &
Yu, 2009).
The motivation for subsampling of large epidemiologic databases is well-established in the
literature. Nothing prevents such data from being subject to dependent censoring. However, each
of the subsampling methods cited thus far (and, to the best of our knowledge, all existing methods
for Cox regression analysis of case-control data) assumes that subjects are censored in a manner
independent of the failure rate. In most practical applications in which the independent censoring
assumption is violated, censoring is actually a mixture of independent (e.g., administrative, or
random loss to follow-up) and dependent (e.g., censoring of pre-treatment death by the initiation of
treatment). In such cases, it is necessary to distinguish between subjects who were independently
and dependently censored. The most popular method for handling dependent censoring is Inverse
Probability of Censoring Weighting (IPCW) proposed by Robins & Rotnitzky (1992). From this
perspective, case-control studies appear to extend nicely to the dependent censoring setting. In
particular, since IPCW entails modelling the dependent censoring hazard, there is an incentive to
over-sample dependently censored subjects.
We propose semiparametric methods for the analysis of failure time data generated by casecontrol sampling and subject to dependent censoring. In the setting of interest, we over-sample
dependently censored subjects, since correcting for dependent censoring through IPCW requires
that the dependent censoring process be modelled. The dependent censoring hazard is modelled
using existing case-control survival analysis methods featuring Inverse Probability of Sampling
Weighting (IPSW). Parameter estimation for the failure time hazard model then proceeds through
weighted estimating equations, with the weights correcting for both the case-control sampling and
the dependent censoring. We derive asymptotic properties of the proposed estimator, in a manner
which accounts for the randomness in the estimated IPCW and IPSW components. Since the
derived asymptotic variance is complicated, we suggest an alternative variance estimator that is
much easier to program and much faster to compute. Finite-sample performance of the proposed
procedures is evaluated through simulation.
In the evaluation of a proposed subsampling method, it is useful to assess not only if the method
can be implemented in a real data set, but, in addition, whether the method gives reasonable answers
in practice. Along these lines, we apply the proposed methods to an analysis of pre-transplant
survival among end-stage liver disease (ESLD) patients. In this application, the full cohort is
available, meaning that the results based on the full-cohort can serve as a target for the results
obtained through the proposed case-control sampling methods. With respect to background, liver
transplantation is the preferred method of treatment. However, there are thousands more patients
needing a liver transplant than there are available donor organs. As such, patients clinically suitable
for receiving a deceased-donor liver transplant are put on a wait-list, with the pertinent time origin
(t = 0) then being the date of wait listing. The receipt of a liver transplant does not censor death
but, naturally, does censor pre-transplant death (i.e., the time until death in the absence of liver
transplantation). In liver transplantation, patients are ranked on the wait-list based on their most
recent Model for End-stage Liver Disease (MELD) score. MELD was chosen as the scoring system
largely because it is a very strong predictor of pre-transplant mortality (e.g., Wiesner et al., 2001;
Kremers et al., 2004; Huo et al., 2005; Merion et al., 2005; Basto et al., 2008; Subramanian et al.,
2010). Since MELD is time-dependent, an analysis of the effect of baseline risk factors (i.e.,
measured at time t = 0) on the pre-transplant death hazard could result in substantial bias if liver
transplantation were treated as independent censoring. This and other related issues in the liver
transplant setting are discussed by Schaubel et al. (2009).
In certain cases, special exceptions are made under which a wait-listed patient may be assigned
a MELD score which is higher than that calculated, in an attempt to re?ect the patient’s perceived
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
368
SCHAUBEL ET AL.
Vol. 42, No. 3
medical urgency. The most frequent occurrence of such MELD exceptions is for patients with
hepatocellular carcinoma (HCC), a form of liver cancer. HCC patients are typically assigned a
MELD score of 22, which is often considerably higher than the score based on their laboratory
measures. To our knowledge, no existing analyses in the liver transplant literature have quanti?ed
whether the MELD score of 22 accurately re?ects the true “MELD-equivalent” wait-list mortality
risk faced by HCC patients. As a primary example in this article, we compare HCC versus
non-HCC patients, with the latter group categorized by baseline MELD score. Since MELD (at
time 0 but, also, after time 0) affects both death and liver transplantation probabilities, liver
transplantation is handled as dependent censoring of pre-transplant death time.
It should be noted that MELD exception scores for HCC are usually assigned at initial wait
listing (time 0). Therefore, given our analytic objective and in the interests of interpretation, it is
appropriate to compare HCC patients to non-HCC patients, with the latter categorized speci?cally
by time 0 MELD score. To use MELD at time > 0 would be tantamount to having the HCC and
comparator patients established at different times, which would render the resulting parameter
estimates uninterpretable. More generally, the perils of adjusting for time-dependent factors (such
as MELD at time >0) are well-described in the causal inference literature (e.g., Hern´an, Brumback,
& Robins, 2000). Another general setting for which adjustment for time-dependent factors is
explicitly avoided concerns the setting wherein there is a mediator. In evaluating a baseline factor
(e.g., treatment assignment), adjusting for the mediator (e.g., though a time-dependent indicator
covariate) would generally distort the resulting treatment effect; it might then be necessary to
treat the mediator as dependent censoring.
The remainder of this article is organized as follows. In Section 2, we describe the proposed
estimation procedures. In Section 3, we derive large sample properties for the proposed estimators.
We conduct simulation studies in Section 4 to investigate the ?nite-sample properties of the
proposed estimators. Section 5 provides an application of the methods to the afore-described
end-stage liver disease data. The article concludes with a discussion in Section 6.
2. METHODS
Let Z1i denote the q1 -vector of time-constant covariates for subject i (i = 1, . . . , n). Let Z2i (t)
i (t) =
be the q2 -vector of time-dependent covariates at time t, Zi (t) = {ZT1i , Z2i (t)T }T , and Z
{Zi (u) : 0 = u = t} denote the history of Zi (·) up to time t. Let Ti and Ci be the potential failure
and censoring times, respectively. We suppose that Ci = C1i ? C2i , where a ? b = min {a, b},
C1i is the censoring time due to mechanisms that are independent of Ti given Zi (0), and C2i
denotes the dependent censoring time; that is, C2i is dependent on Ti given Zi (0). Let Xi = Ti ?
Ci , Yi (t) = I (Xi = t), 1i = I (Ti = Ci ), 2i = I (C2i = C1i , C2i < Ti ), 3i = (1 – 1i )(1 –
2i ), Ni (t) = I (Xi = t, 1i = 1) and NiC (t) = I (Xi = t, 2i = 1), where I(·) is the indicator
function. The observable data are assumed to be n independently and identically distributed copies
of {Ni (·), NiC (·), Yi (·), Zi (·)}. Let ?i indicate whether or not subject i is sampled. The variate ?i
is allowed to depend on 1i , 2i and 3i so that the sampling probability can be different
for subjects who fail, subjects who are dependently censored and those who are independently
censored. Let the cohort be divided into three strata according to the outcome (1 , 2 , 3 )
such that
Lk = {i : ki = 1}, k = 1, 2, 3. Let pk = pr(?i = 1 | i ? Lk ), p = (p1 , p2 , p3 )T and
?i (p) = 3k=1 ki ?i /pk . Note that ?i (p) weights the ith subject by the inverse probability that
the subject is sampled.
We assume that the hazard of failure for individual i is speci?ed by the following proportional
hazards model (Cox, 1972),
?i {t | Zi (0)} = ?0 (t) exp{ßT0 Zi (0)},
The Canadian Journal of Statistics / La revue canadienne de statistique
(1)
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
369
where ?0 (t) is an unspeci?ed baseline hazard function for failure time, and ß0 is a (q1 + q2 )dimensional regression parameter. Note that we are chie?y interested in inferring the role of Zi (0)
on the failure time hazard, as opposed to {Zi (t) : t > 0}, for reasons of interpretation. For example,
it is straightforward to predict survival probability from time 0 using a pre-speci?ed value of Zi (0)
i (t) would be
along with parameter estimates from model (1). To do so using a model based on Z
much more complicated, and would generally require modelling the stochastic properties of the
process Zi (t).
If it were also the case that C2i was independent of Ti given Zi (0) (unlike the setting of interest),
the root of the estimating equation U(ß) = 0, where
then ß0 could be consistently estimated by ß,
U(ß) =
n
t
?i (p){Zi (0) – Z(ß, t)}dNi (t),
(2)
0
i=1
where
follow-up time, Z(ß, t) = S(1) (ß, t)/S (0) (ß, t), S(d) (ß, t) =
n t < 8 is the?dmaximum
T
exp{ß Zi (0)}, with a?0 = 1, a?1 = a, and a?2 = aaT . Estimating equai=1 ?i (p)Yi (t)Zi (0)
tions of the same general structure as (2) and arising from IPSW have been proposed by several
previous authors, for example, Kalb?eisch & Lawless (1988), Binder (1992), Borgan et al. (2000)
and Lin (2000) for the Cox model; Kulich & Lin (2000) for the additive hazards model; and Nan,
Kalb?eisch, & Yu (2009) for the accelerated failure time model.
Since Zi (t) affects both the event and censoring times, and since Zi (t) is not incorporated into
model (1), C2i would generally not be independent of Ti given Zi (0). In this case, the estimate ß
derived from (2) could be substantially biased because (2) does not accommodate the dependence
i (t), the hazard of
between C2i and Ti . We assume that conditional on the covariate history Z
dependent censoring C2i at time t does not further depend on the possibly unobserved failure time
Ti , that is,
C
?C
i {t | Zi (Ti ), Ci = t, Ti = t, Ti } = ?i {t | Zi (t), Ci = t, Ti = t}.
(3)
This fundamental assumption is called “no unmeasured confounders for censoring” (Rubin,
1977; Robins, 1993). Borrowing terminology from the competing risks literature (Kalb?eisch
and Prentice, 2002), assumption (3) allows us to identify the cause-speci?c hazard for C2i . We
assume a time-dependent Cox proportional hazards model for the right-hand side of Equation (3),
C
T
?C
i {t | Zi (t), Xi = t} = ?0 (t) exp{a0 Vi (t)},
(4)
where ?C
0 (t) is an unspeci?ed baseline hazard function for dependent censoring, Vi (t) is a
i (t), and a0 is a s-dimensional regression parameter.
s-vector consisting of functions of Z
We propose the following estimating function,
UR (ß) =
n
i=1
t
Ri (t){Zi (0) – ZR (ß, R, t)}dNi (t),
(5)
0
where
ZR (ß, R, t) =
S(1) (ß, R, t)
S (0) (ß, R, t)
S(d) (ß, R, t) = n-1
n
Ri (t)Yi (t)Zi (0)?d exp{ßT Zi (0)}
i=1
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
370
SCHAUBEL ET AL.
Vol. 42, No. 3
Ri (t) = ?i (p)Wi (t)
C
Wi (t) = ei (t) ?i (t),
t
C
T
where C
i (t) = 0 exp{a Vi (u)}d0 (u) and the function ?i (t) in the weight Wi (t) is a stabilization
factor. We consider three choices of ?i (t). One choice is ?1i (t) = 1. However, when the censoring
C
is heavy, ei (t) could be quite large and lead to instability in the estimation. In this case, the
t
†
choice of ?2i (t) = exp[- 0 exp{aT Vi (0)}dC
0 (u)] or ?3i (t) = exp[-i {t | Zi (0)}] may be more
†
appropriate, where i (t) is based on a time-to-censoring model that uses only the baseline coC
variate values, Zi (0). Hereafter, we denote Wji (t) = ei (t) ?ji (t), j = 1, 2, 3, and correspondingly
,ß
and ß
, the solutions to UR (ß) = 0 with weights W1i (t), W2i (t) and
estimate ß0 with ß
W1
W2
W3
W3i (t), respectively.
C (t)}, where
The weight W1i (t) can be estimated using exp{
i
C
i (t) =
C (t, a) =
0
t
0
C (s, a
)
exp{
aT Vi (s)}d
0
n
t
0
i=1
?
?-1
n
?
?j (p)Yj (s) exp{aT Vj (s)}? ?i (p)dNiC (s),
j=1
, the partial likelihood estimate of a0 , is the root of UC (a) = 0, where
where a
UC (a) =
n
0
i=1
t
{Vi (t) – V(a, p, t)}?i (p)dNiC (t),
(1)
(0)
is an IPSW-based estimating function, with V(a, p, t) = SC (a, p, t)/SC (a, p, t) and
T
(d)
SC (a, p, t) = n-1 ni=1 ?i (p)Yi (t)Vi (t)?d ea Vi (t) .
C (t)}, where
C {t, a
|
The weight W2i (t) can be estimated using
?2i (t) exp{
?2i (t) = exp[-
i
i
Zi (0)}]; that is, ?2i (t) is estimated using the same model (4), but only using baseline covariate
values,
C
| Zi (0)} =
i {t, a
0
t
C (s, a
).
exp{
aT Vi (0)}d
0
C (t)}, where
† (t)}, with
?3i (t) exp{
?3i (t) = exp{-
The weight W3i (t) can be estimated by
i
i
?3i (t) estimated using an additional baseline model for C2i ,
†
†
?i {t | Zi (0), Ci = t, Ti , Ti = t} = ?0 (t) exp{a† Vi (0)},
T
such that we have
† (t) =
i
† (t, a† ) =
0
t
0
T
n
i=1
†
(s, a
† ),
exp{
a† Vi (0)}d
0
0
t
?
?-1
n
T
?
?j (p)Yj (s) exp{a† Vj (0)}? ?i (p)dNiC (s),
j=1
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
371
† is the partial likelihood estimate of a† under the model for dependent censoring with
and a
†
hazard ?i (t). Weight stabilizers analogous to ?3i (t) have been suggested, for example, by Robins
& Finkelstein (2000) and Hern´an, Brumback, & Robins (2000). We propose the stabilizer ?2i (t)
as an alternative. The performance of each of W1i (t), W2i (t) and W3i (t) are compared through
simulations studies described in Section 4.
(j = 1, 2, 3) is computed as the root of the estimating equation in (5),
To summarize, ß
Wj
i (t). After estimating ß0 , the cumulative baseline hazard, 0 (t), can then
with Ri (t) replaced by R
be estimated by
W (t) =
0
n
0
i=1
t
n
-1
T Z (0)}
(s)Y (s) exp{ß
R
Wj
i (s)dNi (s).
R
=1
3. ASYMPTOTIC PROPERTIES
The following conditions are assumed throughout this section.
{Ni (·), NiC (·), Yi (·), Zi (·)}, i = 1, . . . , n are independently and identically distributed.
P {Yi (t) =1} > 0.
t
|Zij (0)| + 0 |dZij (t)| < BZ < 8, where Zij is the jth component of Zi and BZ is a constant.
There exists a neighbourhood B of ß0 such that supu?[0,t],ß?B S(d) (ß, R, u) –
s(d) (ß, R, u) -? 0 in probability for d = 0, 1, 2, where s(d) (ß, R, u) = E S(d) (ß, R, u)
is absolutely continuous, for ß ? B, uniformly in u ? (0, t], E (·) denotes expectation. Moreover, s(0) (ß, R, u) is assumed to be bounded away from zero.
(d)
(e) There exists a neighbourhood BC of a0 such that supu?[0,t],a?BC SC (a, p, u)
(a)
(b)
(c)
(d)
(d)
(d)
(d)
-sC (a, u) -? 0 in probability for d = 0, 1, 2, where sC (a, u) = E {SC (a, p, u)} is ab(0)
solutely continuous, for a ? BC , uniformly in u ? (0, t]. Moreover, sC (a, u) is assumed to
be bounded away from zero.
(f ) The matrices A(ß0 ) and AC (a0 ) are positive de?nite, where
A(ß) =
t
0
AC (a) =
0
t
s(2) (ß, R, u)/s(0) (ß, R, u) – z(ß, R, u)?2 dF (u)
(2)
(0)
sC (a, u)/sC (a, u) – v(a, u)?2 dF C (u)
(1)
(0)
with z(ß, R, u) = s(1) (ß,R, u)/s(0) (ß, R,
u), v(a, u) = sC (a, u)/sC (a, u), F (u) = E
{Ri (u)Ni (u)}, F C (u) = E ?i (p0 )NiC (u) .
(g) 0 (t) < 8, C
0 (t) < 8.
We describe the asymptotic properties of the proposed estimators in the following theorems.
– a0 converges to a mean zero
Theorem 1. Under conditions (a) – (g), as n ? 8, n1/2 a
Normal distribution with covariance AC (a0 )-1 (a0 )AC (a0 )-1 , where AC (a0 ) is as de?ned by
Condition (f) and (a) = E{?i (a, p)?2 }, with ?i (a, p) being asymptotically independent and
identically distributed for i = 1, . . . , n; we defer the de?nition of ?i (a, p) to the Supplementary
Materials document.
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
372
SCHAUBEL ET AL.
Vol. 42, No. 3
– a0 =
In the Supplementary Materials document, we show that n1/2 a
– a0 is essentially a scaled
AC (a0 )-1 n-1/2 ni=1 ?i (a0 , p0 ) + op (1); hence, n1/2 a
sum of n independent and identically distributed random quantities with mean zero and ?nite
variance. By the Multivariate Central Limit Theorem (MCLT) and empirical process theory, the
asymptotic normality is proved.
Theorem 2.
– ß , converges to a mean
Under conditions (a) – (g), as n ? 8, n1/2 ß
W1
0
zero Normal distribution with covariance A(ß0 )-1 (ß0 , R)A(ß0 )-1 , with A(ß) having been de?ned in Condition (f) and where (ß, R) = E{i (ß, R)?2 }, with i (ß, R) being independent
and identically distributed mean 0 variates (i = 1, . . . , n) asymptotically. The explicit de?nition
of i (ß, R) is provided in the Supplementary Materials.
The proof of Theorem 2 (provided in the Supplementary Materials) begins by de C (t) – C (t)} into n1/2 {
C (t; a
C (t; a
C (t; a
) –
, p
, p0 )} + n1/2 {
, p0 ) –
composing n1/2 {
0
0
0
0
0
C
C
C
C
C
1/2
1/2
0 (t; a0 , p0 )} + n {0 (t; a0 , p0 ) – 0 (t)}. Then n {0 (t) – 0 (t)} can be expressed
asymptotically as a sum of independent and identically distributed zero-mean variates, as n ? 8.
i (t) – Ri (t)} can
Combining this result and the Functional Delta Method, we can show that n1/2 {R
be written asymptotically as a sum of independent and identically distributed zero-mean variates, as n ? 8. Finally, through the Functional Delta Method, the asymptotic normality of
– ß ) is obtained.
n1/2 (ß
W1
0
is very complicated and dif?cult to
The expression for the asymptotic covariance of ß
W1
implement numerically. A practical way to estimate the variance of the proposed estimators is to
treat the weights Ri (t) as known rather than estimated. In the setting where the weight function
is known, results derived in the Supplementary Materials show that
– ß ) = A(ß )-1 n- 21
n1/2 (ß
0
0
n
‡
Ui ß0 , R + op (1)
(6)
i=1
t
‡
with Ui (ß0 , R) = 0 {Zi (0) – z(ß, R, t)}Ri (t)dMi (t) and dMi (t) = dNi (t) – Yi (t)di (t). Hence,
– ß ) is asymptotically a scaled sum of independent and identically distributed zeron1/2 (ß
0
is estimated by
mean random quantities with ?nite variance. Therefore, the variance of ß
W1
-1
R)
and
R)
-1 , where ‡ (ß, R) = E{U‡ (ß, R)?2 }, A(
ß)
ß)
ß)
‡ (ß,
‡ (ß,
A(
are calculated by
A(
i
replacing limiting values with their corresponding empirical counterparts.
– ß ) and n1/2 (ß
– ß ).
By similar arguments, the asymptotic normality holds for n1/2 (ß
W2
0
W3
0
– ß ). Therefore,
However, the covariance will be even more complicated than that of n1/2 (ß
W1
0
and ß
. Note that (6)
similarly, we can treat Ri (t) as ?xed to calculate the variance of ß
W2
W3
and ß
can be estimated by
holds when using W2 or W3 , such that the variances of each of ß
W2
W3
-1
R)
-1 , with Ri (t) being replaced by ?i (
ß)
ß)
‡ (ß,
A(
2i (t) and ?i (
3i (t), respectively.
A(
p)W
p)W
Our ?nal asymptotic result pertains to the proposed cumulative baseline hazard function
estimator.
W (t) – 0 (t)} converges to a zero-mean GausTheorem 3. Under conditions (a) – (g), n1/2 {
0
sian process as n ? 8, with an explicit covariance function estimator.
A proof of Theorem 3 is provided in the Supplementary Materials, including de?nitions
pertinent to the limiting covariance function.
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
373
Table 1: Simulation study: data con?gurations.
Parameter
1A
1B
2A
2B
3A
3B
?0
0.1
0.1
0.1
0.1
0.05
0.05
ß1
0.406
0.406
0.406
0.406
-0.406
-0.406
ß2
0.406
0.406
0.406
0.406
0.406
0.406
?C0
0.1
0.1
0.05
0.05
0.05
0.05
a1
-0.5
-0.5
-0.406
-0.406
-0.406
-0.406
a2
0.693
0.693
1.010
1.010
1.010
1.010
a3
0.095
0.095
0.095
0.095
0.095
0.095
%T
51%
51%
62%
62%
27%
27%
%C2
44%
44%
28%
28%
45%
45%
%C1
5%
5%
10%
10%
28%
28%
corr(T, C2 )
0.13
0.13
0.20
0.20
0.27
0.27
Full-cohort
2,500
2,500
2,500
2,500
2,500
2,500
p1
0.15
0.10
0.15
0.10
0.15
0.10
p2
0.15
0.10
0.15
0.10
0.15
0.10
p3
0.15
0.10
0.15
0.10
0.15
0.10
Case-control
375
250
375
250
375
250
4. SIMULATION STUDY
We investigated the ?nite-sample properties of the proposed estimators through a series of simulation studies. In each setting, the study population (i.e., full-cohort; prior to case-control sampling),
was n = 2,500. The baseline covariate, to be used in the failure time hazard model, consisted
of two independent Bernoulli variates, Zi (0) = [Z1i , Z2i (0)]T , where Z1i does not change over
time. Each of Z1i and Z2i (0) took the value 1 with probability 0.5. The independent censoring
times C1i were constant and equal to 12. After generating UT from a Uniform(0,1) distribution,
the failure time T was generated from a Cox model with hazard function
?i (t) = ?0 exp {ß1 Z1i + ß2 Z2i (0)} ,
(7)
T
by solving the equation 0 ?0 (u) exp {ß1 Z1i + ß2 Z2i (0)} du = – log UT for T , so that UT corresponds to the survival function at T . The time-dependent covariate Z2i (t) was generated as
Z2i (0)I (UT = 0.3) + {Z2i (0) + UT × int(t)} I (0.3 < UT = 0.6)
+ {Z2i (0) + UT /2 × int(t)} I (UT > 0.6) ,
(8)
where int(t) is the integer part of t. The dependent censoring time C2i was then generated from a
Cox model with hazard function
C
?C
i (t) = ?0 exp [a1 Z1i + a2 Z2i (t)Z2i (0) + a3 Z2i (t){1 – Z2i (0)}] .
(9)
We evaluated three data (full-cohort) con?gurations and, within each, two different case-control
sampling schemes. For each of the six scenarios, 1,000 replicates were generated. Parameter
speci?cations are given in Table 1. Note that the Wi (t) component of the weight was bounded by
10 in all runs.
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
374
SCHAUBEL ET AL.
Vol. 42, No. 3
Table 2: Simulation results based on 1,000 replications: setting 1A–1B.
Wi (t)
Estimator
BIAS
ESD
ASE
CP
0.152
0.928
Setting 1A: p1 = p2 = p3 = 0.15
1
W1
W2
W3
1
W1
W2
W3
ß1
ß1W1
ß1W2
W
ß1 3
ß2
ß2W1
ß2W2
W
ß2 3
0.059
0.152
0.054
0.169
0.169
0.930
0.015
0.152
0.154
0.955
0.014
0.150
0.152
0.957
-0.087
0.150
0.149
0.905
-0.040
0.171
0.168
0.930
-0.003
0.151
0.153
0.948
-0.002
0.150
0.152
0.950
Setting 1B: p1 = p2 = p3 = 0.10
1
W1
W2
W3
1
W1
W2
W3
ß1
ß1W1
ß1W2
W
ß1 3
ß2
ß2W1
ß2W2
W
ß2 3
0.061
0.190
0.187
0.933
0.060
0.213
0.206
0.925
0.021
0.190
0.189
0.951
0.020
0.188
0.187
0.947
-0.085
0.187
0.183
0.922
-0.041
0.212
0.203
0.927
-0.005
0.185
0.187
0.947
-0.004
0.184
0.186
0.948
By the speci?cations above, the time-dependent covariate Z2i (t) is correlated with the event
time Ti through Equation (8). In addition, Z2i (t) also affects the censoring time C2i via model
(9). Since only the baseline value of Z2i (t), that is, Z2i (0), is adjusted for in the model (7), the
censoring time C2i is dependent on Ti given Z1i and Z2i (0).
For each data generation, we selected failures, dependent censoring subjects and independent
censoring subjects by Bernoulli sampling, with p1 = p2 = p3 = 0.15 (Settings 1A, 2A, 3A)
which, on average, resulted in a case-control sample n = 375 subjects; or p1 = p2 = p3 = 0.10
(Settings 1B, 2B, 3B) which resulted in an average of n = 250 subjects.
For each scenario, we report bias; empirical standard deviation (ESD) across the 1,000 replicates; average asymptotic standard error (ASE); and empirical coverage probability (CP) for 95%
con?dence intervals. Four estimators were evaluated, which differ with respect to the handling
of the IPCW component. The ?rst estimator did not make an adjustment for dependent censoring
(Wi (t) = 1). The second, third and fourth estimators are the three proposed estimators, with W1
being based on the unstabilized correction for dependent censoring; while W2 and W3 refer to the
different versions of the proposed stabilized estimators based on W2i (t) or W3i (t), respectively.
Results of the simulation study are provided in Tables 2, 3 and 4. Some general conclusions are
as follows. The estimator which did not account for dependent censoring is notably biased, in all
con?gurations. In some but not all cases, the bias is partly corrected by the estimator using W1i (t)
(i.e., the unstabilized estimator). In certain cases, the bias is actually more pronounced than for
the Wi (t) = 1 estimator. The bias due to dependent censoring is largely corrected by the stabilized
versions of the proposed estimator, with the ASE being similar to the ESD and, correspondingly,
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
375
Table 3: Simulation results based on 1,000 replications: setting 2A–2B.
Wi (t)
Estimator
BIAS
ESD
ASE
CP
0.139
0.901
Setting 2A: p1 = p2 = p3 = 0.15
1
W1
W2
W3
1
W1
W2
W3
ß1
ß1W1
ß1W2
W
ß1 3
ß2
ß2W1
ß2W2
W
ß2 3
0.096
0.136
0.090
0.149
0.151
0.911
0.042
0.140
0.146
0.953
0.028
0.136
0.143
0.954
-0.104
0.133
0.134
0.876
-0.082
0.150
0.147
0.911
-0.034
0.138
0.141
0.937
-0.021
0.135
0.139
0.947
0.916
Setting 2B: p1 = p2 = p3 = 0.10
1
W1
W2
W3
1
W1
W2
W3
ß1
ß1W1
ß1W2
W
ß1 3
ß2
ß2W1
ß2W2
W
ß2 3
0.101
0.169
0.171
0.100
0.182
0.184
0.916
0.053
0.172
0.177
0.947
0.040
0.168
0.174
0.955
-0.102
0.167
0.164
0.912
-0.083
0.184
0.179
0.926
-0.036
0.171
0.171
0.944
-0.023
0.168
0.170
0.948
coverage probabilities approximating the nominal level. Of the three proposed estimators, the
stabilized version based on W3i (t) has the best performance in the simulation studies.
5. APPLICATION
We applied the proposed methods to analyse pre-transplant mortality for patients with end-stage
liver disease (ESLD). Data were obtained from the Scienti?c Registry of Transplant Recipients
(SRTR). The study population consisted of patients who were initially wait-listed for liver transplantation in the United States at age = 18 between March 1, 2002 and December 31, 2008. Patients
were followed from the date of initial wait-listing until the earliest of death, receiving liver transplantation, loss to follow-up, or last day of the observation period (December 31, 2008). The time
scale was in days. The primary outcome of interest is pre-transplant mortality. Loss to followup, living-donor transplantation and administrative censoring were considered to be independent
censoring. Dependent censoring occurred through deceased-donor liver transplantation.
The Model of End-stage Liver Disease (MELD) score is time-dependent and is updated based
on a frequency that ranges from weekly to yearly and depends on the last reported MELD. In
the current liver allocation system, patients are ordered on the wait-list primarily by descending
MELD. That is, patients with higher MELD are considered to be at greater medical urgency
and, therefore, get higher priority for transplantation. However, for hepatocellular carcinoma
(HCC) patients, the calculated MELD based on laboratory measures has generally been thought
to understate actual medical urgency. As such, a MELD score of 22 is usually assigned to an
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
376
SCHAUBEL ET AL.
Vol. 42, No. 3
Table 4: Simulation results based on 1,000 replications: setting 3A–3B.
Wi (t)
Estimator
BIAS
ESD
ASE
CP
p1 = p2 = p3 = 0.15
1
W1
W2
W3
1
W1
W2
W3
ß1
ß1W1
ß1W2
W
ß1 3
ß2
ß2W1
ß2W2
W
ß2 3
0.065
0.200
0.202
0.945
0.119
0.235
0.231
0.909
0.044
0.220
0.215
0.938
0.015
0.207
0.205
0.952
-0.260
0.216
0.206
0.757
-0.173
0.248
0.240
0.881
-0.065
0.241
0.230
0.935
-0.029
0.231
0.222
0.943
p1 = p2 = p3 = 0.10
1
W1
W2
W3
1
W1
W2
W3
ß1
ß1W1
ß1W2
W
ß1 3
ß2
ß2W1
ß2W2
W
ß2 3
0.064
0.252
0.248
0.935
0.120
0.296
0.280
0.906
0.045
0.275
0.262
0.926
0.017
0.259
0.252
0.944
-0.259
0.264
0.255
0.825
-0.170
0.297
0.292
0.905
-0.064
0.288
0.280
0.941
-0.028
0.279
0.273
0.955
HCC patient if the laboratory MELD is less than 22. The primary objective of our analysis is
to determine which range of (calculated) MELD score is actually consistent with the HCC pretransplant mortality hazard.
The factor of interest is, in a general sense, underlying liver disease at wait listing (the time
origin), classi?ed as hepatocellular carcinoma (HCC) versus not. The objective is to determine
the MELD score category to which pre-transplant HCC mortality is equivalent. Therefore, as will
be detailed later in this section, the pre-transplant death model will have HCC as the reference
category, to which all non-HCC patients (classi?ed by baseline MELD category) are compared.
The use of baseline (i.e., t = 0) MELD score is consistent with HCC diagnosis being applied
at baseline. Note that the use of baseline versus time-dependent MELD depends on the analytic
objectives. In Section 6, we describe a setting where use of time-dependent MELD is indicated.
In many studies, it has been shown that MELD is the dominant risk factor for liver pretransplant mortality. Moreover, as stated in the previous paragraph, MELD also strongly affects
the liver transplant hazard. Therefore, since the death model does not adjust for time-dependent
MELD, pre-transplant mortality and deceased-donor liver transplantation will be correlated, such
that it is necessary to account for liver transplantation as dependent censoring.
It is necessary to account for geography in our analysis. There are 60 Organ Procurement
Organizations (OPO) in the United States, each of which maintains a liver transplant wait-list
pertaining to a population size of a state, on average. Adjustment for OPO is important since
it re?ects geography, which can affect mortality in ways not re?ected by the remainder of the
adjustment covariates. With respect to transplantation, the transplant hazard differs markedly by
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
377
Table 5: Sampling design for the analysis of liver pre-transplant mortality.
Patients
OPO size
Deaths
Transplants
Censored
Total
546
4,475
Full-cohort
HCC
non-HCC
All
418
3,511
=400
418
1,499
490
2,407
401–1,300
3,676
10,816
5,230
19,722
6,076
12,794
10,477
29,347
10,588
28,620
16,743
55,951
>1,300
Total
Sampling fractions
HCC
non-HCC
All
1
1
1
=400
1
1
1
401–1,300
0.30
0.10
0.10
>1,300
0.15
0.10
0.10
Case-control Sample
HCC
non-HCC
All
418
3,511
546
4,475
=400
418
1,499
490
2,407
401–1,300
1,079
1,075
496
2,650
893
1,295
1,018
3,206
2,808
7,380
2,550
12,738
>1,300
Total
OPO owing in part to differences in organ availability. For both the time-to-death and time-totransplant models, we adjusted for OPO through strati?cation, since it may not be appropriate to
assume proportionality with respect to the 60 OPO-speci?c baseline transplant or death hazard
functions.
Speci?cs regarding the full cohort (n = 55,951), sampling fractions and case-control sample
are provided in Table 5.
Since a relatively small fraction of the full cohort was diagnosed with HCC (8%; 4,475 patients), HCC patients were over-sampled. In addition, we wanted to have suf?cient representation
from OPOs of various sizes in the case-control sample. The way we classi?ed the size of the OPO
was based on the number of patients wait-listed in that OPO in the full-cohort. Random sampling
within end-point would lead to inadequate representation of patients from smaller OPOs. Therefore, we assigned greater (lesser) sampling fractions to the small (large) OPOs. After Bernoulli
sampling, the case-control sample consisted of 12,738 patients.
There are additional issues regarding the data structure which must be taken into account. In
particular, a patient who is too sick to receive a transplant can be inactivated (usually a temporary
measure) or removed (permanent) from the wait-list. During these intervals, the patient is ineligible to receive a transplant. Therefore, an appropriate Cox model in this setting is given by the
following,
C
T
?C
ij (t) = Ai (t)?0j (t) exp{a Vi (t)},
DOI: 10.1002/cjs
(10)
The Canadian Journal of Statistics / La revue canadienne de statistique
378
SCHAUBEL ET AL.
Vol. 42, No. 3
Table 6: Analysis comparing HCC (assigned MELD of 22; reference) versus non-HCC (by lab MELD)
pre-transplant mortality.
Patients
Case-control
Case-control
Full-cohort
Unweighted: W = 1
Weighted: W2
Weighted: W2
MELD
߈
SE
p
߈
–
0
SE
–
p
߈
–
0
SE
p
HCC
“22”
0
–
–
–
non-HCC
6–8
-1.07
0.12
<0.0001 -1.02
0.15
<0.0001 -0.81
0.09
<0.0001
9–11
-0.60
0.09
<0.0001 -0.48
0.11
<0.0001 -0.50
0.08
<0.0001
12–13
-0.53
0.09
<0.0001 -0.40
0.11
0.0002 -0.24
0.08
0.002
14–15
-0.25
0.09
0.005
16–17
0.10
0.10
0.32
0.002
0.11
0.98
0.24
0.12
0.0498
0.01
0.24
0.08
0.88
0.08
0.004
18–19
0.26
0.12
0.02
0.63
0.15
<0.0001
0.52
0.09
<0.0001
20–22
0.55
0.11
<0.0001
0.83
0.12
<0.0001
0.81
0.09
<0.0001
23–24
0.90
0.17
<0.0001
1.41
0.21
<0.0001
1.12
0.11
<0.0001
25–29
1.59
0.16
<0.0001
1.93
0.19
<0.0001
1.44
0.11
<0.0001
30–39
2.38
0.16
<0.0001
2.69
0.17
<0.0001
1.86
0.12
<0.0001
40
3.68
0.29
<0.0001
3.65
0.37
<0.0001
1.79
0.20
<0.0001
where j denotes OPO number (j = 1, . . . , 60) and Ai (t) is an indicator of being active on the
wait-list (i.e., not being deactivated or previously removed) as of time t. The time-dependent
covariate vector Vi (t) includes MELD at time t (grouped into intervals: [6,8], [9,11], [12,13],
[14,15], [16,17], [18,19], [20,22], [23,24], [25,29], [30,39] and 40) with HCC patients chosen as
the reference group. These MELD categories have been suggested by several previous analyses;
the use of a categorical version of MELD negates the need to pin down its precise functional
form (which is known to not be linear). The vector Vi (t) also includes the following baseline
covariates: age group, gender, race and blood type, with age less than 40, Female, Caucasian and
blood type O as reference levels, respectively. Note that, for the intervals where the patient was
either inactivated or removed, the transplant hazard was treated as 0, as indicated by equation
(10). As such, we delete patient subintervals with Ai (t) = 0 to ?t model (10). However, since the
inactivated or removed patients are still at risk of pre-transplant death, patient subintervals with
Ai (t) = 0 are re-included in the ?tting of the pre-transplant mortality model.
The pre-transplant death model is given by
?ij (t) = ?0j (t) exp{ßT Zi (0)},
(11)
where j again represents OPO and Zi (0) contains indicators for the afore-listed MELD categories
(HCC serving as the reference), and the adjustment covariates used in the transplant model (10).
Since the IPCW weights could be very large toward the tail of the observation time, we truncated
IPCW weights with 10.
Results of the analysis are shown in Table 6. We present results for two case-control analyses:
the “unweighted” analysis (Wi (t) = 1), in which dependent censoring (liver transplantation) is
treated as independent censoring; and the proposed weighted analysis, which corrects for dependent censoring through W2i (t). We also carried out a full-cohort analysis, again based on W2i (t).
Table 6 shows that when the dependent censoring is treated as independent, HCC appears to
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
379
have pre-transplant mortality equal to that of MELD group 16–17 (see bolded entries in Table 6).
However, accounting for the dependent censoring, HCC is the pre-transplant mortality equivalent
of the MELD 14–15 group, a result consistent with the full-cohort analysis.
Based on the simulation studies in Section 4, the unstabilized weight W1i (t) was out-performed
by both W2i (t) and W3i (t) in terms of empirical bias and variance. In our application, results based
on W1i (t), W2i (t) and W3i (t) were all quite similar. It appeared from the simulations that W3i (t)
performed somewhat better than W2i (t); in each case a reduction in bias (for W3i (t) relative to
W2i (t)) was accompanied by a concomitant reduction in standard deviation. However, in the fullcohort analysis, standard errors for W2i (t) tended to be less than those for W3i (t). Hence, we chose
W2i (t) for the main line analysis. As mentioned previously, the IPCW component of the weight
is susceptible to unduly large values, which implies the use of a cap. We employed a cap of 10,
consistent with the simulation study. In the case-cohort analyses, the cap affected 16,507 records
for W1i (t) (1.55% of the total 1,066,630 records used); 2,616 records (0.25%) for W2i (t) and 3,351
(0.31%) for W3i (t). Full cohort and case-control results based on W1i (t) and W3i (t) are provided
in the Supplementary Materials.
The bias in the analysis which does not account for dependent censoring (Wi (t) = 1) is substantial, particularly upon consideration of the objectives of our analysis. This is clear upon comparing
the estimators based on Wi (t) = 1 and W2i (t), either through relative bias, or simply the absolute
difference in point estimates relative to the standard error. By either of these measures, the bias
is considerable, except for the two most extreme categories: MELD 6–8 and MELD=40. Note
that these two extreme categories are of least interest, since it was widely anticipated that the
lowest (highest) MELD category would have much lower (higher) pre-transplant mortality risk
than HCC. For the unweighted analysis, the HCC-mortality-equivalent is given by MELD 16–17,
whereas for the weighted analysis, it appears that MELD 14–15 and HCC have equal mortality.
From this perspective, the difference between the unweighted and weighted analysis would be
considered very important to the liver transplant ?eld.
6. DISCUSSION
In this article, we propose methods for proportional hazards modelling of survival data subject to
dependent censoring and obtained through case-control sampling. The proposed methods employ a
double-inverse-weighting scheme, through which the proposed estimators adjust for the sampling
bias and overcome dependent censoring. Simulation studies show that the proposed estimators
are approximately unbiased and that our asymptotic results are applicable to ?nite samples. The
proposed estimates can be computed using standard software that accommodates left-truncation,
weighting and offsets (e.g., PROC PHREG in SAS; coxph in R).
The case-control sampling scheme we consider involves either sampling all observed failures (cases), or a randomly selected fraction. The censored subjects (controls) are categorized
as either censored dependently, or independently. There is an incentive to over-sample dependently censored subjects, such that one can estimate the IPCW weights with suf?cient precision.
Although the sampling scheme we study involves random sampling based only on ?nal status
(dead, censored), the methods extend in a straightforward fashion to the setting in which sampling
also depends on ?xed covariates, a property illustrated through the application to the liver data
on Section 5. Although we consider case-control sampling explicitly, the proposed methods also
apply if instead case-cohort sampling is utilized.
One could also view the sampling scheme considered in this article to be a form of outcomedependent sampling (ODS). In an ODS design, one collects covariate information from a sample
by allowing selection probabilities to depend on individuals’ outcomes (e.g., death, survival). An
ODS design concentrates resources on observations carrying the greatest amount of information.
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
380
SCHAUBEL ET AL.
Vol. 42, No. 3
There is a large literature on analysing data arising from ODS; see, for example, Breslow &
Holubkov (1997a), Zhou et al. (2002), Schildcrout & Heagerty (2008), Song, Zhou, & Kosorok
(2009) and Wang et al. (2009). It would be possible to generalize the proposed methods to the ODS
setting, interpreting the version we develop in this article as a censored-data-speci?c instance of
inverse-weighted generalized estimating equations.
We propose three different weights to correct the bias induced by dependent censoring. In
general, when the dependent censoring is light or moderate, the unstabilized weight W1 (t) works
well. However, when censoring is heavy, W1 (t) may be quite large toward the tail of the observation
time resulting in unstable estimates. In this case, stabilized weights, W2 (t) and W3 (t), may be
preferable and usually result in a less biased and more precise estimator than that arising from
using the unstabilized weight W1 (t). With respect to the stabilized weights, we found that W3 (t)
tended to out-perform W2 (t), although not uniformly. An apparent advantage of W2 (t) is that only
one dependent censoring hazard model is required. However, details regarding the programming
and ability to exploit functions currently available (in SAS or R) lead to W3 (t) actually being the
easier of the two stabilized weights to compute, even though a second hazard regression model is
required. Overall, it would appear that W3 (t) would generally be the preferred weight, although the
stabilized weight that actually performs best is likely to be application-dependent. For instance,
in our analysis of the liver disease data, W2 (t) resulted in parameter estimators with much smaller
standard errors than W3 (t), leading to our preference for the former in this particular case.
In simulation studies, we treated the IPCW weights and IPSW weights as ?xed to simplify
the computation, which would result in conservative covariance estimators because those weights
are actually estimated as opposed to being known. However, simulation results suggest that the
proposed ASEs by treating the IPCW weights and IPSW weights as ?xed are quite accurate. The
simpli?cation of the variance calculation would likely lead to invalid inference when the variance
of the estimated weight was substantial, which would be in cases where the subsample was not
large enough or did not contain a suf?cient number of failures or dependently censored subjects.
In such cases, the variability due to estimating the weight parameters would not be a negligible
fraction of the total variability, meaning that the suggested approximation would be inaccurate
and coding of the full variance (or, perhaps application of an appropriate bootstrapping method)
would be required.
The proposed methods require the consistency of the IPCW weight. Therefore, the proportional
hazards model for dependent censoring should be correctly speci?ed. This may be approximately
true when a suf?cient number of both baseline and time-dependent covariates are incorporated.
In our analysis of chronic liver disease patients wait-listed for liver transplantation, we found
that (calculated) Model for End-stage Liver Disease scores of 14–15 were consistent with the
pre-transplant mortality hazard of patients granted an exception due to hepatocellular cancer. The
results based on case-control sampling were consistent with those from the full-cohort analysis.
In contrast, if no adjustment was made for dependent censoring, it appeared that MELD 16–17
and HCC had equivalent pre-transplant mortality. Therefore, our results indicate that the current
MELD exception score of 22 granted to HCC patients overstates the actual medical urgency,
and that an assigned MELD score of 14 or 15 may be more appropriate. If this adjustment were
actually implemented, it would considerably decrease the frequency with which HCC patients
receive liver transplants. The results of our analysis provide currently the most de?nitive evidence
that HCC patients are given inappropriately high priority for liver transplantation in the United
States.
In certain settings, the factor of interest (in addition to the adjustment covariates) may be timedependent. In our analysis of the liver disease data, HCC designations are typically made at time 0,
implying that the appropriate comparison is to patients not classi?ed as HCC at time 0. However,
if we were interested in evaluating exception scores generally, then the most useful analysis which
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
381
addressed the appropriateness of (what would primarily be non-HCC) exception scores would
probably use time-dependent factors in the death model. This is because a hepatologist may
observe a patient clearly getting worse, but the patient’s MELD score not increasing, implying
that the score does not accurately re?ect the patient’s need for liver transplantation. This analysis
is currently under our consideration, and it appears that the proposed methods cannot readily be
extended to accommodate such settings.
APPENDIX
Proofs of Theorem 2 and Theorem 3 are provided in the Supplementary Materials document.
ACKNOWLEDGEMENTS
The authors thank the Associate Editor and Referees for their constructive comments and suggestions, which resulted in a substantially improved manuscript. This work was supported in part by
National Institutes of Health Grant R01 DK070869. The authors also thank the Scienti?c Registry
of Transplant Recipients (SRTR) access to the end-stage liver disease data. The SRTR is funded
by a contract from the Health Resource and Service Administration, U.S. Department of Health
and Human Services.
BIBLIOGRAPHY
Barlow, W. E. (1994). Robust variance estimation for the case-cohort design. Biometrics, 50, 1064–1072.
Basto, S. T., Perez, R., Villela-Nogueira, C., Fernandes, E., Barroso, A., Victor, L., Pereira, B., Ribeiro, J.
& Coelho, H. (2008). MELD as a predictor of long term mortality in liver transplantation waiting list.
Liver Transplantation, 14, S219–S219.
Binder, D. A. (1992). Fitting Cox’s proportional hazards models from survey data. Biometrika, 79, 139–147.
Borgan, Ø., Langholz, B., Samuelsen, S. O., Goldstein, L. & Pogoda, J. (2000). Exposure strati?ed casecohort designs. Lifetime Data Analysis, 6, 39–58.
Breslow, N. E. & Holubkov, R. (1997a). Maximum likelihood estimation of logistic regression parameters
under two-phase outcome-dependent sampling. Journal of the Royal Statistical Society Series B, 59,
447–461.
Chen, K. (2001). Generalized case-cohort sampling. Journal of the Royal Statistical Society Series B, 63,
791–809.
Chen, K. & Lo, S. H. (1999). Case-cohort and case-control analysis with Cox’s model. Biometrika, 86,
755–764.
Cox, D. R. (1972). Regression models and life-tables (with discussion). Journal of the Royal Statistical
Society Series B, 34, 187–220.
Gor?ne, M., Zucker, D. M., & Hsu, L. (2009). Case-control survival analysis with a general semiparametric
shared frailty model: A pseudo full likelihood approach. The Annals of Statistics, 37, 1489–1517.
Gray, R. J. (2009). Weighted analyses for cohort sampling designs. Lifetime Data Analysis, 15, 24–40.
Hern´an, M. A., Brumback, B., & Robins, J. M. (2000). Marginal structural models to estimate the causal
effect of Zidovudine on the survival of HIV-positive men. Epidemiology, 11, 561–570.
Hsu, L. & Gor?ne, M. (2006). Multivariate survival analysis for case-control family data. Biostatistics, 7,
387–398.
Hsu, L., Chen, L., Gor?ne, M. & Malone, K. (2004). Semiparametric estimation of marginal hazard function
from case-control family studies. Biometrics, 60, 936–944.
Huo, T. I., Wu, J. C., Lin, H. C., Lee, F. Y., Hou, M. C., Lee, P. C., Chang, F. Y. & Lee, S. D. (2005).
Evaluation of the increase in model for end-stage liver disease (MELD) score over time as a prognostic
predictor in patients with advanced cirrhosis: risk factor analysis and comparison with initial MELD
and Child-Turcotte-Pugh score. Journal of Hepatology, 42, 826–832.
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
382
SCHAUBEL ET AL.
Vol. 42, No. 3
Kalb?eisch, J. D. & Lawless, J. F. (1988). Likelihood analysis of multi-state models for disease incidence
and mortality. Statistics in Medicine, 7, 149–160.
Kalb?eisch, J. D. & Prentice, R. L. (2002). The statistical analysis of failure time data, 2nd ed., Wiley,
Hoboken, NJ.
Kremers, W. K., van IJperen, M., Kim, W. R., Freeman, R. B., Harper, A. M., Kamath, P. S. & Wiesner,
R. H. (2004). MELD score as a predictor of pretransplant and posttransplant survival in OPTN/UNOS
status 1 patients. Hepatology, 39, 764–769.
Kulich, M. & Lin, D. Y. (2000). Additive hazards regression for case-cohort studies. Biometrika, 87, 73–87.
Li, H., Yang, P., & Schwartz, A. G. (1998). Analysis of age at onset from case-control family studies.
Biometrics, 54, 1030–1039.
Lin, D. Y. (2000). On ?tting Cox’s proportional hazards models to survey data. Biometrika, 87, 37–47.
Lin, D. Y. & Ying, Z. (1993). Cox regression with incomplete covariate measurements. Journal of the
American Statistical Association, 88, 1341–1349.
Merion, R. M., Schaubel, D. E., Dykstra, D. M., Freeman, R. B., Port, F. K. & Wolfe, R. A. (2005). The
survival bene?t of liver transplantation. American Journal of Transplantation, 5, 307–313.
Moger, T. A., Pawitan, A. and Borgan, Ø. (2008). Casde-cohort methods for survival data on families from
routine registers. Statistics in Medicine, 27, 1062–1074.
Nan, B., Kalb?eisch, J. D., & Yu, M. (2009). Asymptotic theory for the semiparametric accelerated failure
time model with missing data. Annals of Statistics, 37, 2351–2376.
Oakes, D. (1981). Survival times: aspects of partial likelihood (with discussion). International Statistical
Review, 49, 235–264.
Pan, Q. & Schaubel, D. E. (2008). Proportional hazards models based on biased samples and estimated
selection probabilities. Canadian Journal of Statistics, 36(1), 111–127.
Prentice, R. L. (1986). A case-cohort design for epidemiologic cohort studies and disease prevention trials.
Biometrika, 73, 1–11.
Prentice, R. L. & Brewlow, N. E. (1978). Retrospective studies and failure time models. Biometrika, 65,
153–158.
Robins, J. M. (1993). Information recovery and bias adjustment in proportional hazards regression analysis
of randomized trials using surrogate markers. Proceedings of the Biopharmaceutical Section, American
Statistical Association, 24–33.
Robins, J. M. & Finkelstein, D. M. (2000). Correcting for noncompliance and dependent censoring in an
AIDS clinical trial with inverse probability of censoring weighted (IPCW) log-rank tests. Biometrics,
56, 779–788.
Robins, J. M. & Rotnitzky, A. (1992). Recovery of information and adjustment for dependent censoring using
surrogate markers. In AIDS Epidemiology—Methodological Issues, Jewell N. P., Dietz K., & Farewell
V. T., editors. Birkh¨auser, Boston, pp. 297–331.
Rubin, D. B. (1977). Inference and missing data. Biometrika, 63, 581–592.
Schaubel, D. E., Wolfe, R. A., Sima, C. S. & Merion, R. M. (2009). Estimating the effect of a timedependent treatment by levels of an internal time-dependent covariate. Journal of the American Statistical
Association, 104, 49–59.
Schildcrout, J. S. & Heagerty, P. J. (2008). On outcome-dependent sampling designs for longitudinal binary
response data with time-varying covariates. Biostatistics, 9, 735–749.
Self, S. G. & Prentice, R. L. (1988). Asymptotic distribution theory and ef?ciency results for case-cohort
studies. Annals of Statistics, 16, 64–81.
Shih, J. H. & Chatterjee, N. (2002). Analysis of survival data from case-control family studies. Biometrics,
58, 502–509.
Song, R., Zhou, H., & Kosorok, M. R. (2009). On semiparametric ef?cient inference for two-stage outcomedependent sampling with a continuous outcome. Biometrika, 96, 221–228.
Subramanian, A., Sulkowski, M., Barin, B., Stablein, D., Curry, M., Nissen, N., Dove, L., Roland, M.,
Florman, S., Blumberg, E., Stosor, V., Jayaweera, D. T., Huprikar, S., Fung, J., Pruett, T., Stock, P. &
The Canadian Journal of Statistics / La revue canadienne de statistique
DOI: 10.1002/cjs
2014
CASE-CONTROL DATA WITH DEPENDENT CENSORING
383
Ragni, M. (2010). MELD score is an important predictor of pretransplantation mortality in HIV-infected
liver transplant candidates. Gastroenterology, 138, 159–164.
Wang, W., Scharfstein, D., Tan, Z. & MacKenzie, E. J. (2009). Causal inference in outcome-dependent
two-phase sampling designs. Journal of the Royal Statistical Society Series B, 71, 947–969.
Wiesner, R. H., McDiarmid, S. V., Kamath, P. S., Edwards, E. B., Malinchoc, M., Kremers, W. K., Krom,
R. A. & Kim, W. R. (2001). MELD and PELD: Application of survival models to liver allocation. Liver
Transplantation, 7, 567–580.
Zhou, H., Weaver, M. A., Qin, J., Longnecker, M. P. & Wang, M. C. (2002). A semiparametric empirical
likelihood method for data from an outcome-dependent sampling scheme with a continuous outcome.
Biometrics, 58, 413–421.
Received 23 December 2012
Accepted 31 March 2014
DOI: 10.1002/cjs
The Canadian Journal of Statistics / La revue canadienne de statistique
Copyright of Canadian Journal of Statistics is the property of Wiley-Blackwell and its content
may not be copied or emailed to multiple sites or posted to a listserv without the copyright
holder’s express written permission. However, users may print, download, or email articles for
individual use.