library(survival) # Section 4.1 computation y <- c(3,5,7,9,18,12,19,20,21,33) stat <- c(1,1,1,0,1,1,1,1,0,0) group <- c(1,1,1,1,1,2,2,2,2,2) w <- 1:10 x <- group == 1 n <- length(y) mean(w) * sum(x) # mean ( sum(w^2) - n * mean(w)^2) * (sum(x^2) - n * mean(x)^2 ) / (n- 1) # general variance formula 5 * 5 * 11 / 12 # Wilcoxon test variance (z <- ( 16 - mean(w) * sum(x) ) / sqrt(5 * 5 * 11 / 12)) 2 * (1 - pnorm(abs(z))) # Section 4.2 computation w <- c(9,7,5,-3,2,0,-2,-4,-7,-7) x <- group == 1 n <- length(y) mean(w) * sum(x) # mean (var.gen <-( sum(w^2) - n * mean(w)^2) * (sum(x^2) - n * mean(x)^2 ) / (n- 1)) # general variance formula (z <- 18 / sqrt(var.gen)) 2 * (1 - pnorm(abs(z))) U.vec <- rep(NA,10^5) for(i in 1:10^5) U.vec[i] <- sum(w*sample(x)) table(U.vec) t.test(U.vec) var(U.vec) t.test(U.vec^2) # Section 4.4 example 1 computation y <- c(2,7,13,15,20,25,34) group <- c(0,0,1,1,1,0,0) stat <- rep(1,7) a <- summary(survfit(Surv(y,stat))) n <- a$n.risk d <- a$n.event n.1 <- c(3,3,3,2,1,0,0) d.1 <- group s.hat<-a$sur (tbl <- cbind(time=a$time,n,d,n.1,d.1,e.1=n.1*d/n,v.1=n.1*(n-n.1)*d*(n-d)/(n^2*(n-1)),s.hat)) sum(tbl[,6]) sum(tbl[,7]) # Section 4.4 example 2 computation y <- c(3,5,7,9,18,12,19,20,21,33) stat <- c(1,1,1,0,1,1,1,1,0,0) group <- c(1,1,1,1,1,2,2,2,2,2) a <- summary(survfit(Surv(y,stat))) n <- a$n.risk d <- a$n.event n.1 <- c(5,4,3,1,1,0,0) d.1 <- c(1,1,1,0,1,0,0) s.hat<-a$sur tbl <- cbind(time=a$time,n,d,n.1,d.1,e.1=n.1*d/n,v.1=n.1*(n-n.1)*d*(n-d)/(n^2*(n-1)),s.hat) (sum(tbl[,5]-tbl[,6]))^2 / sum(tbl[,7]) # log-rank statistic s <- c(1,tbl[-7,8]) sum(s*(tbl[,5]-tbl[,6]))^2 / sum(s^2*tbl[,7]) # rank statistic survdiff(Surv(y,stat) ~ group,rho=0) survdiff(Surv(y,stat) ~ group,rho=1) # Example 1: non parametric test for aml data plot(survfit(Surv(aml$time, aml$status) ~ aml$x), conf.int = F, main = "K-M estimate for AML data") survdiff(Surv(aml$time, aml$status) ~ aml$x,rho=0) survdiff(Surv(aml$time, aml$status) ~ aml$x,rho=1) survdiff(Surv(aml$time, aml$status) ~ aml$x,rho=-1) # Example 2: Simpson's paradox time <- c(3, 8, 18, 20, 21, 22, 24, 25, 10, 13, 15, 16, 17, 19, 25, 32) status <- c(1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1) group <- c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2) sex <- c(1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2) Surv(time[group==1 & sex==1],status[group==1 & sex==1]) Surv(time[group==2 & sex==1],status[group==2 & sex==1]) Surv(time[group==1 & sex==2],status[group==1 & sex==2]) Surv(time[group==2 & sex==2],status[group==2 & sex==2]) survdiff(Surv(time, status) ~ group) survdiff(Surv(time, status) ~ group,subset=(sex==1)) survdiff(Surv(time, status) ~ group,subset=(sex==2)) survdiff(Surv(time, status) ~ group + strata(sex))