## ----echo=FALSE------------------------------------------- options(width=60) options(prompt="R> ") library(knitr) opts_chunk$set(fig.path="flexsurv-") render_sweave() ## ----results='hide'--------------------------------------- library("flexsurv") ## --------------------------------------------------------- head(bc, 2) ## --------------------------------------------------------- fs1 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data = bc, dist = "weibull") ## --------------------------------------------------------- fs1 ## --------------------------------------------------------- survreg(Surv(recyrs, censrec) ~ group, data = bc, dist = "weibull") ## ----results='hide',eval=FALSE---------------------------- # library(eha) # aftreg(Surv(recyrs, censrec) ~ group, data = bc, dist = "weibull") ## --------------------------------------------------------- fs2 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data = bc, dist = "gengamma") fs3 <- flexsurvreg(Surv(recyrs, censrec) ~ group + sigma(group), data = bc, dist = "gengamma") ## ----eval=FALSE------------------------------------------- # fs3 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data = bc, # anc = list(sigma = ~ group), dist = "gengamma") ## ----surv,include=FALSE----------------------------------- cols <- c("#E495A5", "#86B875", "#7DB0DD") # from colorspace::rainbow_hcl(3) plot(fs1, col = cols[2], lwd.obs = 2, xlab = "Years", ylab = "Recurrence-free survival") lines(fs2, col = cols[3], lty = 2) lines(fs3, col = cols[3]) text(x=c(2,3.5,4), y=c(0.4, 0.55, 0.7), c("Poor","Medium","Good")) legend("bottomleft", col = c("black", cols[2], cols[3], cols[3]), lty = c(1, 1, 2, 1), bty = "n", lwd = rep(2, 4), c("Kaplan-Meier", "Weibull", "Generalized gamma (AFT)", "Generalized gamma (time-varying)")) ## ----haz,include=FALSE------------------------------------ plot(fs1, type = "hazard", col = cols[2], lwd.obs = 2, max.time=6, xlab = "Years", ylab = "Hazard") lines(fs2, type = "hazard", col = cols[3], lty = 2) lines(fs3, type = "hazard", col = cols[3]) text(x=c(2,2,2), y=c(0.3, 0.13, 0.05), c("Poor","Medium","Good")) legend("topright", col = c("black", cols[2], cols[3], cols[3]), lty = c(1, 1, 2, 1), bty = "n", lwd = rep(2, 4), c("Kernel density estimate", "Weibull", "Gen. gamma (AFT)", "Gen. gamma (time-varying)")) ## --------------------------------------------------------- median.weibull <- function(shape, scale) { qweibull(0.5, shape = shape, scale = scale) } summary(fs1, fn = median.weibull, t = 1, B = 10000) ## ----eval=FALSE------------------------------------------- # custom.llogis <- list(name = "llogis", pars = c("shape", "scale"), # location = "scale", # transforms = c(log, log), # inv.transforms = c(exp, exp), # inits = function(t){ c(1, median(t)) }) # fs4 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data = bc, # dist = custom.llogis) ## ----eval=FALSE------------------------------------------- # dmakeham3 <- function(x, shape1, shape2, scale, ...) { # dmakeham(x, shape = c(shape1, shape2), scale = scale, ...) # } # pmakeham3 <- function(q, shape1, shape2, scale, ...) { # pmakeham(q, shape = c(shape1, shape2), scale = scale, ...) # } ## ----eval=FALSE------------------------------------------- # dmakeham3 <- Vectorize(dmakeham3) # pmakeham3 <- Vectorize(pmakeham3) ## ----eval=FALSE------------------------------------------- # pmakeham3(c(0, 1, 1, Inf), 1, c(1, 1, 2, 1), 1) ## ----echo=FALSE------------------------------------------- options(warn=-1) ## --------------------------------------------------------- hweibullPH <- function(x, shape, scale = 1, log = FALSE){ hweibull(x, shape = shape, scale = scale ^ {-1 / shape}, log = log) } HweibullPH <- function(x, shape, scale = 1, log = FALSE){ Hweibull(x, shape = shape, scale = scale ^ {-1 / shape}, log = log) } custom.weibullPH <- list(name = "weibullPH", pars = c("shape", "scale"), location = "scale", transforms = c(log, log), inv.transforms = c(exp, exp), inits = function(t){ c(1, median(t[t > 0]) / log(2)) }) fs6 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data = bc, dist = custom.weibullPH) fs6$res["scale", "est"] ^ {-1 / fs6$res["shape", "est"]} - fs6$res["groupMedium", "est"] / fs6$res["shape", "est"] - fs6$res["groupPoor", "est"] / fs6$res["shape", "est"] ## ----echo=FALSE------------------------------------------- options(warn=0) ## --------------------------------------------------------- sp1 <- flexsurvspline(Surv(recyrs, censrec) ~ group, data = bc, k = 1, scale = "odds") ## --------------------------------------------------------- sp2 <- flexsurvspline(Surv(recyrs, censrec) ~ group + gamma1(group), data = bc, k = 1, scale = "odds") ## ----splinehaz,include=FALSE------------------------------ plot(sp1, type = "hazard", col=cols[3], ylim = c(0, 0.5), xlab = "Years", ylab = "Hazard") lines(sp2, type = "hazard", col = cols[3], lty = 2) lines(fs2, type = "hazard", col = cols[2]) text(x=c(2,2,2), y=c(0.3, 0.15, 0.05), c("Poor","Medium","Good")) legend("topright", col = c("black",cols[c(3,3,2)]), lty = c(1,1,2,1), lwd = rep(2,4), c("Kernel density estimate","Spline (proportional odds)", "Spline (time-varying)","Generalized gamma (time-varying)")) ## --------------------------------------------------------- sp3 <- flexsurvspline(Surv(recyrs, censrec) ~ group, data = bc, k = 1, scale = "hazard") sp3$res[c("groupMedium", "groupPoor"), c("est", "se")] cox3 <- coxph(Surv(recyrs, censrec) ~ group, data = bc) coef(summary(cox3))[ , c("coef", "se(coef)")] ## --------------------------------------------------------- sp4 <- flexsurvspline(Surv(recyrs, censrec) ~ group + gamma1(group) + gamma2(group), data = bc, k = 1, scale = "hazard") ## ----eval=FALSE------------------------------------------- # flexsurvspline(Surv(recyrs, censrec) ~ group + gamma1(group) + # gamma2(group) + treat + gamma1(treat), # data = bc, k = 1, scale = "hazard") ## --------------------------------------------------------- res <- t(sapply(list(fs1, fs2, fs3, sp1, sp2, sp3, sp4), function(x)rbind(-2 * round(x$loglik,1), x$npars, round(x$AIC,1)))) rownames(res) <- c("Weibull (fs1)", "Generalized gamma (fs2)", "Generalized gamma (fs3)", "Spline (sp1)", "Spline (sp2)", "Spline (sp3)", "Spline (sp4)") colnames(res) <- c("-2 log likelihood", "Parameters", "AIC") ## ----size="scriptsize"------------------------------------ res ## --------------------------------------------------------- gamma <- sp1$res[c("gamma0", "gamma1", "gamma2"), "est"] 1 - psurvspline(5, gamma = gamma, knots = sp1$knots) ## --------------------------------------------------------- pfn <- unroll.function(psurvspline, gamma = 0:2) 1 - pfn(5, gamma0 = gamma[1], gamma1 = gamma[2], gamma2 = gamma[3], knots = sp1$knots) ## --------------------------------------------------------- hsurvspline.lh <- function(x, gamma, knots){ if(!is.matrix(gamma)) gamma <- matrix(gamma, nrow = 1) lg <- nrow(gamma) nret <- max(length(x), lg) gamma <- apply(gamma, 2, function(x)rep(x, length = nret)) x <- rep(x, length = nret) loghaz <- rowSums(basis(knots, log(x)) * gamma) exp(loghaz) } ## --------------------------------------------------------- hsurvspline.lh3 <- unroll.function(hsurvspline.lh, gamma = 0:2) ## --------------------------------------------------------- custom.hsurvspline.lh3 <- list( name = "survspline.lh3", pars = c("gamma0", "gamma1", "gamma2"), location = c("gamma0"), transforms = rep(c(identity), 3), inv.transforms = rep(c(identity), 3) ) dtime <- log(bc$recyrs)[bc$censrec == 1] ak <- list(knots = quantile(dtime, c(0, 0.5, 1))) ## ----eval=FALSE------------------------------------------- # sp5 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data = bc, aux = ak, # inits = c(0, 0, 0, 0, 0), # dist = custom.hsurvspline.lh3, # method = "L-BFGS-B", lower = c(-Inf, -Inf, -0.5), # upper = c(Inf, Inf, 0.5), # control = list(trace = 1, REPORT = 1))