[R] Kaplan Meier analysis: 95% CI wider in R than in SAS

[Paul Miller]

Hi Enrico,

Not sure how SAS builds the CI but I can look into it. The SAS documentation does have a section on computational formulas at:

http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_lifetest_a0000000259.htm

Although I can't provide my dataset, I can provide the data and code below. This is the R-equivalent of an analysis from "Common Statistical Methods for Clinical Research with SAS Examples."

R produces the follwoing output:

> print(surv.by.vac)
Call: survfit(formula = Surv(WKS, CENS == 0) ~ VAC, data = hsv)

        records n.max n.start events median 0.95LCL 0.95UCL
VAC=GD2      25    25      25     14     35      15      NA
VAC=PBO      23    23      23     17     15      12      35

SAS has the same 95% CI for VAC=GD2 but has a 95% CI of [10, 27] for VAC=PBO. This is just like in the analysis I'm doing currently.

Thanks,

Paul

 
#######################################
#### Chapter 21: The Log-Rank Test ####
#######################################
 
#####################################################
#### Example 21.1: HSV2 Vaccine with gD2 Vaccine ####
#####################################################
 
connection <- textConnection("
GD2  1   8 12  GD2  3 -12 10  GD2  6 -52  7
GD2  7  28 10  GD2  8  44  6  GD2 10  14  8
GD2 12   3  8  GD2 14 -52  9  GD2 15  35 11
GD2 18   6 13  GD2 20  12  7  GD2 23  -7 13
GD2 24 -52  9  GD2 26 -52 12  GD2 28  36 13
GD2 31 -52  8  GD2 33   9 10  GD2 34 -11 16
GD2 36 -52  6  GD2 39  15 14  GD2 40  13 13
GD2 42  21 13  GD2 44 -24 16  GD2 46 -52 13
GD2 48  28  9  PBO  2  15  9  PBO  4 -44 10
PBO  5  -2 12  PBO  9   8  7  PBO 11  12  7
PBO 13 -52  7  PBO 16  21  7  PBO 17  19 11
PBO 19   6 16  PBO 21  10 16  PBO 22 -15  6
PBO 25   4 15  PBO 27  -9  9  PBO 29  27 10
PBO 30   1 17  PBO 32  12  8  PBO 35  20  8
PBO 37 -32  8  PBO 38  15  8  PBO 41   5 14
PBO 43  35 13  PBO 45  28  9  PBO 47   6 15
")

hsv <- data.frame(scan(connection, list(VAC="", PAT=0, WKS=0, X=0)))
hsv <- transform(hsv,
  CENS = ifelse(WKS < 1, 1, 0),
  WKS  = abs(WKS),
  TRT  = ifelse(VAC=="GD2", 1, 0))

library("survival")
surv.by.vac <- survfit(Surv(WKS,CENS==0)~VAC, data=hsv)

plot(surv.by.vac, 
 main = "The Log-Rank Test \n Example 21.1: HSV-Episodes with gD2 Vaccine",
 ylab = "Survival Distribution Function",
 xlab = "Survival Time in Weeks",
 lty = c(1,2))

legend(0.75,0.19, 
 legend = c("gD2","PBO"), 
 lty = c(1,2), title = "Treatment")

summary(surv.by.vac)
print(surv.by.vac)