Sta678 Handout #4: Calculating logrank and trend tests

library(survival)

# BNCT data from Ex. 7.7 (p. 240)
bnct<-matrix(scan('http://www-rohan.sdsu.edu/~jjfan/sta678/bnct.txt'),byrow=T,ncol=3)
dim(bnct)
trt<-bnct[,1]
time<-bnct[,2]
status<-bnct[,3]

# getting logrank test statistic in R
fit<-survdiff(Surv(time,status)~trt)
> fit
       N Observed Expected (O-E)^2/E (O-E)^2/V
trt=1 10       10     2.66    20.300    27.365
trt=2 10        9     7.60     0.258     0.431
trt=3 10        8    16.74     4.566    16.587
 Chisq= 33.4  on 2 degrees of freedom, p= 5.64e-08 
## note that (O-E) for the three treatment groups sum to zero

# verify the chisq test statistic using obs, exp and var of the fit
> fit$obs
[1] 10  9  8
> fit$exp
[1]  2.656485  7.600097 16.743418
> fit$var
           [,1]       [,2]      [,3]
[1,]  1.9706489 -0.9545184 -1.016130
[2,] -0.9545184  4.5471241 -3.592606
[3,] -1.0161305 -3.5926057  4.608736

zj<-(fit$obs-fit$exp)[1:2]
zjt<-matrix(zj,ncol=1)
covinv<-solve(fit$var[1:2,1:2])

> zj %*% covinv %*% zjt
         [,1]
[1,] 33.38033
# chisq test statistic with 2 df's

# trend test
#a_j=j
num<-sum((1:3)*(fit$obs-fit$exp))
den<-sqrt(sum( (1:3) %*% fit$var %*% rbind(1,2,3) ))

> (num/den)^2
[1] 30.0511
# num/den gives a test that is asymptotically normal
# Hence, (num/den)^2 is a chisq test statistic with 1 df

#p-value
> 1-pchisq(30.0511,1)
[1] 4.208098e-08