Extracting baseline hazards from stpm2-object

Question

I need to extract the baseline hazards from a general survival model (GSM) that I've constructed using the rstpm2-package (a conversion of the stpm2 module in stata).

using the data in the rstpm2-package let's use this as an example:

library(rstpm2)

gsm <- stpm2(Surv(rectime,censrec==1)~hormon, data=brcancer, df=3)

sum.gsm <- summary(gsm)

So I've noticed that the summary has an element named bhazard:

sum.gsm@args$bhazard

However it seems to be filled with zeroes and holds one value per patient. As far as I understand the baseline hazard should consist of one hazard for every time-point in the data.

Does anyone have any experience that could be of assistance

Answer 1

You can use the plot and lines methods to plot survival and a number of other predictions. For example:

plot(gsm, newdata=data.frame(hormon=0))

Plot

Note that you need to specify the newdata argument. For more general plots, you can get predictions on a time grid with full covariates with standard errors using:

out <- predict(gsm,newdata=data.frame(hormon=0:1), grid=TRUE, 
               full=TRUE, se.fit=TRUE)

Then you could use ggplot2 or lattice to plot the intervals. For example,

library(ggplot2)
ggplot(out, aes(x=rectime, y=Estimate, ymin=lower, ymax=upper, 
                fill=factor(hormon))) +
    geom_ribbon(alpha=0.6) + geom_line() 

Plot

Edit: to predict survival at specific times, you can include the times in the newdata argument:

predict(gsm, newdata=data.frame(hormon=0,rectime=1:365))

(Note that survival is defined in terms of log(time), hence I have excluded rectime=0.)