Additive Hazard Model

Reading Data into R and importing packages:

library(survival)
library(flexsurv)

## Warning: package 'flexsurv' was built under R version 4.0.3

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

data1=read.csv("C:/Users/Prokarso/Documents/Book3.csv")
data1=data1[,1:5]
data1

Fitting Additive Hazard Model:

data2=data1
data2$Survival.Times=data2$Survival.Times+runif(dim(data2)[1],0,1)*1e-8
attach(data2)
S1=Surv(time=Survival.Times,event=status)
library(addhazard)

## Warning: package 'addhazard' was built under R version 4.0.4

A1=ah(S1~factor(gender)+log.WBC+factor(Rx),data=data2,ties = FALSE)
summary(A1)

## Call:
## ah(formula = S1 ~ factor(gender) + log.WBC + factor(Rx), data = data2, 
##     ties = FALSE)
## 
##                    coef       se lower.95 upper.95     z  p.value    
## factor(gender)  0.02509  0.03023 -0.03415  0.08434 0.830 0.406478    
## log.WBC         0.08286  0.02461  0.03463  0.13109 3.368 0.000759 ***
## factor(Rx)      0.07925  0.02944  0.02155  0.13696 2.692 0.007100 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeff=A1$coef
coeff

## factor(gender)1         log.WBC     factor(Rx)1 
##        0.025093        0.082859        0.079254

h=1:42
for(i in 1:42){
  h[i]=coeff[1]*data2$gender[i]+coeff[2]*data2$log.WBC[i]+coeff[3]*data2$Rx[i]
}

d3=cbind(data2$Survival.Times,h)
d3=d3[which(data2$Rx==1 & data2$status==1),]
colnames(d3)=c("Survival.Times","Hazard Rate")
d3

##       Survival.Times Hazard Rate
##  [1,]             23     0.26758
##  [2,]             22     0.30546
##  [3,]             17     0.32369
##  [4,]             15     0.26983
##  [5,]             12     0.20354
##  [6,]             12     0.33280
##  [7,]             11     0.36843
##  [8,]             11     0.25492
##  [9,]              8     0.37092
## [10,]              8     0.33197
## [11,]              8     0.27149
## [12,]              8     0.37447
## [13,]              5     0.39353
## [14,]              5     0.40821
## [15,]              4     0.46561
## [16,]              4     0.30487
## [17,]              3     0.43661
## [18,]              2     0.51119
## [19,]              2     0.47556
## [20,]              1     0.33635
## [21,]              1     0.51864

d4=d3[order(d3[,1]),]
plot(d4,type = "l",main="Hazard Plot for Placebo")

ss=smooth.spline(d4,spar = 0.001)
plot(ss,type = "l",,main="Hazard Plot for Placebo")

d3=cbind(data2$Survival.Times,h)
d3=d3[which(data2$Rx==0 & data2$status==1),]
colnames(d3)=c("Survival.Times","Hazard Rate")
d3

##       Survival.Times Hazard Rate
##  [1,]             23     0.23804
##  [2,]             22     0.21733
##  [3,]             16     0.32339
##  [4,]             13     0.23863
##  [5,]             10     0.24526
##  [6,]              7     0.36707
##  [7,]              6     0.19140
##  [8,]              6     0.36150
##  [9,]              6     0.27178

d4=d3[order(d3[,1]),]
plot(d4,type = "l",main="Hazard Plot for Treatment")