Lognormal Baseline
rm(list=ls())
# Required packages
library(survival)
library(flexsurv)
library(eha)
## Warning: package 'eha' was built under R version 3.6.2
##
## Attaching package: 'eha'
## The following objects are masked from 'package:flexsurv':
##
## dgompertz, dllogis, hgompertz, Hgompertz, hllogis, Hllogis, hlnorm,
## Hlnorm, hweibull, Hweibull, pgompertz, pllogis, qgompertz, qllogis,
## rgompertz, rllogis
# 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
# AFT model using a lognormal baseline distribution
fit.aft <- aftreg(formula = Surv(time, status) ~ age + sex, data = lung,
dist = "lognormal")
summary(fit.aft)
## Call:
## aftreg(formula = Surv(time, status) ~ age + sex, data = lung,
## dist = "lognormal")
##
## Covariate W.mean Coef Time-Accn se(Coef) Wald p
## age 61.961 0.023 1.024 0.008 0.005
## sex 1.438 -0.519 0.595 0.155 0.001
##
## Baseline parameters:
## log(scale) 6.408 0.593 0.000
## log(shape) -0.051 0.056 0.359
## Baseline life expectancy:
##
## Events 165
## Total time at risk 69593
## Max. log. likelihood -1158.8
## LR test statistic 21
## Degrees of freedom 2
## Overall p-value 2.70206e-05
# PH model using a lognormal baseline distribution
fit.ph <- phreg(formula = Surv(time, status) ~ age + sex , data = lung,
dist = "lognormal")
summary(fit.ph)
## Call:
## phreg(formula = Surv(time, status) ~ age + sex, data = lung,
## dist = "lognormal")
##
## Covariate W.mean Coef Exp(Coef) se(Coef) Wald p
## (Intercept) 52.591 11.567 0.000
## age 61.961 0.016 1.016 0.009 0.076
## sex 1.438 -0.507 0.602 0.166 0.002
##
## log(scale) 82.418 11.721 0.000
## log(shape) -2.037 0.058 0.000
##
## Events 165
## Total time at risk 69593
## Max. log. likelihood -1147.2
## LR test statistic 13.84
## Degrees of freedom 2
## Overall p-value 0.000988434
# Comparison of the three models using AIC
AIC.aft <- -2*(fit.aft$loglik[2]) + 2*length(fit.aft$coefficients)
AIC.ph <- -2*(fit.ph$loglik[2]) + 2*length(fit.ph$coefficients)
AIC.aft
## [1] 2325.5
AIC.ph
## [1] 2304.35
# Comparison
plot(fit.aft)
plot(fit.ph)
Weibull Baseline
rm(list=ls())
# Required packages
library(survival)
library(flexsurv)
library(eha)
# 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
# AFT model using a lognormal baseline distribution
fit.aft <- aftreg(formula = Surv(time, status) ~ age + sex, data = lung,
dist = "weibull")
# PH model using a lognormal baseline distribution
fit.ph <- phreg(formula = Surv(time, status) ~ age + sex , data = lung,
dist = "weibull")
# Comparison of the three models using AIC
AIC.aft <- -2*(fit.aft$loglik[2]) + 2*length(fit.aft$coefficients)
AIC.ph <- -2*(fit.ph$loglik[2]) + 2*length(fit.ph$coefficients)
AIC.aft
## [1] 2302.109
AIC.ph
## [1] 2302.109
# Comparison
plot(fit.aft)
plot(fit.ph)
