Weighted Logrank Test Statistics --- R Calculations ==================================================== ## Consider the example from Sec.1.4, given in Table 7.2, p.210 ## The Dialysis data from Sec. 1.4 & Table 7.2 in Klein and Moeschberger # can now (in 2021) be obtained as the dataset "kidney" in KMsurv library(survival); library(KMsurv) data(kidney) > Dialys.fram = kidney names(Dialys.fram) = c("time","Infec","Gp") Dialys.fram = Dialys.fram[order(Dialys.fram$time),] ## Construct dth-columns equal to reverse cumulative numbers ### of "deaths" = Infections within groups > Tb7p2 = cbind.data.frame(Dialys.fram, dth1 = rev(cumsum(rev( Dialys.fram$Infec*(2-Dialys.fram$Gp)))), dth2 = rev(cumsum(rev( Dialys.fram$Infec*(Dialys.fram$Gp-1)))), rsk1 = rev(cumsum(rev( 2-Dialys.fram$Gp))), rsk2 = rev(cumsum(rev(Dialys.fram$Gp-1)))) ## Construct dth-columns equal to reverse cumulative numbers of deaths ## Now restrict to unique times > Tb7p2 <- Tb7p2[!duplicated(Tb7p2$time),] Tb7p2$dth1 <- rev(diff(c(0,rev(Tb7p2$dth1)))) Tb7p2$dth2 <- rev(diff(c(0,rev(Tb7p2$dth2)))) Tb7p2 <- Tb7p2[Tb7p2$dth1+Tb7p2$dth2 > 0, c(1,4:7)] ## remove "Infec, "Gp"" columns ## This Table now agrees with the first 5 columns of Table 7.2 ## The remaining calculations for logrank numerator are easy > LRnum = sum( Tb7p2$dth1 - (Tb7p2$dth1+Tb7p2$dth2)*Tb7p2$rsk1/( Tb7p2$rsk1+Tb7p2$rsk2)) ### = 3.9636 > LRvar = { dtot <- Tb7p2$dth1+Tb7p2$dth2 Ytot <- Tb7p2$rsk1+Tb7p2$rsk2 sum(Tb7p2$rsk1*Tb7p2$rsk2*dtot*(Ytot-dtot)/(Ytot^2*(Ytot-1))) } ## = 6.2106 ## These are logrank Z1(tau) and variance in Table 7.3, etc. ## Compare ratio 3.9636^2/6.2106 = 2.530 with : > survdiff(Surv(time, Infec) ~ Gp, data=Dialys.fram)$chisq ## = 2.53 ### Similarly for Gehan-Wilcoxon (using dtot, Ytot from above): > GWnum = sum(Ytot*Tb7p2$dth1 - dtot*Tb7p2$rsk1) GWvar = sum(Tb7p2$rsk1*Tb7p2$rsk2*dtot*(Ytot-dtot)/(Ytot-1)) > c(GWnum, GWvar) [1] -9.00 38861.81 ### as in Table 7.3 ## We cannot get the Gehan-Wilcoxon from survdiff, but we CAN get ## the chi-square from the (Fleming-Harrington form of the) ## Peto-Peto modified Wilcoxon (7th entry-line of Table 7.3): > survdiff(Surv(time, Infec) ~ Gp, data=Dialys.fram, rho=1)$chisq [1] 1.386523 ---------------------------------------------------------------- In each of these cases, when you do the R coding yourself, it is convenient to separate the numerators and denominators: e.g., LRnum/sqrt(LRvar) is the "logrank statistic", and its square is the associated chi-square for testing survival differences.