sahirbhatnagar/casebase

plot hazard ratio as a function of time

sahirbhatnagar opened this issue · 11 comments

sahirbhatnagar commented

Since one of the main motivations for casebase is flexbile modeling of time, we should provide some functionality to plot time dependent hazard ratios. Some useful links:

sahirbhatnagar commented

I've tried working with visreg and ggeffects to get hazard ratio (HR) plots from the output of casebase::fitSmoothHazard. The code below compares these two packages along with using the output from the casebase::hazardPlot function. So far, visreg is consistent with casebase::hazardPlot, and ggeffects is slightly off, but im not sure why.

@Jesse-Islam @turgeonmaxime Any thoughts or comments on what to try next? Should we try going forward with visreg ? The advantage I see is that it nicely handles models with more than one or two predictors, i.e., you can set the levels of the other predictors in the model with various sensible defaults. This feature isn't implemented in casebase::hazardPlot which currently only handles one covariate pattern as a function of time.

# load packages and data ---------------------------------------------------

devtools::load_all("~/git_repositories/casebase/")
#> Loading casebase
#> See example usage at http://sahirbhatnagar.com/casebase/
pacman::p_load(splines)
pacman::p_load(ggeffects)
pacman::p_load(visreg)
data("simdat")
head(simdat)
#>   id eventtime status trt
#> 1  1 3.3804880      1   0
#> 2  2 2.8934869      1   0
#> 3  3 2.7613184      1   0
#> 4  4 5.0000000      0   1
#> 5  5 2.6698927      1   0
#> 6  6 0.7441971      1   0

# Fit the model on simdat from casebase -----------------------------------

mod_cb <- casebase::fitSmoothHazard(status ~ bs(eventtime,degree = 8) * trt,
                                    time = "eventtime",
                                    event = "status",
                                    data = simdat,
                                    ratio = 10)
summary(mod_cb)
#> 
#> Call:
#> glm(formula = formula, family = binomial, data = sampleData)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -0.6447  -0.5079  -0.4083  -0.3344   3.1492  
#> 
#> Coefficients:
#>                                Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)                     -3.8344     0.6348  -6.041 1.53e-09 ***
#> bs(eventtime, degree = 8)1       5.9234     3.0051   1.971   0.0487 *  
#> bs(eventtime, degree = 8)2      -4.5849     5.8681  -0.781   0.4346    
#> bs(eventtime, degree = 8)3       8.4348    11.2906   0.747   0.4550    
#> bs(eventtime, degree = 8)4       0.8944    12.5621   0.071   0.9432    
#> bs(eventtime, degree = 8)5      -0.6884    10.6707  -0.065   0.9486    
#> bs(eventtime, degree = 8)6       6.4765     5.4857   1.181   0.2378    
#> bs(eventtime, degree = 8)7       1.0430     2.2042   0.473   0.6361    
#> bs(eventtime, degree = 8)8       2.8906     0.7062   4.093 4.26e-05 ***
#> trt                             -2.1283     1.5317  -1.390   0.1647    
#> bs(eventtime, degree = 8)1:trt   1.7757     6.3580   0.279   0.7800    
#> bs(eventtime, degree = 8)2:trt   4.8295    10.3108   0.468   0.6395    
#> bs(eventtime, degree = 8)3:trt  -3.7247    19.2117  -0.194   0.8463    
#> bs(eventtime, degree = 8)4:trt   4.7817    19.5752   0.244   0.8070    
#> bs(eventtime, degree = 8)5:trt   1.5612    16.4359   0.095   0.9243    
#> bs(eventtime, degree = 8)6:trt   1.0980     7.9107   0.139   0.8896    
#> bs(eventtime, degree = 8)7:trt   2.4548     3.5013   0.701   0.4832    
#> bs(eventtime, degree = 8)8:trt   1.8150     1.5618   1.162   0.2452    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 8142.9  on 13364  degrees of freedom
#> Residual deviance: 7861.6  on 13347  degrees of freedom
#> AIC: 7897.6
#> 
#> Number of Fisher Scoring iterations: 7

# using casebase hazardPlot ------------------------------------------------

t1_hazardplot <- hazardPlot(object = mod_cb, newdata = data.frame(trt = 1),
                       ci.lvl = 0.95, ci = TRUE, lty = 1, line.col = 1, lwd = 2)
#> Warning in bs(eventtime, degree = 8L, knots = numeric(0), Boundary.knots =
#> c(5.14404382556677e-05, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases

t0_hazardplot <- hazardPlot(object = mod_cb, newdata = data.frame(trt = 0), ci = TRUE,
           ci.lvl = 0.95, add = TRUE, lty = 2, line.col = 2, lwd = 2)
#> Warning in bs(eventtime, degree = 8L, knots = numeric(0), Boundary.knots =
#> c(5.14404382556677e-05, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
legend("topleft", c("trt=1","trt=0"),lty=1:2,col=1:2,bty="y", lwd = 2)

# using visreg ------------------------------------------------------------

mod_cb$data$offset <- 0
tt <- visreg(mod_cb,
       xvar = "eventtime",
       by = "trt",
       scale = "linear",
       trans = exp,
       data = mod_cb$data,
       plot = T,
       overlay = TRUE)

# get the predicted hazards by trt status and order the
# times so you can plot them
t1 <- tt$fit[which(tt$fit$trt==1),]
t0 <- tt$fit[which(tt$fit$trt==0),]
t1_visreg <- t1[order(t1$eventtime),]
t0_visreg <- t0[order(t1$eventtime),]



# ggpredict ---------------------------------------------

# this stops predictions at eventtimes <0.5
dat1 <- ggpredict(mod_cb, terms = c("eventtime [all]","trt","offset [0]"))
#> Warning: Some model terms could not be found in model data. You probably need to
#> load the data into the environment.
#> Warning: Some model terms could not be found in model data. You probably need to
#> load the data into the environment.
dat1$x %>% range
#> [1] 0.000 0.393

# this works, but produces slightly different results than visreg or hazardPlot
dat <- ggpredict(mod_cb, terms = c("eventtime [0:5 by=.1]","trt","offset [0]"))
#> Warning: Some model terms could not be found in model data. You probably need to
#> load the data into the environment.

#> Warning: Some model terms could not be found in model data. You probably need to
#> load the data into the environment.
plot(dat)

# get the predicted hazards by trt status. they are already in order
# wih ggpredict
t1_ggpredict <- dat$predicted[which(dat$group==1)]
t0_ggpredict <- dat$predicted[which(dat$group==0)]



# compare hazard ratios for treatment -----------------------------------------

plot(t0_hazardplot$eventtime, t1_hazardplot$predictedhazard / t0_hazardplot$predictedhazard,
     type = "l", col = 1, lty = 1, ylab = "hazard ratio for trt=1 vs. trt=0", lwd = 2, ylim = c(0,1.1))
abline(a=1, b=0, col = "grey80")
lines(t0_visreg$eventtime, t1_visreg$visregFit / t0_visreg$visregFit,
      type = "l", col = 2, lty = 2, lwd = 2)
lines(unique(dat$x), t1_ggpredict / t0_ggpredict,
      type = "l", col = 3, lty = 3, lwd = 2)

legend("bottomright", c("hazardPlot","visreg","ggpredict"),
       lty=1:3,col=1:3,bty="y", lwd = 2)

Created on 2020-03-10 by the reprex package (v0.3.0)

turgeonmaxime commented

One thing I prefer about ggpredict and visreg over hazardPlot, the confidence region are transparent and coloured the same way as the lines, so they're easier to read. On the other hand, that's not hard to fix in hazardPlot.

Of these two, I prefer the syntax of visreg, I think the syntax terms = c("eventtime [all]","trt","offset [0]") looks akward, and of course we don't really know why it doesn't quite give the same answer. I can look for a minimal reproducible example (my gut feeling is that the offset is probably driving this difference).

We can definitely make internal changes so that these functions are easier to use. Two things come to mind:

  • Labelling the original data and the sampled data correctly so that both ggpredict and visreg use the "right" data.
  • We could put the offset at 0 internally, and keep track of it in a separate slot of the fit object (I don't think we ever reuse it, but we definitely don't want to lose track of it).

In other words, I can definitely see a world where all three options can be used: hazardPlot as a simple functionality of our package, and for more complex uses we can make sure our package plays nicely with both visreg and ggpredict.

Jesse-Islam commented

in terms of support, it seems ggpredict has a main developer that commits often in comparison to visreg. visreg is still being added to and has been around longer than ggpredict, by about 6 years.

In terms of package support long term, I think it make sense to just wrap one of their functions. It looks like there are a few warnings associated with ggpredict, while visreg seems to run with no issues. With this in mind I think wrapping visreg would be the best bet in terms of time commitment.

turgeonmaxime commented

I agree that it's important to have in mind the long-term support of these packages before committing any resources to supporting their use. But it looks like both were updated in the last 3-6 months, so they both look well-maintained.

However, I don't think we need to wrap anything. If we fix the offset internally, then users just need to call plot(ggpredict(mod_cb, terms = c("eventtime [all]","trt","offset [0]"))) or visreg(mod_cb, xvar = "eventtime", by = "trt", scale = "linear", trans = exp, data = mod_cb$data, plot = T, overlay = TRUE) directly. The other advantage is we wouldn't need to include either package as a dependency.

  • 👍1
sahirbhatnagar commented

Thanks for the answers.

I'm confident that both visreg and ggeffects use the data slot of the fit object. I will focus on making casebase compatible with visreg for now.

Like you said, the only change is to set mod_cb$data$offset <- 0, but then store it in its own slot because we don't want to lose this information. Next steps:

  • Store offset value in its own slot in the object returned by fitSmoothHazard and see if that breaks anything in the package.

  • Try visreg with family="glmnet" and family="gbm"

  • Create vignette to show time dependent hazard ratios with visreg

In order to build the vignette, do we have to put visreg in Suggests ?

turgeonmaxime commented

In order to build the vignette, do we have to put visreg in Suggests ?

Yes, and we should probably also make it fail "gracefully" if the package isn't installed. Something like this in the setup chunk should be enough:

knitr::opts_chunk$set(
    val = requireNamespace("visreg", quietly = TRUE)
)
sahirbhatnagar commented

I have tried a whole bunch of options for how to move forward with hazard ratio plots and hazard plots. Among the things I tried are:

I think after much trial, visreg is the way to go for the following reasons:

  • it plays nice with glm objects and any spline function as long as it comes from a formula interface
  • it's the easiest package to understand in terms of what's happening behind the scenes, i.e., it's clear what the values of the confounders are when plotting time-dependent hazard ratios for a treatment group
  • very little change, if any, is required to the existing fitSmoothHazard output. HOWEVER, it is crucial to note that visreg uses the residuals from the model which are based on the augmented casebase dataset. not sure what the implications are for this (if any).

Plotting time-dependent hazard ratios will require some work, since it doesn't naturally come out of visreg, however I think this is hard to do programtically anyways (for multi-category treatment variables which interact with time). I will make sure a few simple cases work (e.g. for a binary predictor and continuous predictor interacting with time).

Going forward I need to:

  • Add visreg to suggests
  • Create a wrapper function around visreg to plot hazard functions and hazard ratios
  • Figure out what the best way forward is for the API of this plotting functionality with casebase. We had talked about creating a casebase class in #98 and then create a plot.casebase method. But we don't want to override the plot method for glm objects. NextMethod is an option but it's complicated to understand.
  • Create vignette with hazard plots, hazard ratios and potentially cumulative incidence curves.

Here is an example of how we can use visreg for time dependent hazard ratios and hazard functions:

pacman::p_load(rstpm2) # for brcancer data
pacman::p_load(magrittr)
pacman::p_load(visreg)
devtools::load_all("~/git_repositories/casebase/")
#> Loading casebase
#> 
#> Attaching package: 'testthat'
#> The following objects are masked from 'package:magrittr':
#> 
#>     equals, is_less_than, not
#> See example usage at http://sahirbhatnagar.com/casebase/
data("brcancer")
brcancer <- transform(brcancer, recyear=rectime / 365.24)
brcancer <- transform(brcancer, hormon=factor(hormon))
mod_cb <- fitSmoothHazard(censrec ~ hormon*bs(recyear, df = 5),
                          data = brcancer,
                          time = "recyear")

par(mfrow=c(1,2))
t0_hazardplot <- hazardPlotNew(object = mod_cb, newdata = data.frame(hormon = "0"), ci = F,
                            ci.lvl = 0.95, lty = 2, line.col = 2, lwd = 2)

t1_hazardplot <- hazardPlotNew(object = mod_cb, newdata = data.frame(hormon = "1"),
                               ci.lvl = 0.95, ci = F, lty = 1, line.col = 1, lwd = 2, add = T)

# this uses the augmented casebase dataset which has a non-zero offset.
# need to set it to 0 here
tt <- visreg(mod_cb,
             xvar = "recyear",
             # type = "contrast",
             by = "hormon",
             scale = "linear",
             cond = list(offset = 0),
             trans = exp,
             # data = mod_cb$originalData,
             plot = T,
             alpha = 1,
             overlay = TRUE)

# get the predicted hazards by trt status and order the
# times so you can plot them
t1 <- tt$fit[which(tt$fit$hormon==1),]
t0 <- tt$fit[which(tt$fit$hormon==0),]
t1_visreg <- t1[order(t1$recyear),]
t0_visreg <- t0[order(t1$recyear),]


# dev.off()
par(mfrow = c(1,2))
plot(t0_hazardplot$recyear, t1_hazardplot$predictedhazard / t0_hazardplot$predictedhazard,
     type = "l", col = 1, lty = 1, ylab = "hazard ratio for hormon=1 vs. hormon=0", lwd = 2, ylim = c(0,1.1))
abline(a=1, b=0, col = "grey80")
lines(t0_visreg$recyear, t1_visreg$visregFit / t0_visreg$visregFit,
      type = "l", col = 2, lty = 2, lwd = 2)
legend("bottomleft", c("hazardPlot","visreg"),
       lty=1:2,col=1:2,bty="y", lwd = 2)

# this is type contrast: E(y|x=1, z) - E(y|x=0, z)
# where z are the other covariates not of interest
# I thought it would give me just one line, but its giving me two
kk <- visreg(mod_cb,
             xvar = "recyear",
             type = "contrast",
             by = "hormon",
             scale = "linear",
             cond = list(offset = 0),
             trans = exp,
             plot = T,
             alpha = 1,
             overlay = T)

Created on 2020-03-31 by the reprex package (v0.3.0)

turgeonmaxime commented

I know I've mentioned this before, but I'm not sure we answered it completely: is there any value in wrapping visreg inside one of our functions, as opposed to changing the output of fitSmoothHazard so that we can "just" call visreg(mod_cb)? Of course, if it means changing the output of fitSmoothHazard significantly, or in a way that breaks other functionalities, I understand the advantage of writing a wrapper. But otherwise...?

sahirbhatnagar commented

is there any value in wrapping visreg inside one of our functions, as opposed to changing the output of fitSmoothHazard so that we can "just" call visreg(mod_cb)? Of course, if it means changing the output of fitSmoothHazard significantly, or in a way that breaks other functionalities, I understand the advantage of writing a wrapper.

Yes you're right about the hazard function. We would just need to set the offset to 0 in the data slot returned by fitSmoothHazard which is the augmented casebase dataset.

But otherwise...?

I can't figure out how to plot the hazard ratio using only a simple call to visreg. As you can see in the example above, I have to store the results from visreg and then take the appropriate ratio of the hazards. I thought a contrast plot might do the trick, but it's still plotting two lines, one for each level of the binary predictor, which I dont understand why based on his formula:

image

where x_0 is a reference value and x is the value of interest

turgeonmaxime commented

Ah, good point, I see how the wrapper can save a lot of headache to users.

sahirbhatnagar commented

closed by #117