Loading packages and data:
library(survival)
library(survminer)
## Warning: package 'survminer' was built under R version 4.0.3
## Loading required package: ggplot2
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 4.0.3
library(flexsurv)
## Warning: package 'flexsurv' was built under R version 4.0.3
time=c(2.4,4.3,3.2,4.5,5.9,3.6,6.2,7,5,1.9,3.5,4.8)
status=c(1,1,0,0,0,1,1,1,0,1,1,0)
stageDisease=c(1,1,1,1,2,1,2,2,2,1,2,2)
age=c(45,86,51,68,47,70,74,62,66,53,81,63)
cancer = data.frame(time,status,stageDisease,age)
Fitting Weibull Model:
weibullfit=flexsurvreg(formula=Surv(time,status)~1,data=cancer,dist='weibull')
weibullfit
## Call:
## flexsurvreg(formula = Surv(time, status) ~ 1, data = cancer,
## dist = "weibull")
##
## Estimates:
## est L95% U95% se
## shape 2.808 1.533 5.143 0.867
## scale 5.782 4.425 7.555 0.789
##
## N = 12, Events: 7, Censored: 5
## Total time at risk: 52.3
## Log-likelihood = -17.48606, df = 2
## AIC = 38.97213
plot(weibullfit,xlab="Time",ylab="Survival Time S(t)", main="Survival Function for Weibull Distribution")
Fitting a semi-parametric model
coxreg = coxph(Surv(time, status) ~ as.factor(stageDisease) + age, data = cancer)
coxreg
## Call:
## coxph(formula = Surv(time, status) ~ as.factor(stageDisease) +
## age, data = cancer)
##
## coef exp(coef) se(coef) z p
## as.factor(stageDisease)2 -1.919457 0.146687 1.167835 -1.644 0.100
## age -0.005195 0.994819 0.038566 -0.135 0.893
##
## Likelihood ratio test=3.63 on 2 df, p=0.1626
## n= 12, number of events= 7
summary(coxreg)
## Call:
## coxph(formula = Surv(time, status) ~ as.factor(stageDisease) +
## age, data = cancer)
##
## n= 12, number of events= 7
##
## coef exp(coef) se(coef) z Pr(>|z|)
## as.factor(stageDisease)2 -1.919457 0.146687 1.167835 -1.644 0.100
## age -0.005195 0.994819 0.038566 -0.135 0.893
##
## exp(coef) exp(-coef) lower .95 upper .95
## as.factor(stageDisease)2 0.1467 6.817 0.01487 1.447
## age 0.9948 1.005 0.92239 1.073
##
## Concordance= 0.698 (se = 0.142 )
## Likelihood ratio test= 3.63 on 2 df, p=0.2
## Wald test = 2.79 on 2 df, p=0.2
## Score (logrank) test = 3.64 on 2 df, p=0.2
ggsurvplot(survfit(coxreg),palette = "#2E9FDE",ggtheme = theme_minimal(),data=cancer, title=" Survival Function for Semi Parametric Model")
Fitting a stratified model
res.cox = coxph(Surv(time, status) ~ stageDisease + age, data = cancer)
summary(res.cox)
## Call:
## coxph(formula = Surv(time, status) ~ stageDisease + age, data = cancer)
##
## n= 12, number of events= 7
##
## coef exp(coef) se(coef) z Pr(>|z|)
## stageDisease -1.919457 0.146687 1.167835 -1.644 0.100
## age -0.005195 0.994819 0.038566 -0.135 0.893
##
## exp(coef) exp(-coef) lower .95 upper .95
## stageDisease 0.1467 6.817 0.01487 1.447
## age 0.9948 1.005 0.92239 1.073
##
## Concordance= 0.698 (se = 0.142 )
## Likelihood ratio test= 3.63 on 2 df, p=0.2
## Wald test = 2.79 on 2 df, p=0.2
## Score (logrank) test = 3.64 on 2 df, p=0.2
stageDisease_df = with(cancer,data.frame(stageDisease = c(1, 2),age = rep(mean(age), 2)))
stageDisease_df
fit = survfit(res.cox, newdata = stageDisease_df)
ggsurvplot(fit, conf.int = TRUE, legend.labs=c("Stage of disease=1", "Stage of disease=2"),ggtheme = theme_minimal(),data=cancer, title=" Survival Function for Stratified Model")
Obtaining Cox-Snell residuals in semiparametric model
res_martingale=residuals(res.cox,type="martingale")
res_cox_snell=status-res_martingale
res_cox_snell
## 1 2 3 4 5 6 7 8
## 0.3462846 1.1105662 0.3356583 1.2194156 0.1994882 0.8356805 0.6578044 1.7001132
## 9 10 11 12
## 0.1807399 0.1523947 0.0782757 0.1835786
For the Weibull model, from observation from the graph, apparently the fit does not seem very good.
For the Semi parametric model and Stratified model, we have a p value of 0.2>0.05 in each case. Thus the fits are not significant.
*****************************************************