0

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

Nick Cox
  • 35,529
  • 6
  • 31
  • 47
Jarn Schöber
  • 309
  • 1
  • 8
  • I've removed the Stata tag. I can't see why this is of interest to Stata users, unless they happen to be R users, in which case that tag suffices. – Nick Cox Oct 18 '19 at 09:46
  • I figured that as the package is originally developed for stata, they might hold some knowledge as to how this is done there. – Jarn Schöber Oct 18 '19 at 10:33
  • Perhaps, but my experience is that while the main results are or should bevthe same, the details of any translation from one language to another are often quite different. – Nick Cox Oct 18 '19 at 11:07

1 Answers1

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.)

Mark Clements
  • 120
  • 1
  • 5
  • I think this brings me one step closer to what I really want to do. Do you know if there's a way to extract the baseline for specific timepoints? For example seq(0,365,1)? – Jarn Schöber Jan 28 '20 at 08:59