PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Am Stat Assoc. Author manuscript; available in PMC 2010 August 11.
Published in final edited form as:
J Am Stat Assoc. 2009 September 1; 104(487): 1241–1250.
doi:  10.1198/jasa.2009.tm08033
PMCID: PMC2920213
NIHMSID: NIHMS186503

A Class of Semiparametric Mixture Cure Survival Models with Dependent Censoring

Megan Othus, graduate student, Yi Li, Associate Professor, and Ram C. Tiwari, Associate Director

Abstract

Modern cancer treatments have substantially improved cure rates and have generated a great interest in and need for proper statistical tools to analyze survival data with non-negligible cure fractions. Data with cure fractions are often complicated by dependent censoring, and the analysis of this type of data typically involves untestable parametric assumptions on the dependence of the censoring mechanism and the true survival times. Motivated by the analysis of prostate cancer survival trends, we propose a class of semiparametric transformation cure models that allows for dependent censoring without making parametric assumptions on the dependence relationship. The proposed class of models encompasses a number of common models for the latency survival function, including the proportional hazards model and the proportional odds model, and also allows for time-dependent covariates. An inverse censoring probability reweighting scheme is used to derive unbiased estimating equations. We validate small sample properties with simulations and demonstrate the method with a data application.

Keywords: Estimating equation, proportional hazards model, proportional odds model, right censoring, transformation model

1. INTRODUCTION

Survival models incorporating a cure fraction are an important statistical tool for analyzing cancer survival data. Research for many cancer types, including prostate and breast cancer, non-Hodgkin’s lymphoma, leukemia, melanoma, and head and neck cancer, has shown that a significant proportion of patients with these cancers are cured after therapy; that is, an individual will have little or no risk of experiencing the event of interest, e.g. breast cancer, after treatment. Cure survival models are often developed assuming the study population is a mixture of cured and not cured subjects (Berkson and Gage 1952; Farewell 1982; Kuk and Chen 1992; Peng and Dear 2000; Sy and Taylor 2000). Recently, Lu and Ying (2004) used the mixture formulation to extend a class of semiparametric transformation models proposed by Cheng, Wie, and Ying (1995) to incorporate cure fractions. An alternative to the mixture cure formulation is the promotion cure model, which is popular for Bayesian cure survival analysis and has an attractive biological interpretation (Chen, Ibrahim, and Sinha 1999; Tsodikov, Ibrahim, and Yakovlev 2003). The promotion cure model has one covariate function that describes both the survival and cure components. In contrast, the two components of the mixture model can depend on different parameters, allowing for separate covariate interpretations for the cure function and the survival function for those who are not cured.

A key assumption made in the literature is independence between the survival time and the censoring time, but such an assumption is often violated in observational studies. For example, in our motivating study a non-negligible portion of patients diagnosed with prostate cancer was found to have died from cardiovascular disease. Research has shown that prostate cancer and heart disease share high fat diet as a common risk factor (Wilson et al. 1998; Strom et al. 2008). Therefore, assuming independence between the main end point (death from prostate cancer) and censoring causes (e.g. death from heart disease) is problematic. Dependent censoring will often hamper statistical analysis since classical models will not be valid, and ignoring dependent censoring will typically introduce bias (e.g. Li, Tiwari, and Guha 2007). To account for this, Li et al. (2007) developed a cure model in the presence of dependent censoring, but their model requires a priori assumptions on the dependence structure and does not allow for covariates.

In this project, we propose a semiparametric transformation model which allows for covariates as well as dependent censoring. The key idea is to use an inverse censoring probability reweighting scheme to derive unbiased estimating equations that account for dependent censoring. In this manner, we are able to avoid making parametric assumptions about the dependence structure between the survival time and the censoring time. It is also worth noting that our proposed model, which accommodates time-dependent covariates, is more general than that proposed by Lu and Ying (2004), which only allows for time-independent covariates.

Broadly speaking, this project advances the field in three distinct ways. First, the proposed methods can be used to investigate trends in Surveillance Epidemiology and End Results (SEER) cancer survival data (www.seer.cancer.gov); a number of cancer sites in the SEER database may indicate cure patterns, including prostate, colon, and ovarian cancer (Tai et al. 2005). Secondly, survival studies with cure fractions are subject to heavy censoring that is likely to violate the traditional independent censoring assumption. This project addresses the problem by developing methods that can accommodate dependent censoring, and by providing procedures for model diagnostics and techniques for statistical inference. The proposed estimation procedure successfully avoids parametric assumptions on the joint distribution of the censoring and true failure times, which have caused difficulties in the study of dependent censoring in the existing literature. Finally, the theoretical proofs developed in this project provide an original application of empirical process theory and functional analysis techniques, which will be valuable for semiparametric estimation.

The rest of the paper is stuctured as follows: in Section 2 we define notation and set up the survival model. In Section 3 we define the weight function. Section 4 uses the results of Sections 2 and 3 to define estimating equations and provides theoretical results. Section 5 summarizes a computational algorithm and a weighted bootstrap method for inference procedures on the estimating equations. Section 6 provides a model checking procedure. In Section 7 we present simulation results and data analysis findings. Section 8 contains a short discussion. Sketches of proofs of theorems are contained in the Appendix.

2. A CLASS OF SEMIPARAMETRIC MODELS

Denote the failure time for a patient by T, which is assumed to have the decomposition T = ηT* + (1 − η)∞, where η [set membership] {0, 1} indicates whether the sampled subject is cured (η = 0) or not (η = 1), and T*, with support on [0,∞), denotes the failure time if the patient is not cured. Let C denote the censoring time, X a κ length vector of covariates related to the cure indicator η, which includes 1, and Z(·) a ξ length vector of external time-dependent covariates related to T* that may share time-independent components with X. Define T = min(T,C), Δ = I(TC), a counting process related to the survival progress N(t) = I(Tt, Δ = 1), and the at-risk process Y (t) = I(Tt). In practice, we observe (Ti, Δi, Zi(·), Xi), for i = 1, …, n, which are independent and identically distributed copies of (T, Δ, Z(·), X). To ensure that we have unbiased estimates for our parameters we assume that the support of the censoring random variable is not shorter than the support of the failure time (Fine, Ying, and Wei 1998; Fine 1999).

Using the above notation, we can define a model for T. Following Farewell (1982), we use a logistic model for the cure indicator η:

P(η=1|X,Z)=G(γX),
(1)

where G(x) = exp(x)/(1 + exp(x)). Among patients who are not cured, we assume that given the covariate path Z(t) = {Z(u), ut}, T* possesses the following survival function:

Sg(t;Z¯(t))=g{0texp(βZ(u))dH(u)},
(2)

where H is a fixed but unspecified nondecreasing function with H(0) = 0, and g is a known continuous and strictly decreasing function such that g(0) = 1 and limt→∞ g(t) = 0. This model has been considered in Kosorok and Song (2007) and Zeng and Lin (2007). Zeng and Lin (2007) suggested extending their model to cure applications in future research; in contrast to our proposal, their work assumed conditional independence between T and C. The survival function for the noncured subjects, Sg(·), is commonly referred to as the latency survival function. The formulation of (2) is comprehensive: when g(x) = 1/(1 + x) the resulting model is a time-dependent construction of the proportional odds model (Bennett 1983), and when g(x) = exp(−x) the model reduces to the time-dependent proportional hazards model (Cox 1972).

Combining (1) and (2), we can describe the survival function for T as

P(Tt|Z¯(t),X)=g{0texp(βZ(u))dH(u)}G(γX)+G¯(γX),
(3)

where G(x) = 1 − G(x). Let W(·) denote the external time-dependent covariate process associated with censoring and let w(t) = {W(u), ut} denote the covariate path associated with W(·). We assume that P(Tt|Z(t), X) = P(Tt|Z(t), X, w (t)), i.e. that w (t) does not contain any extra information for T beyond what is already in Z(t) and X.

3. A REWEIGHTING SCHEME FOR DEPENDENT CENSORING

We wish to make minimal assumptions on the dependence between T and C, and so it will not be feasible to construct likelihood-based inference procedures. Instead, we will derive unbiased estimating equations for β, the parameter associated with the latency survival function, and γ, the parameter associated with the cure function. To this end, we first define some notation related to the censoring process. Define the crude hazard for censoring as

λc(t|W¯(t))=limdt0(dt)1P(C(t,t+dt),Δ=0|T˜t,W¯(t))
(4)

which can also be regarded as a cause specific hazard. Conveniently, (4) can be estimated with observed data even when T and C are dependent (Fleming and Harrington 1991, Theorem 1.3.1).

We make the following assumption on the crude hazard for censoring:

λc(t|W¯(t),T,T>t)=λc(t|W¯(t),T>t).

This assumption implies that, conditional on the covariate process w(t), the crude hazard for censoring at time t does not further depend on the possibly unobserved failure time. This assumption has been described as “no unmeasured confounders for censoring” and the assumption would fail if a covariate related to both T and C were not included in w(t) (Robins and Finkelstein 2000). This assumption has been well discussed for semiparametric survival problems in the literature with examples and information on how to select covariates for w(t) (Robins 1992; Robins and Rotnitzky 1992; Robins 1993; Rotnitzky and Robins 1995; Robins and Finkelstein 2000; Rotnitzky and Robins 2003; Miloslavsky, Keles, Laan, and Butler 2004).

Flexible models on λc(t|w(t)) can be considered. For example, if αE is the parameter vector associated with W(·), we could write the crude hazard function as

λc(t|W¯(t))=λ0(t)exp(αEW(t)),
(5)

where λ0(t) is left unspecified.

Let Λc(t|W¯(t))=0tλc(u|W¯(u))du and K(t|w(t)) = Πut (1 − dΛc (u|w(u))). Throughout we will write K (t|w(t)) = K (t). Note that if all the components of W(·) are time-independent and T and C are independent, then K(t) = P(C > t|W) and K(t) can be interpreted as the probability that a person with covariate W is not censored up to time t. This interpretation does not hold in general as the crude hazard λc(t) is not equal to the net hazard of censoring when T and C are dependent.

Similar arguments as in Robins (1993) lead to

E[I(T˜(t,t+dt),Δ=1)/K(T˜)]=P(T(t,t+dt))
(6)

and

E[I(T˜t)/K(t)]=P(Tt).
(7)

These equations imply that the reweighted sample of patients observed to fail at time t is representative of all patients who fail at time t and that the reweighted sample of patients who have not failed by time t is representative of the general survival probability of all patients at time t.

4. ESTIMATING EQUATIONS AND THEORETICAL RESULTS

Equations (6) and (7) are important as they hold even when T and C are dependent, and we can use them to derive estimating equations for the unknown regression parameters β and γ, and the function H. Specifically, (6) implies that E[dN(t)/K(t−)+dP(Tt|Z(t), X)] = 0; or, equivalently, EW[E[dN(t)+K(t−) dP(Tt|Z(t), X)|W]] = 0; while (7) implies that EW [E[Y(t) − K(t−)P(Tt|Z(t), X)|W]] = 0: Together these constitute a basis for estimation. Here, P(Tt|Z(t), X) is specified in Equation (3) and dP(·) is the derivative of Equation (3) with respect to t.

Let α denote the parameters associated with K(·); define θ = (α, β, γ, H) and write the true value of θ as θ0 = (α0, β0, γ0, H0). Denote the estimating equations associated with α as Unα (α). In the following, the subscript i marks individuals. Define

μβ,i(θ)=0νZi(u)[dNi(u)+Ki(u)G(γXi)dSg(u;Z¯i(u))]μγ,i(θ)=0νXi[Yi(u)Ki(u)[G(γXi)Sg(u;Z¯i(u))+G¯(γXi)]]duμH,i(θ)(t)=0tdNi(u)+Ki(u)G(γXi)dSg(u;Z¯i(u)).

A natural set of estimating equations for β, γ and H using (3) and our results from (6), and (7) is:

Unβ(θ)=defi=1nμβ,i(θ)
(8)

Unγ(θ)=defi=1nμγ,i(θ)
(9)

UnH(θ)(t)=defi=1nμH,i(θ)(t),t[0,ν].
(10)

In the above, choose ν to be a constant such that P(T > ν) > č > 0 for some fixed č. The constraint on ν is imposed to avoid tail instability associated with large values of t. In practice the weight function K(·) is unknown and can be replaced with an empirical estimator. This is easily done by replacing Λc(u|w(u)) in the definition of K(t) with a consistent estimate [Lambda with circumflex]c(u|w(u)). Since the left term of UnH only jumps at observed failure times, the estimate for H will be a step-function that jumps at the observed failure times. Let Un(θ) = (Unα, Unβ, Unγ, UnH)(θ) and define U(θ) = EUn(θ).

To make the definition of Unα more concrete, we provide an example where the crude hazard for censoring (Equation (4)), follows a proportional hazards model (Equation (5)). In this case we can decompose α into two parts (αE, αΛ) where αE represents the Euclidean parameter of the model and αΛ denotes the baseline cumulative hazard function for the model. The estimating equations associated with (αE, αΛ) are the usual partial likelihood score equations and the Breslow formula for the cumulative hazard:

UnαE(α)=i=1n0νWi(u)(dNci(u)Yi(u)exp(αEWi(u))dΛ0(u))UnαΛ(α)(t)=i=1n0tdNci(u)Yi(u)exp(αEWi(u))dΛ0(u),t[0,ν]

where Nci(t) = I(Tit, Δi = 0) for i = 1, …, n, and Λ0(·) is the baseline cumulative hazard function. For this proportional hazards case we would write Unα = (UnαE, UnαΛ)′.

Define [theta w/ hat] = ([alpha], [beta], [gamma with circumflex], Ĥ) as zeroes of Un(θ) = (Unα, Unβ, Unγ, UnH)(θ). In general one is interested in the zeroes of (8) and (9), which we denote [var phiv with circumflex] = ([beta], [gamma with circumflex]). We note that the asymptotic variance formula for n1/2([var phiv with circumflex][var phi]) is analytically so complicated that it will be of limited computational utility. Instead, we propose an application of a weighted bootstrap, which is introduced in Section 5 and is shown to yield a consistent estimate for the variance. We will write the weighted bootstrap version of Un(θ) by Un(θ)* and denote the weighted bootstrap estimate of θ from Un(θ)* as [theta w/ hat]*.

Conditions for and sketches of proofs for the following asymptotic results are contained in the Appendix.

Theorem 1

Under regularity conditions (c.1)–(c.4), the survival model in Equation (3) is identifiable.

This result is more general than that for cure survival models provided by Li, Taylor, and Sy (2001) since time-dependent cure proportional odds and proportional hazards models are included in our result as special cases. Further, in contrast to Li et al. (2001), we do not require the value zero to be in the support of the covariates, which may not be a realistic assumption in many situations. As pointed out by a reviewer, a key component of the identifiability result is the presence of covariates in both components of the mixture model. In fact, Li et al. (2001) show that if the latency survival function is nonparameteric and independent of covariates and if the cure fraction is taken to be a constant (and so does not depend on covariates) the mixture cure survival model is not identifiable. The dependence of our result on covariates is analogous to identifiability results in the competing risks literature. A number of authors have put forward identifiable competing risks models; and all of the results depend on the presence of covariates in the models (Heckman and Honore 1989; Abbring and van den Berg 2003; Lee 2006).

Theorem 2

Under regularity conditions (c.1)–(c.7), [theta w/ hat] is a consistent estimator for θ0.

Theorem 3

Under regularity conditions (c.1) – (c.7) and (b.1)–(b.5), n1/2([theta w/ hat]θ0) converges weakly to a zero-mean Gaussian process. Conditional on the observed data, n1/2([theta w/ hat]* − [theta w/ hat]) has the same limiting distribution as n1/2([theta w/ hat]θ0).

5. INFERENCE PROCEDURES AND WEIGHTED BOOTSTRAP SCHEME

5.1. Algorithm

With the large sample theory established, we propose an iterative algorithm to solve estimating equations (8)(10). To start, let t0 < t1 < … < tk be the ordered observed failure times such that tk ≤ ν and t0 = 0. In practice, ν can be set equal to the 90th or 95th percentile of the observed failure times.

  • Step 1: Choose initial values for β, γ, and H, denoted [beta](0), [gamma with circumflex](0), Ĥ(0), respectively. In general, the initial values for β and γ can be obtained from estimates that ignore dependent censoring or potential cure fractions and one can let Ĥ(0) be the constant value 0.
  • Step 2: Compute estimates for K(·). For example, if the crude censoring distribution follows a time-dependent proportional hazards model, the calculation can be done following Robins (1993).

    Suppose that we have estimates for β, γ, and H from the (m − 1)th iteration; denote these estimates by [beta](m − 1), [gamma with circumflex](m − 1) and Ĥ(m − 1).

  • Step 3: Recall that Ĥ(·) is a step function with jumps at tj, j = 1, … k. Denote the jump size at tj as Ĥ{tj}. For j = 1, …, k obtain estimates Ĥ(m){tj}, in order, as solutions to UnH(tj) = 0 with β and γ set equal to [beta](m − 1) and [gamma with circumflex](m − 1), respectively. As Ĥ(m)(t) is a step function, 0tjeZi(u)β^(m1)dH^(m)(u)=l=1jeZi(tl)β^(m1)H^(m){tl}. Hence Ĥ(m) {tj} can be obtained recursively given the estimates of Ĥ(m) {tl}, lj − 1. A unique solution is guaranteed to exist due to the monotonicity of g.
  • Step 4: To estimate [beta](m), solve for β such that Unβ = 0 with γ set equal to [gamma with circumflex](m−1) and H set equal to Ĥ(m). The proof of Theorem 2 establishes that the derivative of Unβ is negative definite, therefore the Newton-Raphson solution will be unique.
  • Step 5: To find [gamma with circumflex](m), solve for γ such that Unγ = 0, with β = [beta](m) and H = Ĥ(m). As with the solution [beta](m), the derivative of Unγ is negative definite (see the proof of Theorem 2), so the Newton-Raphson method can be used to find a unique solution.
  • Step 6: Repeat Steps 3 – 5 until predetermined convergence criteria are met.

5.2. Weighted Boostrap

To conduct inference on the parameter estimates we use a weighted bootstrap (Wellner and Zhan 1996), which is applicable even in this case with an infinite dimensional nuisance parameter. Define μαi to be such that i=1nμαi are the estimating equations for α. Recall that the estimates for θ = (α, β, γ, H) are approximate zeros of Unj=i=1nμji for j = α, β, γ, H. Choose a positive distribution B and generate (bn1, …, bnn) such that B and (bn1, …, bnn) meet conditions (b.1) – (b.5) in the Appendix. Denote the zeroes of i=1nμjibni, j = α, β, γ, H as [theta w/ hat]* = ([alpha]*, [beta]*, [gamma with circumflex]*, Ĥ*). As stated in Theorem 3, the conditional distribution of n1/2([theta w/ hat][theta w/ hat]*) is asymptotically equivalent to the unconditional distribution of n1/2([theta w/ hat]θ0). We can approximate the distribution of [theta w/ hat]* by generating a large number of samples, (bn1,…, bnn), from B and for each sample calculating a realization of [theta w/ hat]*. The empirical distribution can be used to compute statistics for inference.

6. MODEL DIAGNOSTICS

In practice, it is of substantial interest to verify whether the functional forms imposed on the regression models are plausible. To this end we propose a two-step model checking procedure. First, assumptions on the form of K need to be verified. If the assumed model for K is correct, the asymptotic distribution of Ξ(t)=n1/2i=1nq˜(Wi,t)μαi(α^) will be a zero-mean Gaussian process, where q(·, ·) is a weight function that is assumed to be known and bounded. Furthermore, if we let (bn1,…,bnn) and [alpha]* be defined as in the previous section, we can approximate the distribution of Ξ(t) with Ξ˜(t)=n1/2i=1nq˜(Wi,t)(μαi(α^)(1bni)μαi(α^)+μαi(α^*)). That the distribution [Xi](t) approximates Ξ(t) follows from the argument in Appendix C of Peng and Huang (2008).

A formal lack-of-fit test can be based upon the statistic TK = ||Ξ (t)||, where ||· || denotes the supremum norm on [0,ν]. We can generate a large number of sample paths from [Xi](·) by repeatedly sampling new (bn1,…, bnn). A p-value for TK can be estimated by the empirical proportion of Tk||[Xi](t)||. Similar arguments as in Lin, Wei, and Ying (1993) show that this test is consistent against the alternative that the model specification of K is false. A corresponding visual diagnostic method is found by plotting a number of realizations of [Xi](·) and observing how extreme the pattern of Ξ(·) looks in comparison to the simulated paths.

Once we are satisfied with the fit of K we can check the form of (3). Define Mi(t, θ) = Yi(t)/Ki(t−) −G(γXi)Sg(t, Z(t)) − G(γXi). Consider the function Γ(t)=n1/2i=1nw(Zi,Xi,t)Mi(t,θ^), where the weight function w(·, ·, ·) is assumed to be known and bounded. If (3) is correctly specified and K(·) is consistently estimated, Γ(t) converges to a zero-mean Gaussian process. We can approximate the distribution of Γ(t) with the distribution of

Γ˜(t)=n1/2i=1nw(Zi,Xi,t)(Mi(t,θ^)(1bni)Mi(t,θ^)+Mi(t,θ^*)),
(11)

where bni and [theta w/ hat]* are defined as in the previous section. If M(·, θ0) is specified correctly, Γ(t) should have mean zero for any t. Following the arguments for Ξ we can establish that the distribution of Γ̃ is asymptotically the same as the distribution of Γ.

As in the model-checking procedure for K(·), a formal test for lack-of-fit can be made based upon TS = ||Γ(t)||. As with the diagnostic method for K(·), we can calculate a p-value for TS (the corresponding test can be shown to be consistent against the general alternative that (3) is not correctly specified) and develop a goodness-of-fit plot.

7. NUMERICAL STUDIES

7.1. Simulations

We evaluated the finite sample performance of the proposed methods with simulations. Two such studies are summarized in Table 1. We considered proportional odds and proportional hazards models. Dependent censoring can occur when a covariate that affects both the survival and censoring times is not included in the model. To mimic this situation we generated two covariates for each subject: a Bernoulli(0.5) random variable and a Uniform(0,1) random variable. The censoring times followed a proportional hazards model with both the Bernoulli and Uniform random variables as covariates. The latency survival times followed either the proportional hazards or proportional odds models with both the Bernoulli and Uniform random variables as covariates. The cure functions only depended on the Bernoulli random variable. For the proportional hazards model we took H(t) = t and for the proportional odds model we took H(t) = 2t.

Table 1
Simulation results

In all of the analyses we only included one covariate in the cure function and in the latency survival function, the Bernoulli random variable. To model the crude hazard for censoring we used a proportional hazards model with the Bernoulli random variable as the only covariate. Bias, monte carlo standard error estimates, bootstrap based standard error estimates, and estimates of the root mean squared error (RMSE) for several different parameter configurations are presented in Table 1. We found the “true” values of the misspecified latency survival parameters by generating a very large dataset (sample size = 1,000,000) with the two covariates, but fitting the model on the dataset with just one covariate.

To investigate the biases that occur when dependent censoring is ignored we compared our method to the method of Lu and Ying (2004) (we will refer to this method as LY hereafter), which does not compensate for dependent censoring. Summaries for bias, monte carlo standard errors, and estimates of RMSE can also be found in Table 1. The proposed method has the most significant gains over LY in estimating the parameter associated with the cure proportion. The proposed method and LY have comparable bias in estimating the latency survival parameter. In practice, if correct estimates of the cure proportion are of primary interest, the proposed method can reduce bias even when accounting for only some of the dependent censoring.

7.2. SEER Prostate Cancer Data

We applied the proposed method to prostate cancer data from the SEER database, using data released in 2004. Specifically, we looked at males diagnosed with prostate cancer from the Detroit, Michigan area whose cancer stage was classified as local or regional. The Detroit registry was one of the original nine SEER registries and so has longer follow-up than the more recent registries. To avoid potential confounding related to whether there was any prior cancer diagnosis, we restricted our sample to those men whose prostate cancer was their first cancer diagnosis. The Federal Drug Administration (FDA) approved the prostate-specific antigen (PSA) test for use in monitoring patients diagnosed with prostate cancer in 1986. The use of the PSA test for screening of prostate cancer increased dramatically starting in 1990, and so to avoid potential biases associated with the new screening method we only considered cases diagnosed up to 1989. Our random variable of interest was time from diagnosis of prostate cancer to death from prostate cancer.

One question of interest was whether survival or cure fractions differed in this dataset by race. Due to the small number of subjects in some racial categories in this dataset we restricted our attention to the two largest racial groups, black and white. There were 15,524 cases in the SEER database that met the criteria; 3,837 were classified as black and 11,687 were classified as white.

In this cohort, 3,717 people died of prostate cancer while 4,349 people died of heart disease. These numbers are summarized in Table 2. Heart disease and prostate cancer share the common risk factor of high fat diet (Wilson et al. 1998; Strom et al. 2008). High fat diets have been shown to increase androgen levels in the body, increase cholesterol levels, and are associated with increased consumption of Omega-6 fatty acids. Increased levels of androgens have been shown to be associated with prostate cancer cell proliferation (Fleshner et al. 2004; Pandini et al. 2005; Xu et al. 2006; Asirvatham et al. 2006). High cholesterol and Omega-6 fatty acids have been shown to be associated with prostate carcinogenesis (Homma et al. 2004; Hughes-Fulford et al. 2005; Zhuang et al. 2005; Bravi et al. 2006; Ritch et al. 2007). Further complicating the relationship between heart disease and prostate cancer is the fact that many men diagnosed with prostate cancer are treated with androgen therapy to reduce hormone levels. This treatment has been shown to be associated with an increased risk of death from heart disease (Keating et al. 2006). Given these connections and the high rate of death from heart disease in this population, there is some evidence that there may exist dependent censoring in this dataset.

Table 2
Summary of causes of death among SEER data

Figure 1 shows the naïve (assuming independent censoring) Kaplan-Meier survival plots for the two subgroups. We can observe that both curves plateau above 0.5, indicating that a large proportion of the population may not die of prostate cancer. While for the first several years the survival curves are indistinguishable, the white subgroup’s survival tail is noticeably higher than the tail of the black group.

Figure 1
Naïve Kaplan-Meier curves (ignoring potential dependent censoring) for survival comparing white subgroup (dashed line) to the black subgroup (bold line).

There are several situations when physicians consider patients to be cured of prostate cancer. First, surgery (radical prostatectomy) and radiotherapy (external-beam radiotherapy, brachytherapy, or both) in early stage prostate cancer patients can entirely eradicate the disease, inducing a cure among a significant proportion of patients (Hanlon and Hanks 2000; Jani and Hellman 2003). Also, improvements in prostate cancer treatment have made it so that many prostate cancer patients who are not cured are treated for chronic disease rather than acute disease, and their probability of dying of prostate cancer is not elevated above that of the general population (American Cancer Society 2007). Given these considerations, a model that explicitly accounts for a cure proportion and allows for inference on the cure proportion can be of use for this prostate cancer application.

We considered both proportional hazards and proportional odds models. Race was a covariate for the survival and cure functions, defined RS and RC, respectively. The models also contain a time-dependent covariate written PSAS, which takes the value one after the year 1985 and the value zero otherwise. This covariate can help to account for changes in survival induced by the policy change in 1986 when the FDA approved the PSA test to help monitor progress of prostate cancer. We modeled the crude hazard for censoring with a proportional hazards model containing PSAS, RS, and covariates for grade, year of diagnosis, and age at diagnosis. We took B ~Exponential(1) and standard error estimates are based on 1000 bootstrap replicates.

Figure 2 provides examples of model-checking plots for the time-dependent proportional hazards and proportional odds cure models. The plot of Equation (11) against time for the proportional odds cure model does not appear to have mean zero and the p-value from 1000 bootstrap replicates with B ~ Exponential(1) is 0.008, which indicates that the proportional odds cure model may not be appropriate. The corresponding plot for the proportional hazards cure model shows an approximately mean zero trend and the p-value from 1000 bootstrap replicates with B ~ Exponential(1) is 0.416. These data suggest the proportional hazards cure model fits the data better than the proportional odds cure model, so we will present results for the proportional hazards cure model. Results for the time-dependent model and a univariate model containing only RS and RC can be found in Table 3.

Figure 2
Model checking plots. Twenty-five random sample paths are plotted in grey with the observed paths marked in black.
Table 3
SEER data results for a proportional hazards cure model. Parameter estimates and associated standard errors from the regression models are provided. Estimates with the subscript S can be interpreted as the log-hazard ratio for non-cured subjects. Estimates ...

The race covariate is significant in the univariate model for both the survival and cure functions. The signs of the estimates indicate that black race is associated with worse survival and a smaller probability of being cured. These covariates remain significant in the time-dependent model, though the time-dependent covariate indicating PSA approval is not significant.

We contrast our univariate results with alternative methods in Tables 4 and and5.5. All models contain one covariate denoting race. The standard proportional hazards model (ignoring possible cure fractions and dependent censoring) estimate is almost identical to the competing risks proportional hazards model estimate with deaths from other causes treated as censored events (Fine and Gray 1999). The three parameter values presented in Table 4 have different interpretations. The parameter estimate from our model has the interpretation of the log-hazard ratio for subjects who are not cured. The standard proportional hazards model parameter has the interpretation of the log-hazard ratio assuming the population is composed entirely of persons who are not cured. The log-hazard ratio for the subdistribution for death from prostate cancer are not easily interpretable, as noted in Fine and Gray (1999).

Table 4
Parameter estimates for proportional hazards functions for three different models. Black race is the reference category for each of the models.
Table 5
Estimates of the proportion of men cured of prostate cancer: (1) from the tail of the Kaplan-Meier survival curve stratified by race and (2) from transforming the regression results from our proportional hazards cure model.

Under the assumption of independent censoring, the tail end of the Kaplan-Meier curve can estimate the proportion of cured subjects in a sample (Maller and Zhou 1992). Table 5 compares estimates of the tail of the Kaplan-Meier survival function stratified by race to estimates of the cure proportion found from transforming the regression results from the our proportional hazards cure model. The proportional hazards cure model has a smaller cure fraction than the tail end of the Kaplan-Meier plot indicates, and the confidence intervals for the estimates do not overlap.

As a referee pointed out, competing risks models (Prentice et al. 1978; Fine and Gray 1999) are a useful way to model survival data with the potential for dependent censoring. However, this type of analysis may not be ideal for this specific application. First, competing risks models can characterize the crude risk of failure for a specific failure type, but do so without removing dependent risks (Pepe and Mori 1993). We wish to evaluate the effect of race on prostate cancer survival, but these competing risks models will not allow for the direct covariate interpretation we are interested in. Second, a subject who is cured of prostate cancer could die of heart disease (a competing risk). Therefore the probability of being cured is not readily modeled using the existing competing risk framework, in contrast to the model we proposed.

8. DISCUSSION

This paper provides a unified framework for a class of transformation cure models that allow for dependent censoring. We utilize an inverse censoring probability formula to develop unbiased estimating equations. Our method allows one to estimate cure fractions more accurately in cases where there is dependence between censoring and survival times. We also develop a weighted bootstrap for inference procedures and supply a method to check whether the model specification is adequate.

Extending our estimation procedure to be able to incorporate internal time-dependent covariates is a subject of future research. This extension would allow for a greater range of application in dependent censoring situations. A common cause of dependent censoring is excluding covariates related to both the survival and censoring times. This is an issue in any model, including models with only time-independent covariates. As shown in our simulation study, our method can decrease bias in situations with only time-independent covariates.

Rigorous testing for whether there is sufficient follow-up to identify cure proportions is difficult even without dependent censoring. A test for sufficient follow-up in the presence of dependent censoring may deserve a separate paper. A number of tests for sufficient follow-up in cure survival situations have been put forward (Maller and Zhou 1994, 1995; Klebanov and Yakovlev 2007). These tests, combined with subject matter specific knowledge on whether a mixture cure formulation is suitable, can provide information to allow an investigator to decide whether it is appropriate to use this model.

Acknowledgments

The work of Ram Tiwari was conducted while he was at the National Cancer Institute (NCI). The views expressed in this article are his own and do not necessarily represent those of the NCI or the FDA. This work was supported in part by National Institute of Health Grants CA09337-25 and ES07142-24. The authors thank Eric (Rocky) Feuer and Angela Mariotto for their extensive and insightful comments on this manuscript.

APPENDIX: TECHNICAL DETAILS

Before we present the asymptotic results we will review and define some notation. All expectations are with respect to the true distributions of all random variables involved. Denote the Euclidean norm as || · ||, the supremum norm on [0, ν] as || · ||, and the norm corresponding to the parameters of α as || · ||α. For example, if we assumed the crude hazard for censoring followed a proportional hazards form, our model would have two parameters: αE, a Euclidean parameter associated with W(·), and αΛ, associated with the cumulative baseline hazard. In this case a sensible choice for ||αα0||α would be ||αEα0E|| + ||αΛα||. Define ||θθ0|| = ||αα0||α+||ββ0||+ ||γγ0|| + ||HH0||. Let a single superscript dot denote derivatives and a double superscript dot denote second derivatives.

Consider the following regularity conditions:

  • (c.1) the covariates Xi and the covariate processes Zi(·) and Wi(·) are uniformly bounded, have nondegenerate variance-covariance matrices, and are cadlag for i = 1, …, n;
  • (c.2) [var phi]0 = (β0, γ0) is a point in a bounded, open, convex subset of Rξ+κ;
  • (c.3) the transformation function g is twice differentiable. The function g, its derivative ġ, and its second derivative g are all Lipschitz;
  • (c.4) H0(t) is differentiable and is of bounded variation over [0, ν];
  • (c.5) E(n−1 Un[var phi](θ)|[var phi]=[var phi]0) is nondegenerate with probability 1;
  • (c.6) max Ni(ν) > 0, i = 1,…, n;
  • (c.7) K(·) is uniformly bounded away from 0 on [0, ν]. The estimating equation Unα meets conditions 1, 3, and 4 of Theorem 2. Unα is Fréchet differentiable at α0 and the derivative is continuously invertible;
  • (b.1) the vectors (bn1,…, bnn)′ are exchangeable for all n = 1, 2,…;
  • (b.2) bni ≥ 0 for all n and all i and i=1nbni=n for all n;
  • (b.3) the L2,1 norm of bn1 is bounded, i.e. 0PB(bn1u)duM*<;
  • (b.4) limψ→∞ lim supn→∞t≥ψ t2 P(bn1t) = 0;
  • (b.5) n1i=1n(bni1)2c2>0 in PB-probability.

Conditions (c.1) – (c.5) are standard assumptions in survival analysis literature. (c.6) is used in the proofs of Theorems 1 and 2. (c.7) requires that the estimates for the weight function, K, meet the same regularity conditions that we impose on β, γ, and H. These conditions have been previously shown for a variety of survival functions (e.g. van der Vaart 1998, Section 25.12 for the proportional hazards model).

Conditions (b.1) – (b.5) are met by most positive distributions by defining bni=bi/b¯n, where bi is an observed realization from the bootstrap distribution, B, and b′ denotes the average of n observed realizations. If B ~ exp(1), B meets all the regularity conditions.

The following are sketches of the proofs. More complete details can be found in a technical report available from the first author.

Proof of Theorem 1

We want to show that

G¯(γX)+G(γX)g{0texp(βZ(u))dH(u)}=G¯(γ˜X)+G(γ˜X)g{0texp(β˜Z(u))dH˜(u)}
(12)

for all t [set membership] [0, ν] with probability 1 implies that γ = [gamma with tilde], β = beta, H = H. Equation (12) implies that there exists a time-independent constant v* (possibly random) such that

g{0texp(β˜Z(u))dH˜(u)}1g{0texp(βZ(u))dH(u)}1=υ*

for all t [set membership] [0, ν]. Rewrite beta = β + [var pi] and H˜(t)=H(t)+a0tω(u)dH(u) for some nonrandom function ω(u) and a [set membership]R. Express υ* as

g{0texp(β˜Z(u))dH˜(u)}1g{0texp(βZ(u))dH(u)}1=g{0t(1+aω(u))exp(Z(u)(β+ϖ))dH(u)}1g{0texp(βZ(u))dH(u)}1
(13)

for all t [set membership] [0, ν]. It follows that (13) is constant with respect to t if and only if [var pi] = 0 and a = 0. To see that this is the case, use the mean value theorem and write the right hand side of (13) as

1+g¯(t)0texp(βZ(u))[Z(u)ϖ+a]dH(u),
(14)

where

g¯(t)=g˙{0t{1+a¯ω(u)}exp(β¯Z(u))dH(u)}g{0texp(βZ(u))dH(u)}1,

[beta with macron above] = β + sbeta and ā = sa for some constant s [set membership] [0, 1]. Equation (14) is independent of t if and only if [var pi] = 0 and a = 0. Hence beta = β and H = H. Therefore, a reëxamination of (12) yields that G(γX) = G([gamma with tilde]X) with probability 1. Since the logistic function G(·) is invertible, it follows that γX = [gamma with tilde]X or (γ[gamma with tilde])′X = 0 with probability 1. As X has a nondegenerate variance-covariance matrix,γ = [gamma with tilde].

Proof of Theorem 2

First consider the convergence of Ĥ (t, [var phi]0, α0), which is solved recursively from (10) with [var phi] and K(·) fixed at the true values, [var phi]0 and K0(·). Let H be the set of all nondecreasing, nonnegative, finite step functions with H(·) ϵ H such that H(0) = 0 and H(·) jumps only at the observed failure times from the data, t1 < … < tk ≤ ν. We want to show that ||Ĥ(t, [var phi]0, α0) − H0(t)|| → 0 almost surely as n → ∞. Define the mapping Q on H as

Q(H)(t)=i=1n0tdNi(u)/Ki(u)+G(γ0Xi)dSg(Z¯i(t),t,β0).

For an arbitrary but fixed ε′ > 0 consider H1 ϵ H and H2 ϵ H such that ||H1H2|| ≥ ε′. Denote the jump size of Hl at tj as Hl{tj}, l = 1, 2. Write Q(H1)(t) − Q(H2)(t) = ∑tjt wjaj where wj=i=1nG(γ0Xi)g˙(si)exp(β0Zi(tj)),0 ≤ si|0texp(β0Zi(u))dH1(u)0texp(β0Zi(u))dH2(u)|, and aj=H1{tj} − H2{tj}. Since ||H1(t) − H2(t)|| ≥ ε′/2, there must exist tj0 such that aj0 = |H1{tj0} − H2{tj0}| ≥ ε′/2k. Note that since Z(·) and X are uniformly bounded there must exist some finite c such that n/cwjnc for j = 1,…, k.

Suppose that |∑tjtj0 wjaj| ≤ ε′n/4kc. It follows that |∑tjtj0−1 wjaj | ≥ |wj0aj0| − ε′n/4kc ≥ ε′n/2kc − ε′n/4kc = ε′n/4kc. Therefore, either |∑tjtj0 wjaj| ≥ ε′n/4kc or |∑tjtj0−1 wjaj| ≥ ε′n/4kc. Hence

i=1nG(γ0Xi)[Sg(t,Z¯i(t),β0,H1)Sg(t,Z¯i(t),β0,H2)]εn/4kc¯ε/4c¯.

Select H[set membership] H such that H′ (tj) = H0(tj) for i = 1,…, k. The law of large numbers and the continuity of H0 imply that ||Q(H′)(t)|| → 0 almost surely as n → ∞. We defined Ĥ(t, ϕ0, α0) such that Q(Ĥ(t, ϕ0, α0)) = op(1) for all t [set membership] [0, ν], which implies that ||Q(H′)(t) − Q{Ĥ(ϕ0, α0)}(t)|| → 0 almost surely as n → ∞. Therefore Ĥ(·, ϕ0, α0) is in an ε′ neighborhood of H′ in [0, ν] under || · || with probability 1. Using the facts that Ĥ(·, ϕ0, α0) and H0 are monotone, and that ε′ > 0 can be made arbitrarily small, we conclude that ||Ĥ(·, ϕ0, α0) − H0(·) || → 0 in [0, ν] almost surely.

We apply the inverse function theorem to prove consistency for [var phi] (Foutz 1977). That the estimating equations for β and γ are asymptotically unbiased follows from the consistency of [K with circumflex](·) and the convergence of Ĥ(·, [var phi]0, α0). Tedious but straight-forward calculations can show that En−1 Un[var phi](θ)|ϕ=[var phi]0 converges uniformly in a neighborhood of the true parameters and is negative definite. Hence ϕ is consistent. Finally, the consistency of Ĥ(·, θ) follows immediately from the convergence of H(·, ϕ0, α0), the consistency of ϕ and K (consistency of K is assumed in (c.7)), and the continuity of Ĥ(·, ϕ) with respect to ϕ and K(·).

Proof of Theorem 3

We use Theorem 3.3.1 and Lemma 3.3.5 of van der Vaart and Wellner (1996) and Theorem 3.1 from Wellner and Zhan (1996). In the following, define Uj = Eμj for j = α, β, γ, H and U = (Uα,…,UH), where the μj are defined as in Section 4. Let Hβ and Hγ to be finite sets with same cardinality as the parameters and let HH = {ht = 1[0,t] : t [set membership] [0, ν]}. The following four conditions need to be verified:

  1. For j = α, β, γ, H and for any δn → 0: sup{n(UnjUj)(θ)n(UnjUj)(θ0)j/1+nθθ0:θθ0δn}=op(1), where || · ||j denotes the norm associated with parameter j.
  2. U is Fréchet differentiable at θ0 with derivative denoted U0, and U0 has a continuous inverse.
  3. [theta w/ hat] is consistent for θ0 and Un([theta w/ hat]) = op(n−1/2).
  4. If Dnj is the envelope function of the class 𝔻nj{μj(θ)(h)μj(θ0)(h)1+nθθ0:h𝓗j,θθ0<δn}, for every δn → 0 and j = α, β, γ, H,
    limψlim supnsuptψt2P(D˜nj(X,W,Z,Δ)>t)=0.

Conditions 3 and 4 clearly hold in this situation. Theorem 2 and (c.5) establish the consistency of [theta w/ hat], and the second requirement of Condition 3 holds by definition of [theta w/ hat]. Let [M with tilde] = max(||Z(t)||, ||X||, 2) and let Dnβ = Dnγ= DnH = 2[M with tilde]; it follows that Condition 4 is met for j = β, γ, H. Condition 4 is met for j = α by (c.7).

For j = β, γ, H, Condition 1 can be shown using standard arguments by building up Donsker classes of terms (Parner 1998; van der Vaart and Wellner 1996) and by using the continuity of K, G, and g and the pointwise convergence of θ to θ0. Condition 1 holds for j = α by (c.7).

The last step is to verify Condition 2. Write the Fréchet derivative, U0, of U at θ0 as a partitioned matrix: (BCDE), where B corresponds to the the Fréechet derivative of Uα with respect to α evaluated at θ0; C to the Fréchet derivative of Uα with respect to β, γ and H evaluated at θ0; D to the Fréchet derivative of Uβ, Uγ, UH with respect to α evaluated at θ0; and E to the Fréchet derivative of Uβ, Uγ, UH with respect to β, γ and H evaluated at θ0. The inverse of this partitioned matrix can be written as (B1+B1CF1DB1B1CF1F1DB1F1) where F = E − DB−1C (Stapleton 1995). Therefore, to prove that the partitioned matrix is continuously invertible, one only needs to show that B and F are continuously invertible. Since Uα is not a function of β, γ or H, C is the zero matrix and we only need to show that B and E are continuously invertible. B is continuously invertible by (c.5).

To verify that E is continuously invertible, write E as a partitioned matrix and use the same argument as above. Specifically, write: E=(U˙pU˙qU˙rU˙s), where Up is the derivative of (Uβ, Uγ) with respect to β and γ evaluated at θ0; Uq is the derivative of (Uβ, Uγ) with respect to H evaluated at θ0; Ur is the derivative of UH with respect to β and γ evaluated at θ0; and Us is the derivative of UH with respect to H evaluated at θ0. We need to show that Up and U˙sU˙rU˙p1U˙q are continuously invertible. Up is finite dimensional and so continuous invertibility is a consequence of it being one-to-one. Us can be shown to be continuously invertible by direct computation, and U˙rU˙p1U˙q can be shown to be compact by invoking Helly’s selection theorem and the dominated convergence theorem. Direct calculation can show that (U˙sU˙rU˙p1U˙q)[A]=0 if and only if ||A|| = 0. We can therefore conclude that U˙sU˙rU˙p1U˙q is one-to-one. The sum of a continuously invertible operator and a compact operator is continuously invertible if the sum is one-to-one (Rudin 1987). Therefore U˙sU˙rU˙p1U˙q is continuously invertible and the result follows.

Contributor Information

Megan Othus, Department of Biostatistics, Harvard University and Department of Biostatistics and Computational Biology, Dana Farber Cancer Institute, Boston, MA, 02115 (ude.dravrah.saf@suhtom)

Yi Li, Department of Biostatistics, Harvard University, Department of Biostatistics and Computational Biology, Dana Farber Cancer Institute, Boston, MA, 02115 (ude.dravrah.ymmij@iliy)

Ram C. Tiwari, Statistical Science and Policy, Office of Biostatistics, Center for Drug Evaluation & Research, Federal Drug Administration (FDA), Silver Springs, MD, 20993 (vog.shh.adf@irawit.mar)

References

  • Abbring J, van den Berg G. The Identifiability of the Mixed Proportional Hazards Competing Risks Model. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2003;65:701–710.
  • American Cancer Society. Prostate Cancer: Detailed Guide. 2007
  • Asirvatham A, Schmidt M, Gao B, Chaudhary J. Androgens Regulate the Immune/Inammatory Response and Cell Survival Pathways in Rat Ventral Prostate Epithelial Cells. Endocrinology. 2006;147:257–271. [PubMed]
  • Bennett S. Analysis of Survival Data by the Proportional Odds Model. Statistics in Medicine. 1983;2:273–277. [PubMed]
  • Berkson J, Gage R. Survival Curve for Cancer Patients Following Treatment. Journal of the American Statistical Association. 1952;47:501–515.
  • Bravi F, Scotti L, Bosetti C, Talamini R, Negri E, Montella M, Franceschi S, La Vecchia C. Self-reported History of Hypercholesterolaemia and Gallstones and the Risk of Prostate Cancer. Annals of Oncology. 2006;17:1014–1017. [PubMed]
  • Chen M, Ibrahim J, Sinha D. A New Bayesian Model for Survival Data with a Surviving Fraction. Journal of the American Statistical Association. 1999;94:909–919.
  • Cheng S, Wie L, Ying Z. Analysis of Transformation Models with Censored Data. Biometrika. 1995;82:835–845.
  • Cox D. Regression Models and Life-tables (with Discussion) Journal of the Royal Statistical Society, Series B (Methodological) 1972;34:187–220.
  • Farewell V. The Use of Mixture Models for the Analysis of Survival Data with Long-term Survivors. Biometrics. 1982;38:1041–1046. [PubMed]
  • Fine J. Analysing Competing Risks Data with Transformation Models. Journal of the Royal Statistical Society. Series B (Statistical Methodology) 1999;61:817–830.
  • Fine J, Gray R. A Proportional Hazards Model for the Subdistribution of a Competing Risk. Journal of the American Statistical Association. 1999;94:496–509.
  • Fine J, Ying Z, Wei L. On the Linear Transformation Model for Censored Data. Biometrika. 1998;85:980–986.
  • Fleming T, Harrington D. Counting Processes and Survival Analysis. Wiley; 1991.
  • Fleshner N, Bagnell P, Kltoz L, Venkateswaran V. Dietary Fat and Prostate Cancer. The Journal of Urology. 2004;171:19–24. [PubMed]
  • Foutz R. On the Unique Consistent Solution to the Likelihood Equations. Journal of the American Statistical Association. 1977;72:147–148.
  • Hanlon A, Hanks G. Failure Patterns and Hazard Rates for Failure Suggest the Cure of Prostate Cancer by External Beam Radiation. Urology. 2000;55:725–729. [PubMed]
  • Heckman J, Honore B. The Identifiability of the Competing Risks Model. Biometrika. 1989;76:325–330.
  • Homma Y, Kondo Y, Kaneko M, Kitamura T, Nyou W, Yanagisawa M, Yamamoto Y, Kakizoe T. Promotion of Carcinogenesis and Oxidative Stress by Dietary Cholesterol in Rat Prostate. Carcinogenesis. 2004;25:1011–1014. [PubMed]
  • Hughes-Fulford M, Tjandrawinata R, Li C, Sayyah S. Arachidonic Acid, an Omega-6 Fatty Acid, Induces Cytoplasmic Phospholipase A2 in Prostate Carcinoma Cells. Carcinogenesis. 2005;26:1520–1526. [PubMed]
  • Jani A, Hellman S. Early Prostate Cancer: Clinical Decision-Making. The Lancet. 2003;361:1045–1053. [PubMed]
  • Keating N, O’Malley A, Smith M. Diabetes and Cardiovascular Disease During Androgen Deprivation Therapy for Prostate Cancer. Journal of Clinical Oncology. 2006;24:4448–4455. [PubMed]
  • Klebanov L, Yakovlev A. A New Approach to Testing for Sufficient Follow-up in Cure-rate Analysis. Journal of Statistical Planning and Inference. 2007;137:3557–3569.
  • Kosorok M, Song R. Inference Under Right Censoring for Transformation Models with a Change-point Based on a Covariate Threshold. The Annals of Statistics. 2007;35:957–989.
  • Kuk A, Chen C. A Mixture Model Combining Logistic Regression with Proportional Hazards Regression. Biometrika. 1992;79:531–541.
  • Lee S. Identification of a Competing Risks Model with Unknown Transformations of Latent Failure Times. Biometrika. 2006;93:996.
  • Li C, Taylor J, Sy J. Identifiability of Cure Models. Statistics & Probability Letters. 2001;54:389–395.
  • Li Y, Tiwari R, Guha S. Mixture Cure Survival Models with Dependent Censoring. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2007;69:285–306.
  • Lin D, Wei L, Ying Z. Checking the Cox Model with Cumulative Sums of Martingale-based Residuals. Biometrika. 1993;80:557–572.
  • Lu W, Ying Z. On Semiparametric Transformation Cure Models. Biometrika. 2004;91:331–343.
  • Maller R, Zhou S. Estimating the Proportion of Immunes in a Censored Sample. Biometrika. 1992;79:731–739.
  • Maller R, Zhou S. Testing for Sufficient Follow-Up and Outliers in Survival Data. Journal of the American Statistical Association. 1994;89:1499–1506.
  • Maller R, Zhou S. Testing for the Presence of Immune or Cured Individuals in Censored Survival Data. Biometrics. 1995;51:1197–1205. [PubMed]
  • Miloslavsky M, Keles S, Laan M, Butler S. Recurrent Events Analysis in the Presence of Time-dependent Covariates and Dependent Censoring. Journal of the Royal Statistical Society Series B. 2004;66:239–257.
  • Pandini G, Mineo R, Frasca F, Roberts C, Marcelli M, Vigneri R, Belfiore A. Androgens Up-regulate the Insulin-like Growth Factor-I Receptor in Prostate Cancer Cells. Cancer Research. 2005;65:1849–1857. [PubMed]
  • Parner E. Asymptotic Theory for the Correlated Gamma-frailty Model. Annals of Statistics. 1998;26:183–214.
  • Peng L, Huang Y. Survival Analysis with Quantile Regression Models. Journal of the American Statistical Association. 2008;103:637–649.
  • Peng Y, Dear K. A Nonparametric Mixture Model for Cure Rate Estimation. Biometrics. 2000;56:237–243. [PubMed]
  • Pepe M, Mori M. Kaplan-Meier, Marginal or Conditional Probability Curves in Summarizing Competing Risks Failure Time Data? Statistics in Medicine. 1993;12:737–751. [PubMed]
  • Prentice R, Kalbeisch J, Peterson A, Jr, Flournoy N, Farewell V, Breslow N. The Analysis of Failure Times in the Presence of Competing Risks. Biometrics. 1978;34:541–554. [PubMed]
  • Ritch C, Wan R, Stephens L, Taxy J, Huo D, Gong E, Zagaja G, Brendler C. Dietary Fatty Acids Correlate With Prostate Cancer Biopsy Grade and Volume in Jamaican Men. The Journal of Urology. 2007;177:97–101. [PubMed]
  • Robins J. Estimation of the Time-dependent Accelerated Failure Time Model in the Presence of Confounding Factors. Biometrika. 1992;79:321–334.
  • Robins J. Information Recovery and Bias Adjustment in Proportional Hazards Regression Analysis of Randomized Trials Using Surrogate Markers. Proceedings of the Biopharmaceutical Section, American Statistical Association. 1993:24–33.
  • Robins J, Finkelstein D. Correcting for Noncompliance and Dependent Censoring in an AIDS Clinical Trial with Inverse Probability of Censoring Weighted (IPCW) Log-Rank Tests. Biometrics. 2000;56:779–788. [PubMed]
  • Robins J, Rotnitzky A. Recovery of Information and Adjustment for Dependent Censoring Using Surrogate Markers. In: Jewell N, Dietz K, Farewell V, editors. AIDS Epidemiology: Methodological Issues. Birkhauser; 1992. pp. 297–331.
  • Rotnitzky A, Robins J. Semiparametric Regression Estimation in the Presence of Dependent Censoring. Biometrika. 1995;82:805–820.
  • Rotnitzky A, Robins J. Inverse Probability Weighted Estimation in Survival Analysis. In: Armitage P, Colton T, editors. The Encyclopedia of Biostatistics, Second Edition. Wiley; 2003. pp. 2619–2625.
  • Rudin W. Functional Analysis. McGraw-Hill; 1987.
  • Stapleton JH. Linear Statistical Models. Wiley; 1995.
  • Strom S, Yamamura Y, Forman M, Pettaway C, Barrera S, Digiovanni J. Saturated Fat Intake Predicts Biochemical Failure After Prostatectomy. International Journal of Cancer. 2008;122:2581–2585. [PubMed]
  • Sy J, Taylor J. Estimation in a Cox Proportional Hazards Cure Model. Biometrics. 2000;56:227–236. [PubMed]
  • Tai P, Yu E, Cserni G, Vlastos G, Royce M, Kunkler I, Vinh-Hung V. Minimum Follow-up Time Required for the Estimation of Statistical Cure of Cancer Patients: Verification Using Data from 42 Cancer Sites in the SEER Database. BMC cancer. 2005;5:48. [PMC free article] [PubMed]
  • Tsodikov A, Ibrahim J, Yakovlev A. Estimating Cure Rates from Survival Data: An Alternative to Two-component Mixture Models. Journal of the American Statistical Association. 2003;98:1063–1079. [PMC free article] [PubMed]
  • van der Vaart A. Asymptotic Statistics. Cambridge University Press; 1998.
  • van der Vaart A, Wellner J. Weak Convergence and Empirical Processes. Springer; 1996.
  • Wellner J, Zhan Y. Bootstrapping Z-estimators. University of Washington Department of Statistics Technical Report. 1996:308.
  • Wilson P, D’Agostino R, Levy D, Belanger A, Silbershatz H, Kannel W. Prediction of Coronary Heart Disease Using Risk Factor Categories. Circulation. 1998;97:1837–1847. [PubMed]
  • National Cancer Institute, DCCPS, Surveillance Research Program, Cancer Statistics Branch, Surveillance, Epidemiology, and End Results (SEER) Program, Public-Use Data (1973–2001) 2004. Released April 2004, based on November 2003 submission www.seer.cancer.gov.
  • Xu Y, Chen S, Ross K, Balk S. Androgens Induce Prostate Cancer Cell Proliferation through Mammalian Target of Rapamycin Activation and Post-transcriptional Increases in Cyclin D Proteins. Cancer Research. 2006;66:7783–7792. [PubMed]
  • Zeng D, Lin D. Maximum Likelihood Estimation in Semiparametric Regression Models with Censored Data. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2007;69:507–564.
  • Zhuang L, Kim J, Adam R, Solomon K, Freeman M. Cholesterol Targeting Alters Lipid Raft Composition and Cell Survival in Prostate Cancer Cells and Xenografts. Journal of Clinical Investigation. 2005;115:959–968. [PubMed]