1
$\begingroup$

I am experimenting Piece-Wise Exponential Model for survival data. What is the best way to fit the model? Is there an existing package?

$\endgroup$
  • $\begingroup$ Have you taken a look at Generalized Additive Models? (web.stanford.edu/~hastie/Papers/gam.pdf) $\endgroup$ – Jon Sep 22 '16 at 16:39
  • $\begingroup$ Take a look at the pch package on CRAN cran.r-project.org/web/packages/pch/pch.pdf $\endgroup$ – boscovich Sep 22 '16 at 18:52
  • $\begingroup$ Thanks. I am actually interested in a piecewise constant baseline proportional hazard model. I want the effect of covariates stays the same cross intervals, only the baseline hazard to be a piecewise spline. In other words, is there a better way to estimate the piecewise baseline hazard for Cox model? $\endgroup$ – A-L Sep 22 '16 at 22:10
1
$\begingroup$

in order to get a satisfying answer, you should provide some details about the context and problem for which you want to run a piecewise Exponential model. There are several packages which might address your problem and each of them has its own peculiarity. You may want to look at the CRAN Task View on Survival Analysis where you can have several references.

By far, the most know R package to run survival analysis is survival. This is a huge package which contains dozens of routines. As you pointed out in the comment, you can run a Cox proportional model through the function coxph(). There are caveats when you do this though. For instance, coxph() estimates a semi-parametric Cox Proportional Model as given by the following equation:

$$ \lambda(t) = \lambda_0(t) \exp\left[\boldsymbol\beta^T(t) \mathbf{z}(t)\right] $$ where $\lambda_0(t)$ is the baseline hazard on which you don't take any assumptions. This is estimated non-parametrically then, while the $\exp\left[\boldsymbol\beta^T(t) \mathbf{z}(t)\right]$ is estimated through Maximum Likelihood so that it gives you the effects of every covariate you put into the model. If you want to model your data this way, then coxph() is the way to go.

survival also provides a function which estimates the same model, but assuming a fully-parametric schema. That is, you assume a probability distribution on $\lambda_0(t)$. This can be achieved through survreg(). Though, I suggest to use a more flexible package called flexsurv which extends what survival does in the context of parametric models applied to survival data.

All these packages treat any time-dependent covariate as a piece-wise constant function. This means that, you do not observe a change in the value of a given covariate except when the event of interest does occur.

Hope this helped.
Franz

| cite | improve this answer | |
$\endgroup$
1
$\begingroup$

You basically just need to transform the data to a suitable format. Then you can estimate the piece-wise constant baseline hazard using penalized splines. If you assume that the true, underlying hazard is smooth, you can improve the approximation by icreasing the number of intervals/reducing interval lengths (the cut argument below). More info (with R examples) on estimating the baseline hazard with PEMs/PAMMs can be found here: https://adibender.github.io/pammtools/articles/baseline.html

An example with R below:

## install the pammtools package for data transformation and example data
devtools::install_github("adibender/pammtools")
library(pammtools)
library(survival)
library(ggplot2)
theme_set(theme_bw())
library(mgcv)

# load example data
data(tumor)

## transform to piece-wise exponential data (PED)
ped_tumor <- as_ped(Surv(days, status) ~ ., data = tumor, cut = seq(0, 3000, by = 100))

# Fit the Piece-wise-exponential Additive Model (PAM) instead of PEM
# with piece-wise constant hazards, term s(tend), and constant covariate effects:
pam <- mgcv::gam(ped_status ~ s(tend) + age + sex + complications, 
  data = ped_tumor, family = poisson(), offset = offset)
summary(pam)
## plot baseline hazard: 
ped_tumor %>% 
  make_newdata(tend=unique(tend), age = c(0), sex = c("male"), complications=c("no")) %>% 
  add_hazard(pam) %>% 
  ggplot(aes(x = tstart, y = hazard)) + 
  geom_stepribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2) +
  geom_step()
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I used your code to perform a similar analysis. I obtained hazard rate from add_hazard. How would you interpret hazard = 0.02 given that age = 20, sex = female, and tend = 100? I interpreted this as the expected number of deaths the female cohort at age 20 would experience death during the interval 0-100 days is 0.02 but that doesn't make much sense...Thank you. $\endgroup$ – Meo Sep 7 '18 at 21:55
  • $\begingroup$ @Meo add_hazard returns the hazard rate, which is the probability of having an event (given covariate specification) in interval 0-100, divided by the length of the interval. You usually interpret the hazard rate by comparison with hazard rates with other covariate specifications (e.g. males vs. femaies to see which covariate specification has higher hazard), for example by calculating the hazard ratio. $\endgroup$ – adibender Sep 9 '18 at 16:23
  • $\begingroup$ Hi adibender, I was trying to replicate the result for a paper where they plotted the hazard rates but using a different dataset. The confidence interval in my codes exceeded 1, which is why I am having a hard time explaining it in my report. Hence my question. Thank you so much for your reply and your R package- very helpful! $\endgroup$ – Meo Sep 10 '18 at 22:14
  • $\begingroup$ @Meo In general, the hazard rate ca be > 1, thus upper CI with values of > 1 can also occurr. Glad I could help. If you have other questions or want to collaborate on survival projects you can also drop me an email (see contact details on the package homepage) $\endgroup$ – adibender Sep 11 '18 at 6:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.