Cox Regression when survival doesn't go to 0?

Question

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.

Answer 1

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:

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)

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 

Answer 2

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 <- pmin(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")