Adjusted survival curves for cox model with time dependent covariate 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.


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

 A: 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. 
