Survival Analysis Assignment Problem Set 3 Question 1
Prokarso Chattopadhyay Roll 36
Reading Data into R and importing packages:
library(survival)
library(flexsurv)## Warning: package 'flexsurv' was built under R version 4.0.3library(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.3data1=read.csv("C:/Users/Prokarso/Documents/Book3.csv")
data1=data1[,1:5]
data1Fitting Cox Proportional Hazard Model:
attach(data1)
survv=Surv(time=Survival.Times,event=status)
s1=coxph(survv~gender+Rx+log.WBC, ties = "breslow")
summary(s1)## Call:
## coxph(formula = survv ~ gender + Rx + log.WBC, ties = "breslow")
##
## n= 42, number of events= 30
##
## coef exp(coef) se(coef) z Pr(>|z|)
## gender 0.2632 1.3010 0.4494 0.586 0.55817
## Rx 1.3909 4.0184 0.4566 3.046 0.00232 **
## log.WBC 1.5936 4.9215 0.3300 4.829 1.37e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## gender 1.301 0.7686 0.5392 3.139
## Rx 4.018 0.2489 1.6419 9.834
## log.WBC 4.922 0.2032 2.5775 9.397
##
## Concordance= 0.851 (se = 0.041 )
## Likelihood ratio test= 43.75 on 3 df, p=2e-09
## Wald test = 31.74 on 3 df, p=6e-07
## Score (logrank) test = 44.57 on 3 df, p=1e-09H_0=basehaz(s1)$hazard ##obtaining cumulative base hazard function
h_0=H_0
for(i in 2:length(H_0))
{
h_0[i]=H_0[i]-H_0[i-1] ## obtaining hazard function
}
wm=median(log.WBC)Creating Hazard Functions:
h_1=h_0*exp(s1$coefficients["gender"]+(s1$coefficients["log.WBC"]*wm)+s1$coefficients["Rx"])
plot(unique(sort(Survival.Times)),h_1,type = "l",main="Hazard Functions vs Time (Cox)",xlab = "Time", ylab = "Hazard Function h(t)" )
legend(0.1,600,legend=c("Gender = Male, Group= Placebo","Gender = Female, Group= Placebo","Gender = Male, Group= Treatment","Gender = Female, Group= Treatment"),col=c("black","blue","red","purple"),lty=c(1,1,1,1))
h_2=h_0*exp((s1$coefficients["log.WBC"]*wm)+s1$coefficients["Rx"])
lines(unique(sort(Survival.Times)),h_2,type = "l",col="blue")
h_3=h_0*exp(s1$coefficients["gender"]+(s1$coefficients["log.WBC"]*wm))
lines(unique(sort(Survival.Times)),h_3,type = "l",col="red")
h_4=h_0*exp((s1$coefficients["log.WBC"]*wm))
lines(unique(sort(Survival.Times)),h_4,type = "l",col="purple")
Fitting Additive Hazard Model:
library(timereg)## Warning: package 'timereg' was built under R version 4.0.3s2=aareg(survv~factor(gender)+log.WBC+factor(Rx))
s2## Call:
## aareg(formula = survv ~ factor(gender) + log.WBC + factor(Rx))
##
## n= 42
## 16 out of 17 unique event times used
##
## slope coef se(coef) z p
## Intercept -0.524 -0.1080 0.0398 -2.71 0.00664
## factor(gender)1 0.123 0.0197 0.0162 1.21 0.22500
## log.WBC 0.210 0.0441 0.0141 3.12 0.00181
## factor(Rx)1 0.125 0.0320 0.0149 2.15 0.03160
##
## Chisq=16.45 on 3 df, p=0.000918; test weights=aalencoeffmat=s2$coefficient
wm=median(log.WBC)
a0=data1[,c(3,4,5)]
a0["intercept"]=rep(1, each=length(a0$gender))
a0=a0[,c(4,1,2,3)]
a1=c(1,1,wm,1)
h1=coeffmat%*%a1
plot(h1,type="l",xlab="Time",ylab = "Hazard Function h(t)",main = "Hazard Function for Additive Hazard Model",cex.main=1.1,ylim=c(-0.5,1))
legend(0,1,legend=c("Gender=Male,Group=Placebo","Gender=Female,Group=Placebo","Gender=Male,Group=Treatment","Gender=Female,Group=Treatment"),col=c("black","red","blue","purple"),lty=c(1,1,1,1))
a2=c(1,0,wm,1)
h2=coeffmat%*%a2
lines(h2,col="red")
a3=c(1,1,wm,0)
h3=coeffmat%*%a3
lines(h3,col="blue")
a4=c(1,0,wm,0)
h4=coeffmat%*%a4
lines(h4,col="purple")
Creating Parametric Hazard Models:
(i) Using Weibull Distribution:
S1=sort(Survival.Times)
s3=flexsurvreg(survv~factor(gender)+log.WBC+factor(Rx),dist="weibull")
s3## Call:
## flexsurvreg(formula = survv ~ factor(gender) + log.WBC + factor(Rx),
## dist = "weibull")
##
## Estimates:
## data mean est L95% U95% se exp(est)
## shape NA 2.258 1.702 2.995 0.325 NA
## scale NA 235.605 118.508 468.403 82.604 NA
## factor(gender)1 0.476 -0.156 -0.508 0.195 0.179 0.855
## log.WBC 2.930 -0.792 -0.998 -0.586 0.105 0.453
## factor(Rx)1 0.500 -0.722 -1.117 -0.327 0.201 0.486
## L95% U95%
## shape NA NA
## scale NA NA
## factor(gender)1 0.602 1.215
## log.WBC 0.369 0.557
## factor(Rx)1 0.327 0.721
##
## N = 42, Events: 30, Censored: 12
## Total time at risk: 541
## Log-likelihood = -89.69631, df = 5
## AIC = 189.3926alpha=2.258
malpha=(-1)*alpha
lambda=1/(235.605^alpha)
b1=-0.156
b2=-0.792
b3=-0.722
t=sort(Survival.Times)
h0w=alpha*(t^(alpha-1)) #baseline hazard
hw1=h0w*exp((malpha*b1)+(malpha*b2*wm)+(malpha*b3))
hw2=h0w*exp((malpha*b2*wm)+(malpha*b3))
hw3=h0w*exp((malpha*b1)+(malpha*b2*wm))
hw4=h0w*exp((malpha*b2*wm))
plot(S1,hw1,type="l",xlab="Time",ylab="Hazard Function h(t)",main="Hazard Function for Weibull Distribution")
lines(S1,hw2,col="blue")
lines(S1,hw3,col="red")
lines(S1,hw4,col="green")
legend(0,200000,legend=c("Gender=Male,Group=Placebo","Gender=Female,Group=Placebo","Gender=Male,Group=Treatment","Gender=Female,Group=Treatment"),col=c("black","blue","red","purple"),lty=c(1,1,1,1))
(i) Using Log LOgistic Distribution:
logisfit=flexsurvreg (survv~factor(gender)+log.WBC+factor(Rx), dist="llogis")
logisfit ## Call:
## flexsurvreg(formula = survv ~ factor(gender) + log.WBC + factor(Rx),
## dist = "llogis")
##
## Estimates:
## data mean est L95% U95% se exp(est)
## shape NA 2.9726 2.2117 3.9954 0.4485 NA
## scale NA 151.2566 75.9014 301.4248 53.2142 NA
## factor(gender)1 0.4762 -0.0408 -0.4381 0.3566 0.2027 0.9601
## log.WBC 2.9302 -0.7278 -0.9431 -0.5125 0.1098 0.4830
## factor(Rx)1 0.5000 -0.7494 -1.1674 -0.3314 0.2133 0.4726
## L95% U95%
## shape NA NA
## scale NA NA
## factor(gender)1 0.6452 1.4285
## log.WBC 0.3894 0.5990
## factor(Rx)1 0.3112 0.7179
##
## N = 42, Events: 30, Censored: 12
## Total time at risk: 541
## Log-likelihood = -92.31117, df = 5
## AIC = 194.6223Survival.Times1=sort(Survival.Times)
alphallogis=2.9726
### The beta estimates are:
betal=-0.0408 #for gender
beta2=-0.7278 #for log.WBC
beta3=-0.7494 #for group
ma=alphallogis*(-1)
# Hazard for gender=male, group = placebo
l1=exp((ma*betal) + (ma*beta2*wm) + (ma*beta3)) #calculating lambda for male, placebo
n1=alphallogis*l1*(Survival.Times1^(alphallogis-1)) # The numerator
d1= 1+ (l1*(Survival.Times1^(alphallogis))) #the denominator
h1=n1/d1 # the hazard function
h1## [1] 2.97193613 2.97193613 1.48625770 1.48625770 0.99085822 0.74314731
## [7] 0.74314731 0.59451889 0.59451889 0.49543280 0.49543280 0.49543280
## [13] 0.49543280 0.42465685 0.37157483 0.37157483 0.37157483 0.37157483
## [19] 0.33028878 0.29725993 0.29725993 0.27023632 0.27023632 0.27023632
## [25] 0.24771663 0.24771663 0.22866151 0.19817332 0.18578749 0.17485881
## [31] 0.17485881 0.15645263 0.14863000 0.13511818 0.13511818 0.12924348
## [37] 0.12924348 0.11890400 0.09289375 0.09289375 0.08742941 0.08493143plot (Survival.Times1,h1, type='l',xlab='Time', col='red', ylab='Hazard function h(t)',main="Hazard Function for Log Logistic Distribution")
legend(15,3,legend=c("Gender=Male,Group=Placebo","Gender=Female,Group=Placebo","Gender=Male,Group=Treatment","Gender=Female,Group=Treatment"),col=c("black","blue","red","purple"),lty=c(1,1,1,1))
# Hazard for gender=female, group = placebo
l2=exp((ma*beta2*wm)+(ma*beta3)) #calculating lambda for female,placebo
n2=alphallogis*l2*(Survival.Times1^(alphallogis-1)) # The numerator
d2= 1+(l2*(Survival.Times1^(alphallogis))) #the denominator
h2=n2/d2 # the hazard function
h2## [1] 2.97185055 2.97185055 1.48625225 1.48625225 0.99085713 0.74314696
## [7] 0.74314696 0.59451875 0.59451875 0.49543273 0.49543273 0.49543273
## [13] 0.49543273 0.42465681 0.37157481 0.37157481 0.37157481 0.37157481
## [19] 0.33028877 0.29725992 0.29725992 0.27023631 0.27023631 0.27023631
## [25] 0.24771663 0.24771663 0.22866151 0.19817332 0.18578749 0.17485881
## [31] 0.17485881 0.15645263 0.14862999 0.13511818 0.13511818 0.12924348
## [37] 0.12924348 0.11890400 0.09289375 0.09289375 0.08742941 0.08493143lines(Survival.Times1,h2,col='blue')
# Hazard for gender=male, group = treatment
l3=exp((ma*betal) + (ma*beta2*wm)) #calculating lambda for male,treatment
n3=alphallogis*l3*(Survival.Times1^ (alphallogis-1)) # The numerator
d3= 1+ (l3*(Survival.Times1^(alphallogis))) #the denominator
h3=n3/d3 # the hazard function
lines(Survival.Times1,h3, type='l',col='red')
# Hazard for gender=female, group = treatment
l4=exp((ma*beta2*wm)) #calculating lambda for male,treatment
n4=alphallogis*l4*(Survival.Times1^ (alphallogis-1)) # The numerator
d4= 1+(l4*(Survival.Times1^(alphallogis))) #the denominator
h4=n4/d4 # the hazard function
lines(Survival.Times1,h4, type='l',col='purple')
Comments: The Cox Model has an irregular hazard function, which is neither increasing nor decreasing, and the four groups have more or less distinctive hazard functions, although overlapping at certain points.