Survival Analysis – Are Cumulative Hazard Rates More Accurate Than Survival Estimates in Cox Models?

cox-modelsurvivaltime-varying-covariate

I am using a time-varying covariate to represent subjects that can enter and leave groups at arbitrary times. Because group membership is not fixed, I posit that the true time-to-event is more accurately represented by cumulative hazard rates than by survival estimates.

The calculation of hazard rates has the interpretation of the number of events per unit time among those at risk, while survival estimates imply some probability of an event among a fixed cohort of subjects moving forward in time. Intuitively, I cannot describe what a survival curve is measuring when group membership is changing; the probability of an individual in a group surviving to the next time point given prior survival to the current time point seems unsupportable.

Do cumulative hazard rates have an advantage for these reasons?

If not, what is the interpretation of a survival curve in the presence of a time-dependent covariate?

Reproducible example:

library(dplyr)
library(survminer)
library(survival)

set.seed(42)

id <- rep(c(1:10), 100)
n <- length(id)
event <- rbinom(n, 1, 0.5)
group <- sample(c(1,2),n, TRUE)
time <- abs(rnorm(n,0,1)* 10 * group)
df <- data.frame(id, event,group, time)

df <- df %>% arrange(id) %>% group_by(id) %>%
      mutate(t0 = lag(time, default = 0, order_by = id),
             t1 = lead(t0 + time, default = time, order_by = id))

surv.obj <- with(df, Surv(t0, t1, event, type = "counting") )
survfit.obj <- survfit(surv.obj ~ group, df)

ggsurvplot(survfit.obj, fun = "cumhaz", df, censor = F)
ggsurvplot(survfit.obj, fun = "pct", df, censor = F)

Best Answer

There isn't any conceptual advantage of cumulative hazard over instantaneous hazard rates. The cumulative hazard is just the integral of the instantaneous hazard over time. With respect to your conceptual issue,

I cannot describe what a survival curve is measuring when group membership is changing; the probability of an individual in a group surviving to the next time point given prior survival to the current time point seems unsupportable

Therneau and Grambsch on pages 271-2 suggest that you can

imagine some hypothetical "cohort" of subjects who follow a given path [of covariate values], losing members to death along the way

as a way to think about modeling time-varying covariates.

Part of the problem is that you use the word "group" in two different ways. In your case, the "group membership" is a time-varying covariate. But in the phrase "a fixed group of subjects moving forward in time," what you are talking about is an initial cohort of subjects (as in Therneau and Grambsch) whose "group membership" can change over time.

You do raise an important point in making predictions with time-varying covariates. It's very easy to hypothesize time courses for covariates that aren't realistic. Furthermore, making a prediction about a covariate value at a given time implies that an individual has survived up to that time. Survivorship bias is a serious risk.

Related Solutions

Cox Model – How to Analyze Cox Survival Data with Time-Varying Continuous Variables and Fixed Covariates in R

I figure you are looking for the model in equation 4 in this article. Then this code will do

> library(survival)
> set.seed(747)
> 
> N <- 1e4 # number of subjects
> k <- 0.3 # assuming the time varying covariates is proporitonal to time x2=k*t
> lambda <- 0.01
> betat <- 0.3 
> beta <- -0.6
> rateC <- 0.01
> 
> #####
> # simulate 
> x1 <- as.numeric(runif(N) > .5) # covariate
> 
> # simulate outcome (NB: OP forgot parentheses)
> u <- runif(n = N)
> Tlat <- log(1 + betat*k*(-log(u))/(lambda*exp(beta*x1))) / (betat * k)
> 
> # censoring and status indicator
> Cen<- rexp(n = N, rate = rateC)
> 
> # make data.frame
> df <- data.frame(
+   time = pmin(Tlat, Cen), 
+   status = Tlat <= Cen, 
+   id = 1:N, 
+   x1 = x1, 
+   x2 = rep(k, N))
> 
> # fit model -- cannot use coxph for x2 due to the non-parametric intercept
> fit <- coxph(Surv(time, status) ~ x1, data = df)
> fit # coefficient for x1 is right
Call:
coxph(formula = Surv(time, status) ~ x1, data = df)

      coef exp(coef) se(coef)   z      p
x1 -0.6083    0.5442   0.0234 -26 <2e-16

Likelihood ratio test=675  on 1 df, p=0
n= 10000, number of events= 7891 
> 
> # # won't work as we try to estimate a term that is equal for all individuals 
> # # at each time t
> # fit <- coxph(Surv(time, status) ~ x1 + tt(x2), data = df, 
> #              tt = function(x, t, ...) x * t)
> 
> # The non-parametric cumulative hazard though match with what we expect
> base <- basehaz(fit, centered = FALSE)
> plot(base$time, base$hazard, type = "l")
> lines(
+   base$time, 
+   lambda / (betat * k) * (exp(betat * k * base$time) - 1), col = "red")

> #####
> # estimate parametric model
> neg_loglike_func <- function(b){
+   # compute intermediates
+   lp1 <- exp(b[1] + b[2] * df$x1)    # time-invariant terms
+   lp2 <- exp(b[3] * df$x2 * df$time) # time-varying terms
+   log_haz <- log(lp1) + log(lp2)                  # instant hazard
+   cumhaz <- lp1 * (lp2 - 1)/(b[3] * df$x2 + 1e-8) # cumhaz
+   
+   # compute log likelihood terms
+   ll_terms <- ifelse(df$status == 1,  log_haz - cumhaz, -cumhaz)
+   -sum(ll_terms) # neg log likelihood
+ }
> 
> # fit model
> b <- c(0.01, -0.02, 0.01) # starting valus
> fit <- optim(b, neg_loglike_func)
> 
> exp(fit$par[1]) # lambda estimate
[1] 0.01035017
> fit$par[2:3]    # beta and betat estimate
[1] -0.6099467  0.2966224

You may be able to avoid the optim function above by e.g., using one of the parametric survival models in the rstpm2 package or polspline.

Before edits

I cannot see how what you do is related to the cited article. Further, I am struggling to see how x2 enters into the model. Maybe post equations for the model you are trying to simulate from?

What may help you is that all of your observation end up dying at the end it seems which seems odd given that you have # censoring times in the Cen object. I show this below

library(survival)
set.seed(747)

N=1000 # number of subjects
k=0.3 # assuming the time varying covariates is proporitonal to time x2=k*t
lambda=0.01
betat=0.3 
beta=-0.6
rateC=0.01


#####
# Time fixed variable
x1 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5))

#Exponential latent event times
u <- runif(n=N)
Tlat <- log(1+betat*k*(-log(u)) / (lambda*exp(beta*x1))) / betat*k
Cen<- rexp(n=N, rate=rateC) #censoring times

# follow-up times and event indicators
time <- pmin(Tlat, Cen)
status <- as.numeric(Tlat <= Cen)
# data set
df.tfixed<-data.frame(id=1:N, time=time, status=status, x1=x1)


#####
# time dependent continuous variable
ntp<-sample(1:6,N,replace=T) #number of follow up time points
mat<-matrix(ncol=3,nrow=1)
i=0
for(n in ntp){
  i=i+1
  ft<-runif(n,min=0,max=df.tfixed$time[i])
  ft<-sort(ft)
  seq<-rep(ft,each=2)
  seq<-c(0,seq,df.tfixed$time[i])
  matid<-cbind(matrix(seq,ncol=2,nrow=n+1,byrow=T),i)
  mat<-rbind(mat,matid)
}
df.td<-data.frame(mat[-1,])
colnames(df.td)<-c("start","stop","id")
df.td$x2<-k*df.td$start


#combine the two data frames 
df<-merge(df.td,df.tfixed,by="id")
df$status=0
df$status[cumsum(as.vector(table(df$id)))]<-1

fit <- coxph(Surv(start,stop, status) ~ x1+x2, data=df)
fit$coef
>         x1         x2 
> -0.6040089  1.8471387 

# All observations dies exaclty once -- no censoring
all(tapply(df$status == 1, df$id, sum) == 1)
> [1] TRUE

It comes down to this line df$status[cumsum(as.vector(table(df$id)))]<-1.

Best Answer

Related Solutions

Cox Model – How to Analyze Cox Survival Data with Time-Varying Continuous Variables and Fixed Covariates in R

Before edits

Related Question