flexsurv-examples.R
In flexsurv: Flexible Parametric Survival and Multi-State Models

### R code from vignette source 'flexsurv-examples.Rnw'

###################################################
### code chunk number 1: flexsurv-examples.Rnw:39-59
###################################################
library(flexsurv)

hgengammaPH <- function(x, dummy, mu=0, sigma=1, Q){
    dummy * hgengamma(x=x, mu=mu, sigma=sigma, Q=Q)
}

HgengammaPH <- function(x, dummy, mu=0, sigma=1, Q){
    dummy * Hgengamma(x=x, mu=mu, sigma=sigma, Q=Q)
}

custom.gengammaPH <- list(name="gengammaPH", 
                          pars=c("dummy","mu","sigma","Q"), location="dummy",
                          transforms=c(log, identity, log, identity),
                          inv.transforms=c(exp, identity, exp, identity),
                          inits=function(t){
                              lt <- log(t[t>0])
                              c(1, mean(lt), sd(lt), 0)
                          })
fs7 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data=bc, 
                   dist=custom.gengammaPH, fixedpars=1)


###################################################
### code chunk number 2: flexsurv-examples.Rnw:84-127
###################################################

fs2 <- flexsurvreg(Surv(recyrs, censrec) ~ group + sigma(group), 
                   data=bc, dist="gengamma")
B <- 5000
t <- seq(0.1, 8, by=0.1)

hrAFT.est <-
    summary(fs2, t=t, type="hazard",
            newdata=data.frame(group="Medium"),ci=FALSE)[[1]][,"est"] /
    summary(fs2, t=t, type="hazard",
            newdata=data.frame(group="Good"),ci=FALSE)[[1]][,"est"]
pars <- normboot.flexsurvreg(fs2, B=B, newdata=data.frame(group=c("Good","Medium")))
hrAFT <- matrix(nrow=B, ncol=length(t))
for (i in seq_along(t)){ 
    haz.medium.rep <- do.call(hgengamma, c(list(t[i]), as.data.frame(pars[[2]])))
    haz.good.rep <- do.call(hgengamma, c(list(t[i]), as.data.frame(pars[[1]])))
    hrAFT[,i] <- haz.medium.rep / haz.good.rep
}
hrAFT <- apply(hrAFT, 2, quantile, c(0.025, 0.975))

hrPH.est <-
    summary(fs7, t=t, type="hazard",
            newdata=data.frame(group="Medium"),ci=FALSE)[[1]][,"est"] /
    summary(fs7, t=t, type="hazard",
            newdata=data.frame(group="Good"),ci=FALSE)[[1]][,"est"]
pars <- normboot.flexsurvreg(fs7, B=B, newdata=data.frame(group=c("Good","Medium")))
hrPH <- matrix(nrow=B, ncol=length(t))
for (i in seq_along(t)){ 
    haz.medium.rep <- do.call(hgengammaPH, c(list(t[i]), as.data.frame(pars[[2]])))
    haz.good.rep <- do.call(hgengammaPH, c(list(t[i]), as.data.frame(pars[[1]])))
    hrPH[,i] <- haz.medium.rep / haz.good.rep
}
hrPH <- apply(hrPH, 2, quantile, c(0.025, 0.975))

plot(t, hrAFT[1,], type="l", ylim=c(0, 10), col="red", xlab="Years",
     ylab="Hazard ratio (Medium / Good)", lwd=1, lty=2)
lines(t, hrAFT[2,], col="red", lwd=1, lty=2)
lines(t, hrPH[1,], col="darkgray", lwd=1, lty=2)
lines(t, hrPH[2,], col="darkgray", lwd=1, lty=2)
lines(t, hrAFT.est, col="red", lwd=2)
lines(t, hrPH.est, col="darkgray", lwd=2)
legend("topright", lwd=c(2,2), col=c("red","darkgray"), bty="n",
       c("Generalized gamma: standard AFT", "Generalized gamma: proportional hazards"))


###################################################
### code chunk number 3: flexsurv-examples.Rnw:132-146 (eval = FALSE)
###################################################
## nd <- data.frame(group=c("Good","Medium"))
## hr_aft <- hr_flexsurvreg(fs2, t=t, newdata=nd)
## hr_ph <- hr_flexsurvreg(fs7, t=t, newdata=nd)
## 
## plot(t, hr_aft$lcl, type="l", ylim=c(0, 10), col="red", xlab="Years",
##      ylab="Hazard ratio (Medium / Good)", lwd=1, lty=2)
## lines(t, hr_aft$ucl, col="red", lwd=1, lty=2)
## lines(t, hr_aft$est, col="red", lwd=2)
## 
## lines(t, hr_ph$lcl, col="darkgray", lwd=1, lty=2)
## lines(t, hr_ph$ucl, col="darkgray", lwd=1, lty=2)
## lines(t, hr_ph$est, col="darkgray", lwd=2)
## legend("topright", lwd=c(2,2), col=c("red","darkgray"), bty="n",
##        c("Generalized gamma: standard AFT", "Generalized gamma: proportional hazards"))


###################################################
### code chunk number 4: flexsurv-examples.Rnw:168-170
###################################################
summary(fs2, type="rmst", t=100, tidy=TRUE)
summary(fs2, fn=rmst_gengamma, t=100, tidy=TRUE)


###################################################
### code chunk number 5: flexsurv-examples.Rnw:174-176
###################################################
est <- fs2$res[,"est"]
rmst_gengamma(t=100, mu=est[1], sigma=est[2], Q=est[3])


###################################################
### code chunk number 6: flexsurv-examples.Rnw:180-181
###################################################
summary(fs2, type="median")


###################################################
### code chunk number 7: flexsurv-examples.Rnw:211-231
###################################################
if (require("TH.data")){

GBSG2 <- transform(GBSG2,
                   X1a=(age/50)^-2,
                   X1b=(age/50)^-0.5,
                   X4=tgrade %in% c("II","III"),
                   X5=exp(-0.12*pnodes),
                   X6=(progrec+1)^0.5
                   )
(progc <- coxph(Surv(time, cens) ~ horTh + X1a + X1b + X4 + 
                  X5 + X6, data=GBSG2))
(prog3 <- flexsurvspline(Surv(time, cens) ~ horTh + X1a + X1b + X4 + 
                           X5 + X6, k=3, data=GBSG2))
predc <- predict(progc, type="lp")
progc <- cut(predc, quantile(predc, 0:3/3))
predf <- model.matrix(prog3) %*% prog3$res[-(1:5),"est"]
progf <- cut(predf, quantile(predf, 0:3/3))
table(progc, progf)

}


###################################################
### code chunk number 8: flexsurv-examples.Rnw:251-255
###################################################
set.seed(1)
nsim <- 10000
onsetday <- runif(nsim, 0, 30) 
deathday <- onsetday + rgamma(nsim, shape=1.5, rate=1/10)


###################################################
### code chunk number 9: flexsurv-examples.Rnw:265-269
###################################################
datt <- data.frame(delay = deathday - onsetday,
                   event = rep(1, nsim),
                   rtrunc = 40 - onsetday)
datt <- datt[datt$delay < datt$rtrunc, ] 


###################################################
### code chunk number 10: flexsurv-examples.Rnw:279-282
###################################################
fitt <- flexsurvreg(Surv(delay, event) ~ 1, data=datt, rtrunc = rtrunc, dist="gamma")
fitt
summary(fitt, t=1, fn = mean_gamma)