rm(list=ls())
# Required packages
library(survival)
library(flexsurv)
# lung cancer data set
data(lung)
?lung
head(lung)
## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0
# Survival model using a Weibull distribution
fit.weibull <- flexsurvreg(formula = Surv(time, status) ~ 1, data = lung,
dist = "weibull")
fit.weibull
## Call:
## flexsurvreg(formula = Surv(time, status) ~ 1, data = lung, dist = "weibull")
##
## Estimates:
## est L95% U95% se
## shape 1.3168 1.1652 1.4882 0.0822
## scale 417.7587 372.0394 469.0963 24.7045
##
## N = 228, Events: 165, Censored: 63
## Total time at risk: 69593
## Log-likelihood = -1153.851, df = 2
## AIC = 2311.702
# Survival model using a Weibull distribution
fit.lnorm <- flexsurvreg(formula = Surv(time, status) ~ 1, data = lung,
dist = "lognormal")
fit.lnorm
## Call:
## flexsurvreg(formula = Surv(time, status) ~ 1, data = lung, dist = "lognormal")
##
## Estimates:
## est L95% U95% se
## meanlog 5.6633 5.5104 5.8162 0.0780
## sdlog 1.0976 0.9828 1.2258 0.0619
##
## N = 228, Events: 165, Censored: 63
## Total time at risk: 69593
## Log-likelihood = -1169.269, df = 2
## AIC = 2342.538
# Survival model using a Generalised Gamma distribution
fit.ggama <- flexsurvreg(formula = Surv(time, status) ~ 1, data = lung,
dist = "gengamma")
fit.ggama
## Call:
## flexsurvreg(formula = Surv(time, status) ~ 1, data = lung, dist = "gengamma")
##
## Estimates:
## est L95% U95% se
## mu 6.0766 5.8918 6.2613 0.0943
## sigma 0.7270 0.5971 0.8851 0.0730
## Q 1.1267 0.6763 1.5771 0.2298
##
## N = 228, Events: 165, Censored: 63
## Total time at risk: 69593
## Log-likelihood = -1153.69, df = 3
## AIC = 2313.38
# Comparison of survival functions
plot.new()
plot.window(xlim = c(0,1000), ylim = c(0,1))
axis(1, cex.axis = 1.5); axis(2, cex.axis = 1.5); box(); title(ylab="Survival", xlab="Time", cex.lab = 1.5)
lines(fit.weibull, col="red", lwd.ci=0, lty.ci=1)
lines(fit.lnorm, col="blue", lwd.ci=0, lty.ci=1)
lines(fit.ggama, col="green", lwd.ci=0, lty.ci=1)
legend("topright", legend = c("weibull", "lnorm", "gengamma"),
lty = 1, col = c("red","blue","green"))
# Comparison of hazard functions
plot.new()
plot.window(xlim = c(0,1000), ylim = c(0,0.007))
axis(1, cex.axis = 1.5); axis(2, cex.axis = 1.5); box(); title(ylab="Hazard", xlab="Time", cex.lab = 1.5)
lines(fit.weibull, col="red", lwd.ci=0, lty.ci=1, type = "hazard")
lines(fit.lnorm, col="blue", lwd.ci=0, lty.ci=1, type = "hazard")
lines(fit.ggama, col="green", lwd.ci=0, lty.ci=1, type = "hazard")
legend("topright", legend = c("weibull", "lnorm", "gengamma"),
lty = 1, col = c("red","blue","green"))
# Comparison of the three models using AIC
AIC(fit.weibull)
## [1] 2311.702
AIC(fit.lnorm)
## [1] 2342.538
AIC(fit.ggama)
## [1] 2313.38
# Comparison of survival functions and KM estimator
plot(fit.weibull, col="red", lwd.ci=0, lty.ci=1, ylab="Survival", xlab="Time", cex.lab = 1.5,cex.axis = 1.5)
lines(fit.lnorm, col="blue", lwd.ci=0, lty.ci=1)
lines(fit.ggama, col="green", lwd.ci=0, lty.ci=1)
legend("topright", legend = c("weibull", "lnorm", "gengamma", "KM"),
lty = 1, col = c("red","blue","green","black"))
