Using `flexsurv` to fit a piecewise exponential model

[elong0527]

Could anyone provide some hint if I can fit a piecewise exponential model using flexsurv?

I feel I can do it using flexsurvspline, but not sure how to define the parameters for a piecewise linear model when scale=hazard

[chjackson]

There's no way of doing this using flexsurvspline, I think, as that fits a cubic spline basis. At the moment you'd have to use flexsurvreg and construct a custom piecewise exponential distribution with a fixed number of pieces. Easiest way would be via the hazard function and cumulative hazard functions.

Here's a quick example, which isn't tested much. Ideally there should be something built in to construct the functions automatically given the vector of knots.

hpwexp <- Vectorize(function(t, rate0, rate1, rate2){
  rate <- c(rate0, rate1, rate2)
  knots <- sort(knots)
  rate[findInterval(t, knots) + 1]
})

Hpwexp <- Vectorize(function(t, rate0, rate1, rate2){
  k0 <- c(0,knots)
  rate <- c(rate0, rate1, rate2)
  int <- findInterval(t, k0)
  dk <- diff(k0)
  cumrate <- c(0, cumsum(dk*rate[-length(rate)]))    # cumulative rate up to each knot
  cumrate[int] + rate[int]*(t - k0[int])             # add on the remainder for each t
})

knots <- c(2,6)
nk <- length(knots)
custom_pwexp <- 
  list(name = "pwexp", 
       pars = paste0("rate", 0:nk),
       location = "rate0",
       transforms = rep(list(log), nk+1),
       inv.transforms = rep(list(exp), nk+1),
       inits =  function(t,aux) {lam <- sr.exp.inits(t,aux); rep(lam, nk+1)})

pwexpfs <- flexsurvreg(Surv(recyrs, censrec) ~ 1, data=bc, dist=custom_pwexp)

plot(pwexpfs, ci=FALSE)
plot(pwexpfs, type="hazard", ci=FALSE, ylim=c(0,0.5))

[elong0527]

Thanks for the clarification.

I will cross check the function with the pch package.

Agree, it would be nice to add a built in function to handle piecewise exponential model that become popular in the study design literatures.

[Philip-Cooney]

I think the code should be modified to set the location to a "dummy" variable so that covariates can be introduced (as per Chris's example).

The optimization doesn't really work as it takes quite a while and doesn't give finite confidence intervals. Are there additional distribution or density functions that should be provided to facilitate convergence?

hpwexp <- Vectorize(function(t, dummy, rate0, rate1, rate2){
  rate <- c(rate0, rate1, rate2)
  knots <- sort(knots)
  dummy*(rate[findInterval(t, knots) + 1])
})

Hpwexp <- Vectorize(function(t,dummy, rate0, rate1, rate2){
  k0 <- c(0,knots)
  rate <- c(rate0, rate1, rate2)
  int <- findInterval(t, k0)
  dk <- diff(k0)
  cumrate <- c(0, cumsum(dk*rate[-length(rate)]))    # cumulative rate up to each knot
  dummy*(cumrate[int] + rate[int]*(t - k0[int]))             # add on the remainder for each t
})

knots <- c(0.45,1.869)
nk <- length(knots)

custom_pwexp <- 
  list(name = "pwexp", 
       pars = c("dummy", paste0("rate", 0:nk)),
       location = "dummy",
       transforms = rep(list(log), nk+2), #Add one more for dummy param
       inv.transforms = rep(list(exp), nk+2), #Add one more for dummy param
       inits =  function(t,aux) {lam <- 0.2; c(1,rep(lam, nk+1))})

pwexpfs <- flexsurvreg(Surv(recyrs, censrec) ~ as.factor(group), data=bc, dist=custom_pwexp)
plot(pwexpfs)

[chjackson]

For putting covariates in, I'd consider parameterising as (bhaz, hr2, hr3, ...) where bhaz is the hazard in period 1, and hr2, hr3 are the hazard ratios between period 2, 3 etc. and period 1. bhaz is the location parameter and putting covariates on that gives a proportional hazards model. This may or may not facilitate convergence - if the CIs can't be determined that's often just because the model isn't identifiable given the data.