require(survival) # Log rank: GVHD (Slide 13) Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/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)