#Sta678 Handout #5 (continued): Cox proportional hazards model 

library(survival)
# read in pbc data

#################################################
# Testing: continuous variable or dummy variable?
#################################################
fit<-coxph(Surv(time,delta)~ageyr+edema+log(bili)+log(protime)
+log(albumin))
fitd<-coxph(Surv(time,delta)~ageyr+factor(edema)+log(bili)+log(protime)
+log(albumin))
summary(fit)
summary(fitd)

#LRT
> 232-231
[1] 1
# chisq with 1 df. The continuous variable is sufficient

## Or equivalently
> fit$loglik
[1] -866.9573 -751.4697
> fitd$loglik
[1] -866.9573 -751.0135
> 2*(751.4697-751.0135)
[1] 0.9124

########################################################################
# compare score test and logrank test for a single categorical covariate
########################################################################
fit1<-coxph(Surv(time,delta)~factor(edema))  
#without using function factor, edema will be treated as numerical (or continuous)
> summary(fit1)
Likelihood ratio test= 62.1  on 2 df,   p=3.23e-14
Wald test            = 91.8  on 2 df,   p=0
Score (logrank) test = 131  on 2 df,   p=0
#why 2 df here?

fit2<-survdiff(Surv(time,delta)~edema)
            N Observed Expected (O-E)^2/E (O-E)^2/V
edema=0   354      116   145.47      5.97      62.3
edema=0.5  44       26    13.05     12.84      14.0
edema=1    20       19     2.47    110.44     113.1
 Chisq= 131  on 2 degrees of freedom, p= 0 
#why 2 df here?

## small differences may be observed due to different ways of tie handling

########################################################################
# compare score test and trend test for a single continuous covariate
# with dose a_j given by x
########################################################################
fit2<-survdiff(Surv(time,delta)~edema)

# trend test with a_j=edema
num<-sum(c(0,0.5,1)*(fit2$obs-fit2$exp))
den<-sqrt(sum( c(0,0.5,1) %*% fit2$var %*% rbind(0,0.5,1) ))

> (num/den)^2
[1] 102.0246
# chisq with 1 df

fit3<-coxph(Surv(time,delta)~edema)
summary(fit3)
#Score (logrank) test = 102  on 1 df,   p=0