library(survival)

# Log rank: GVHD (Slide 13)
Data <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/gvhd.txt")
survdiff(Surv(Time, Status) ~ Group, Data)

# Log rank: pbc (Slide 15)
survdiff(Surv(time, status!=0) ~ stage, pbc)

# "By hand" (Slide 16)
lrt <- survdiff(Surv(time, status!=0) ~ stage, pbc)
w <- lrt$obs[1:3] - lrt$exp[1:3]
V <- lrt$var[1:3, 1:3]
t(w) %*% solve(V) %*% w
w <- lrt$obs[2:4] - lrt$exp[2:4]
V <- lrt$var[2:4, 2:4]
t(w) %*% solve(V) %*% w
w <- lrt$obs - lrt$exp
V <- lrt$var
t(w) %*% MASS::ginv(V) %*% w

# Peto-prentice: GVHD (Slide 19)
survdiff(Surv(Time, Status) ~ Group, Data, rho=1)

# Gehan test is same as Wilcoxon in absence of censoring
# Gehan test not implemented in survival package
x <- rexp(40)
g <- rep(c("A", "B"), each=20)
wilcox.test(x~g, exact=FALSE, correct=FALSE)$p.value
npsm::gehan.test(x, rep(1, length(x)), g)$p.value
