2
$\begingroup$

I have found two ways to get an adjusted survival curve from a cox model with time dependent covariates (A and B). I am wondering which one (or if I am completely wrong how?) would be the appropriate way to graph survival of individuals stratified on a categorical variable that changes over the course of observation.

Approach A Approach B

Since the graphs do look inherently different I am wondering which one is correct or if they are correct at all. Also would an interpretation like

individuals with score 4 had the lowest probability of survival with a distinct curve deviation from lower values

be valid?

Everything should be in accordance with https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf

library(survival)
library(survminer)
library(tidyverse)

set.seed(199)

df1 <- data.frame(ID = rep(seq(1, 100, by = 1), 2),
                  score = factor(sample(1:4, 200, replace = TRUE)),
                  timetoFU = sample(300:3000, 200, replace = TRUE),
                  status = 0
                  )
df1 <- df1 %>% group_by(ID) %>% arrange(ID, timetoFU) %>% mutate(obs_n = row_number(), time_max = last(timetoFU)+sample(200:400, 1)) %>% ungroup()


temp <- df1 %>% filter(obs_n == 1) %>% 
  mutate(status = replace(status, score == 1, sample(c(0,1), 1, prob = c(0.6, 0.4)))) %>% 
  mutate(status = replace(status, score == 2, sample(c(0,1), 1, prob = c(0.4, 0.6)))) %>% 
  mutate(status = replace(status, score == 3, sample(c(0,1), 1, prob = c(0.3, 0.7)))) %>% 
  mutate(status = replace(status, score == 4, sample(c(0,1), 1, prob = c(0.2, 0.8))))

td_df <- tmerge(temp, temp, id = ID, outcome = event(time_max, status))
td_df <- tmerge(td_df, df1, id = ID, td_score = tdc(timetoFU, score))

#A
cox_fit_A <- coxph(formula = Surv(tstart, tstop, status) ~ td_score, data = td_df, id = ID)
summary(cox_fit_A)

adjusted_surv_A <- ggadjustedcurves(cox_fit_A, data = td_df, variable = "td_score") 

adjusted_surv_A

#B
cox_fit_B <- coxph(formula = Surv(tstart, tstop, status) ~ td_score + strata(td_score), data = td_df, id = ID)
summary(cox_fit_B)

adjusted_surv_B <- ggadjustedcurves(cox_fit_B, data = td_df) 

adjusted_surv_B
$\endgroup$
3
$\begingroup$

The difference between these two plots is whether or not you constrain the model to satisfy proportional hazards -- the difference is in the Cox model, not in the adjusted curve

The first model has four curves for the four groups. The hazard for individuals in different groups is proportional; you can see this because the curves all have steps at the same time values.

The second model also has four curves for the four groups, but they are supplied as strata not as predictors, so there is no constraint at all -- the fitted curve for an individual is simply the Kaplan-Meier survival curve for that group.

These are both graphs showing estimated survival for someone whose td_score does not change over time, under the assumption that td_score affects survival and not the other way around.

It's not clear in general whether constraining the estimates to satisfy proportional hazards is good or bad. I think your code that generates the data does satisfy proportional hazards, which would argue for the first graph, but in general you won't know that. The only general rule is that the more data you have, the better the case for dropping the proportional hazards constraint: the constraint will decrease variance and may increase bias, and in larger samples you care more about the bias.

Things would be more complicated if there was another variable as well as td_score: the different adjustment methods provided by ggadjustedcurves would average over these other variables in different ways.

$\endgroup$
  • $\begingroup$ Thank you very much for the clarification, I think I understand now why the two curves look different. However you mentioned that both graphs show only estimated survival for individuals whos score doesn't change over time. Does that mean others are deleted for the calculation of the survival curves. So far I thought individuals would jump curves if their score changes. And could you also clarify what would be exactly different in the calculation of adjusted survival curves if there is another variable? $\endgroup$ – Rootless17b May 29 '20 at 21:45
  • 1
    $\begingroup$ No, all the data are used, but the predictions are for people who don't change their score over time. You could also get predictions for people who do change over time in various ways -- the survival package vignette describes how to do this -- but there are infinitely many ways to change and only four ways to stay the same. $\endgroup$ – Thomas Lumley May 30 '20 at 6:33
  • $\begingroup$ Thanks Thomas!! $\endgroup$ – Rootless17b May 31 '20 at 10:04

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.