Suppose there are two population subgroups: the frail (with greater risk of death) and the robust (with lower risk of death). The code below simulates data from this hypothetical scenario (click the button to see the code).
library(tidyverse)
library(survival)
simulated <- tibble(id = 1:1e3) |>
mutate(
# Make equal-sized subgroups that are frail and robust
subgroup = rep(c("frail","robust"), n() / 2),
# Define exponential rate parameter in each subgroup
rate = case_when(
subgroup == "frail" ~ 10,
subgroup == "robust" ~ 1
),
# Generate an Exponential outcome
t = rexp(n(), rate = rate),
# For simplicity, assume no censoring
event = 1
)Below is a visualization of the time until death in each subgroup of this simulation. The frail subgroup dies in a short time, while the robust subgroup sometimes lives much longer.
Now imagine you see these data, but you do not know which people are frail and which are robust. You just see the mixture distribution
Without knowing about the frailty issue, you might estimate an Exponential survival model with no predictors.
expo_model <- survreg(
Surv(t, event) ~ 1,
data = simulated,
dist = "exponential"
)You might then visualize the estimated hazard function and survival function. I wrote the viz_survreg() function to make it easier to visualize these curves.
source("https://ilundberg.github.io/eventhistory/viz_survreg.R")expo_model |> viz_survreg()
To know if this is a good model, you might compare the survival curve to a Kaplan-Meier estimate. I wrote the km_compare() function to make it easier to visualize these curves.
source("https://ilundberg.github.io/eventhistory/km_compare.R")expo_model |> km_compare()
The comparison suggests something is wrong: the model overestimates survival at early time points, then underestimates survival at late time points. If you did not know the data generating process, you might wonder if this could arise because of a non-constant hazard function.
Knowing that a Weibull can model a hazard that changes with time, you might then proceed to a Weibull model.
weibull_model <- survreg(
Surv(t, event) ~ 1,
data = simulated,
dist = "weibull"
)Predictions from this model show a hazard that starts high and decreases with time,
weibull_model |> viz_survreg()
and in fact you can reject the null of the Exponential model (see summary(weibull_model) and note that the p-value on Log(scale) = 0 is very close to 0). You might also be encouraged by seeing that the predicted survival function closely matches Kaplan-Meier.
weibull_model |> km_compare()
The illustration above seems paradoxical.
What has happened here is a consequence unmeasured heterogeneity. Some of the people are robust, and others are frail. Because the frail die first, the proportion of survivors who are robust rises with time.
The consequence of this shifting population composition is the appearance of a decreasing hazard function. When frail and robust are unmeasured, a survival model fitted to the full population will estimate one hazard that decreases over time.
How much should we worry about unmeasured heterogeneity? On one hand, all properties of the model rely on the assumption that outcomes are actually generated by the model that is selected. On the other hand, a model that only approximates the truth may still yield good population-average or subgroup survival curves even if it does not correspond to the hazard of any individual person.
If you knew the form of the unmodeled heterogeneity, you could define a model that incorporates this unmodeled heterogeneity. In certain circumstances, a correctly-specified model can overcome latent heterogeneity.
For example, suppose we knew that simulated is actually a mixture of two exponential distributions. We could define a data generating process with this mixture.
\[ \begin{aligned} U &\sim\text{Bernoulli}(\text{probability} = \pi) \qquad &&\text{unmeasured frailty} \\ T &\sim \text{Exponential}(\text{rate} = \lambda_1^U\lambda_2^{1-U}) &&\text{survival time} \end{aligned} \]
where the survival time is generated with rate \(\lambda_1\) when \(U = 1\) and \(\lambda_2\) when \(U = 0\).
We can write down the probability density function of the observed data \(t_1,\dots,t_n\).
\[ \begin{aligned} f(t) &= f(T = t\mid U = 1)\text{P}(U = 1) + f(T = t\mid U = 0)\text{P}(U = 0) \\ &= \underbrace{f(t\mid\lambda = \lambda_1)}_\text{Exponential PDF}\pi + \underbrace{f(t\mid\lambda = \lambda_2)}_\text{Exponential PDF}(1 - \pi) \\ \end{aligned} \]
We can write this log likelihood function in code where parameters will be \(\log(\lambda_1)\), \(\log(\lambda_2)\), and \(\text{logit}(\pi)\).
log_likelihood <- function(parameters, data) {
lambda_1 <- exp(parameters[1])
lambda_2 <- exp(parameters[2])
pi <- plogis(parameters[3])
log_likelihood <- data |>
mutate(
likelihood_i = dexp(t, rate = lambda_1) * pi +
dexp(t, rate = lambda_2) * (1 - pi)
) |>
summarize(log_likelihood = sum(log(likelihood_i))) |>
pull(log_likelihood)
return(log_likelihood)
}We can optimize this log likelihood using optim.
optim.out <- optim(
par = c(0,0,0),
fn = log_likelihood,
control = list(fnscale = -1),
data = simulated
)This can successfully recover our data generating process! Below are the estimated parameters.
estimated_parameters <- tibble(
lambda_1 = exp(optim.out$par[1]),
lambda_2 = exp(optim.out$par[2]),
pi = plogis(optim.out$par[3])
)# A tibble: 1 × 3
lambda_1 lambda_2 pi
<dbl> <dbl> <dbl>
1 1.06 10.2 0.509We can plot a survival curve with these estimated parameters.
curves <- tibble(t = seq(.01,3,.01)) |>
mutate(
A = pexp(t, rate = estimated_parameters$lambda_1, lower.tail = FALSE),
B = pexp(t, rate = estimated_parameters$lambda_2, lower.tail = FALSE),
Pooled = A * estimated_parameters$pi + B * (1 - estimated_parameters$pi)
) |>
select(t, A, B, Pooled)
curves |>
pivot_longer(cols = -t, names_to = "Latent Subgroup", values_to = "Survival") |>
ggplot(aes(x = t, y = Survival, color = `Latent Subgroup`, linetype = `Latent Subgroup`)) +
geom_line() +
labs(x = "Time")
We can also see that the pooled survival curve is incredibly close to the Kaplan-Meier estimates.
survfit(Surv(t, event) ~ 1, data = simulated) |>
broom::tidy() |>
filter(time <= 3) |>
select(time, estimate) |>
mutate(method = "Kaplan-Meier") |>
bind_rows(
curves |>
rename(time = t, estimate = Pooled) |>
select(time, estimate) |>
mutate(method = "Latent Mixture of\nTwo Exponentials")
) |>
ggplot(aes(x = time, y = estimate, color = method, linetype = method)) +
geom_line() +
labs(x = "Time", y = "Survival", color = "Method", linetype = "Method")
In practice, it is uncommon to know the parametric form of latent heterogeneity. Thus, models that incorporate latent heterogeneity often rely on assumptions (such as a mixture of two Exponentials) that may themselves be imperfect.
The more important thing for practice is to recognize that latent heterogeneity is likely to exist in any survival model. When a hazard declines over time, it is important to recognize that this population-level result may arise from two equivalent sources
Data can be equally consistent with both interpretations. Nonetheless, one can still produce good survival curve estimates for the full population or in subgroups without distinguishing the two interpretations above.