postscript("l03.ps",paper="letter")
par(mfrow=c(2,3))
larynx<-as.data.frame( scan("larynx.dat",
   what=list(stage=0,time=0,y1=0,y2=0,status=0)))
#**********************************************/
#* Life table data is from                    */
#* http://www.census.gov/ipc/www/idbprint.html*/
#**********************************************/
census<-as.data.frame(scan("census.dat",
   what=list(b=0,e=0,weight=0)))
census$time<-(census$b+census$e)/2
census$status<-rep(1,length(census$time))
#*************************************************************/
#* Calculate a life table (aka actuarial estimate) from      */
#* larynx cancer survival times. SAS will generate life table*/
#* estimator if act or lt is used in place of life below.    */
#* The h in the plot request refers to the hazard function,  */
#* and not the integrated hazard function.                   */
#*************************************************************/
source("l03f.R")
larynx.lifetable<-actuarial(larynx$time,larynx$status,
   interval=2*(0:10))
census.lifetable<-actuarial(census$time,census$status,
   weight=census$weight, interval=c(0,1,5*seq(17)))
plot(larynx.lifetable$time,larynx.lifetable$survival,type="l",
   main="Actuarial Estimation for Larynx Data",ylab="Survival")
plot(larynx.lifetable$time,larynx.lifetable$haz,type="l",
   main="Actuarial Estimation for Larynx Data",ylab="Hazard")
plot(census.lifetable$time,census.lifetable$haz,type="l",
   main="Actuarial Estimation for Census Data",ylab="Hazard")
census.lifetable
#*************************************************************/
#* Here are some testing examples, with CIs superimposed.    */
#*************************************************************/
stages23<-larynx[(larynx$stage==2)|(larynx$stage==3),]
stages14<-larynx[(larynx$stage==1)|(larynx$stage==4),]
stages14$fs<-(stages14$status==0)*(abs(stages14$time-5.8)>1.4)+
   (stages14$status==1)*(abs(stages14$time-2.5)>1.5)
library(survival)
# The following code does the test; it's also done by the
# following macro, and so we don't have to do it here.
# survdiff(Surv(time, status) ~ stage, data=stages23, rho=0)
fun.survci<-function(ds,time,status,type,byv,title){
   print(survdiff(Surv(time, status) ~ byv, data=ds, rho=0))
   plot(survfit(Surv(time,status)~byv,data=ds, conf.type=type),
   conf.int=T,main=title,lty=seq(length(unique(byv))))
}
fun.survci(stages23,time,status,"log",stages23$stage,
   "Larynx Cancer Surv. Stage 2,3 by Stage")
fun.survci(stages14,time,status,"log",stages14$stage,
   "Larynx Cancer Surv. Stage 1,4 by Stage")
fun.survci(stages14,time,status,"log",stages14$fs,
   "Larynx Cancer Surv. Stage 4,4 by Artifical Variable")