How to Fit Piece-Wise Exponential Model in r?

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

• Have you taken a look at Generalized Additive Models? (web.stanford.edu/~hastie/Papers/gam.pdf) – Jon Sep 22 '16 at 16:39
• Take a look at the pch package on CRAN cran.r-project.org/web/packages/pch/pch.pdf – boscovich Sep 22 '16 at 18:52
• 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? – A-L Sep 22 '16 at 22:10

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

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
library(pammtools)
library(survival)
library(ggplot2)
theme_set(theme_bw())
library(mgcv)

data(tumor)

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

# 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")) %>%

• @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. – adibender Sep 9 '18 at 16:23