#  R code to estimate the survivor  
#  function for the VA lung cancer trial 
#  of 137 male patients with inoperable 
#  lung cancer.  This code is posted as  
#    vakm.sas   

# Variables
#     Treatment: 1=standard, 2=test (chemotherapy)
#     Celltype: 1=squamous, 2=smallcell, 
#               3=adeno, 4=large
#     Survival in days
#     Status: 1=dead, 0=censored
#     Karnofsky score
#     Months from Diagnosis
#     Age in years
#     Prior therapy: 0=no, 10=yes
      

# Enter the data into a data frame.

va <- read.table("c:/stat565/va.dat",
      header=F, col.names=c("rx", "cellt", "time", 
       "status", "karno", "months", "age", "prior_rx")) 

va$prior_rx <- va$prior_rx/10;    
va

library(survival)

vafit <- survfit(Surv(time, status) ~ rx, data = va,
                 conf.int=.95, conf.type="log-log",
                 type="kaplan-meier", se.fit=T)
vafit

#options 
#Specify the type of survival curve estimator. 
# Possible values are 
#    "kaplan-meier", "fleming-harrington" or "fh2" 
# if a formula is given and 
#     "aalen" or "kaplan-meier" 
# if the first argument is a coxph object (only the 
# first two characters are necessary). The default 
# is "aalen" when a coxph object is given, and 
# it is "kaplan-meier" otherwise. 
             
# conf.type: specifies the confidence interval type. 
# Possible values are: 
# "none"     for no confidence intervals
# "plain"    for standard intervals, 
# "curve +- k*se(curve)", where k is determined from conf.int, 
# "log" for intervals based on the cumulative hazard or 
#                      log(survival) (the default), 
#  "log-log" for intervals based on the log hazard or   
#               log(-log(survival)). 
#  The last type will never extend past 0 or 1. 


plot(vafit, lty = 2:3, lwd=4, cex=3, xlab="Time (Days)",
      ylab="Survival probability", main=" ")
legend(400, .8, c("Standard", "Chemotherapy"), lty = 2:3)

# Print out estimates of survivor curve and 
# confidence limits

summary(vafit)