I've been trying for multiple days to visualize my Cox regression with time-dependent covariates with no luck. I'd like to isolate one of two time-dependent covariates to visualize its effect on survival probability, risk, hazard ratio, and/or expected failures while holding the other variable constant at its mean.
As an example, I'm using data from the pbc
dataset from the survival
package in R. The data is set up as follows using the tmerge
function:
data("pbc", package = "survival")
head(pbc)[, c(1:5, 11, 12)]
head(pbcseq)[, c(1, 4:5, 7, 12, 13)]
pbc <- pbc %>% mutate(bili = log(bili), protime = log(protime)) # log = e^x = y; y = bili
pbcseq <- pbcseq %>% mutate(bili = log(bili), protime = log(protime))
temp <- subset(pbc, id <= 312, select = c(id:sex)) # baseline
pbc2 <- tmerge(temp, temp, id = id, death = event(time, status)) #set range
pbc2 <- tmerge(pbc2, pbcseq, id = id, bili = tdc(day, bili),
protime = tdc(day, protime))
fit <- coxph(Surv(tstart, tstop, death == 2) ~ bili + protime, pbc2)
What I'd like to do is make a new dataframe to predict over using the output from the model (fit
). I'd like to hold the protime
variable constant at its mean while predicting the effect of bili
on survival probability, the hazard rate, and/or the number of failures expected. I'd then like to plot these effects.
I've tried using a survfit
object to plot survival. I've tried using the expand.grid
function to create a new dataset, but because the new data is not formatted using time intervals I get an error. I've tried the make_newdata
function (from pammtools
) to create a new dataset and then the predict()
function to apply the model over the new dataset, but it will not take a survfit
object and I'm not sure how to get my values of interest (survival, hazard, etc) by just using the fit
model.
I would be SO grateful if someone can help me!!