4
$\begingroup$

Is it appropriate to use cox regression to study a process that doesn't result in the survival function going to 0 when time = infinity? Cox regression is interested in hazard ratios, and it makes no attempt to actually estimate/model the hazard function itself, so intuitively it makes sense to me that this should be OK. Am I missing something? (If it helps, I can find a lot of examples in the literature of people doing this, but that doesn't mean it's right...)

As an example, suppose I want to study divorce rates... In this case, plenty of people never get divorced, so marriage survival will never go to 0.

$\endgroup$
  • $\begingroup$ A search term regarding your question is "cure rate". $\endgroup$ – Michael M Sep 6 '17 at 20:28
  • $\begingroup$ See also frailty models, discussed e.g. in the works by Aalen cited in this answer. Some marriages are 'frailer' than others, and clearly it would be of interest to investigate the factors contributing to frailty vs resilience in this context. $\endgroup$ – David C. Norris Sep 9 '17 at 11:39
4
+50
$\begingroup$

There is quite some literature on survival analysis under population heterogeneity, but like you notice yourself, I rarely see such models being used - or even considered - in applied research. I'll give some brief intuition, and hopefully others can add mathematically-heavier explanations if that's what you're looking for.

So let's say our observations look like this - quarter of the sample is "cured" and does not experience the event during observation period, while the rest are "susceptible" and slowly dying/failing/divorcing:

plot

The entire blue fraction will be censored at $t=50$, and you can have additional censoring in both groups as usual (not seen in the plot).

Now, Cox PH models assume that survival time is independent of censoring. However, if you fit such a model to the pictured sample, the censored group will have much higher time-to-event a different hazard function because it contains the entire "cured" fraction. Simulations (e.g. 2) show that the HR estimate can be strongly biased under these conditions.
The problem could be solved by excluding the "cured" observations, but those usually are not distinguishable from censored, but "susceptible" individuals. Another analysis option is to only consider the binary outcome event vs. no event, but this again means excluding all censored individuals.

A conceptually simple solution is to treat the sample as a mixture of "cured" and "susceptible" distributions with weights $p$ and $1-p$, and fit a mixture model with survival $S(t) = 1$ for the "cured" fraction and some decreasing $S(t)$ for the "susceptible" fraction. These are known as mixture cure models; there are also non-mixture solutions, but I'm unsure if those are popular in practice.

Some references, in the order from most-layman to most mathematical:

Cure Models as a Useful Statistical Tool for Analyzing Survival
Mixture and non-mixture cure fraction models based on the generalized modified Weibull distribution with an application to gastric cancer data
What Cure Models Can Teach us About Genome-Wide Survival Analysis
A General Approach for Cure Models in Survival Analysis


EDIT

Given the good comments by @JarleTufto below, I should clarify this answer. I do not propose that merely censoring at the end of observation period is bad, or that the hazard after final $t$ is somehow relevant, but that problems are caused by underlying heterogeneity in population - i.e. fractions with different hazard functions $h(t)$ are censored differently. The corresponding Cox model assumption is stated as: $$h(t|covariates) = h(t|C_i > t, covariates)$$

(e.g. Relaxing the independent censoring assumption in the Cox proportional hazards model using multiple imputation)

Let's take the simulation code from the answer below, and add a "cured" fraction with $h(t)=0$ (also added seed and increased sample size for easier reproducibility):

set.seed(1234)

# simulated survival times from the model
n <- 5000
x <- rnorm(n)
beta <- 0.5
# variation of inversion method
u <- runif(n)
eventtime <- 1/(1 + exp(-beta*x)*log(u)) - 1 
eventtime[eventtime < 0] <- Inf

# simulate independent right censoring points
censoringtime <- runif(n,0,20)

# compute the observed data, that is, the censoring indicator and
# the time of whichever event comes first
delta <- eventtime < censoringtime
time <- pmin(eventtime, censoringtime)

# add "cured" fraction, censored at max observation time:
n2 <- 1000
time2 <- c(time, rep(20, n2))
delta2 <- c(delta, rep(FALSE, n2))
x2 <- c(x, rnorm(n2))

So the observed (event or censoring) timepoints look like:

qplot(time2)

plot2

Using only the "susceptible" fraction, we get the expected estimate of $\beta=0.5$:

# Fit the cox proportional hazards model
library(survival)
model <- coxph(Surv(time,delta) ~ x)
model

    coef exp(coef) se(coef)    z      p
x 0.5067    1.6598   0.0198 25.6 <2e-16

Likelihood ratio test=661  on 1 df, p=0
n= 5000, number of events= 2923 

However, using the full sample, the HR is underestimated:

model2 <- coxph(Surv(time2,delta2) ~ x2)
model2

    coef exp(coef) se(coef)    z      p
x2 0.4192    1.5207   0.0192 21.8 <2e-16

Likelihood ratio test=478  on 1 df, p=0
n= 6000, number of events= 2923 
$\endgroup$
  • $\begingroup$ Accepted for "Cox PH models assume that survival time is independent of censoring." I knew Cox models had that assumption, but I hadn't realized it would be violated in this case. Thanks! $\endgroup$ – John Chrysostom Sep 8 '17 at 11:31
  • $\begingroup$ -1 @John This answer is wrong. Independence between survival times $T$ and right censoring times $C$ means that $P(T\le t\cap C\le c)=P(T\le t)P(C\le c)$. Right censoring happens when $T>C$ so there will always be a tendency that those who are right censored have longer lifetimes $T$. $\endgroup$ – Jarle Tufto Sep 8 '17 at 13:19
  • 1
    $\begingroup$ @John With random independent censoring, the censoring times $C$ are random (sometimes unobserved) variables that are independent of the survival times $T$. In general, even when $T$ and $C$ are independent random variables, the distribution of $T$ conditional on the event $T>C$ (the survival time $T$ becoming right censored) is different from the unconditional distribution of $T$. So even when independence holds it follows that the "censored group will have a much higher time-to-event" as observed by juod but this do not imply that the assumption of independence is violated. $\endgroup$ – Jarle Tufto Sep 8 '17 at 20:08
  • 1
    $\begingroup$ @JarleTufto you are certainly right, and my phrasing was indeed sloppy. But my conclusion still holds for the situation I described and I believe is relevant to OP's example. I rephrased it - but of course, my semi-applied sources might be misrepresenting the basic idea as well. Please see if you agree with the edit. $\endgroup$ – juod Sep 9 '17 at 15:06
  • 1
    $\begingroup$ Yes, if $S(t;x)$ goes to the same limiting value for different values of the covariates $x$, for example because some part of the population has acquired immunity, this would be inconsistent with the model (the limiting value of $S(t;x)$ needs to depend on $x$ under the Cox model) such that the proportional hazard assumption would no longer hold, something which would show up in plots of the Schoenfeld residuals. But that "survival doesn't go to zero" do not imply that this will be the case. And I'm assuming independent censoring here. Dependent censoring would further complicate matters. $\endgroup$ – Jarle Tufto Sep 9 '17 at 18:13
2
$\begingroup$

It may perfectly fine to apply the Cox model in the situation you describe. The Cox model makes no assumptions about the baseline hazard $h_0(t)$ other than that it is non-negative for all $t$. So the baseline hazard may well go to zero fast enough to make the cumulative hazard $\int_0^t h_0(u) du$ go to a finite value and $S(t)$ to some positive value as $t\rightarrow \infty$.

Unlike parametric survival models, no inference is made about the hazard function beyond the last censoring event/last observed failure. Instead, just as the simpler Kaplan-Meier estimator of $S(t)$, the Cox model gives you a non-parametric estimate of the baseline survival function $S_0(t)$ only up to the last censoring event/last observed failure. So the theoretical behaviour of the hazard function after this time point is of no consequence whatsoever.

For example, suppose that the true unknown baseline hazard has the form $$ h_0(t) = \frac1{(1+t)^2}. $$ The survival function for an individual with covariate vector $x$ is then \begin{align} S(t;x) &= e^{-e^{\beta x}\int_0^t h_0(u)du} \\&= e^{-e^{\beta x}(1-\frac1{t+1})}. \end{align} The cumulative baseline hazard $1-\frac1{t+1}$ then tends to a limiting value (of one) and the survival function $S(t;x)$ to the positive probability $e^{-e^{\beta x}}$ as $t\rightarrow\infty$. Note that the limiting value of $S(t;x)$ in general must depend on the covariates $x$, if this is not the case the proportional hazard assumption will be violated.

In the following, survival times $T$ from this model (some of which are infinite) and independent right censoring times $C$ uniformly distributed on the interval from 0 to 20 are simulated and the Cox proportional hazard model is fitted to the observed data (the censoring indicator and the minimum of $T$ and $C$).

# simulated survival times from the model
n <- 300
x <- rnorm(n)
beta <- .5
# variation of inversion method
u <- runif(n)
eventtime <- 1/(1 + exp(-beta*x)*log(u)) - 1 
eventtime[eventtime < 0] <- Inf

# simulate independent rigth censoring points
censoringtime <- runif(n,0,20)

# compute the observed data, that is, the censoring indicator and
# the time of whichever event comes first
delta <- eventtime < censoringtime
time <- min(eventtime, censoringtime)

# Fit the cox proportional hazards model
library(survival)
model <- coxph(Surv(time,delta) ~ x)

Both the regression coefficient $\beta=0.5$,

> model
Call:
coxph(formula = Surv(time, delta) ~ x)

    coef exp(coef) se(coef)    z       p
x 0.5463    1.7269   0.0892 6.12 9.2e-10

Likelihood ratio test=39.6  on 1 df, p=3.19e-10
n= 300, number of events= 164     Call:

and the baseline survival function $S_0(t)$ are estimated just fine.

# Compare the estimated and true baseline survival functions
plot(survfit(model, newdate=data.frame(x=0)))
curve(exp(-(1-1/(t+1))),xname="t",add=TRUE,col="red")

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.