R Textbook Examples
Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence
by Judith D. Singer and John B. Willett
Chapter 10: Describing Discrete-time Event Occurrence Data

The examples on the page were written in R version 2.4.1.  This page requires package survival.

library(survival)

Table 10.1 on page 327.

teachers<-read.table("http://www.ats.ucla.edu/stat/examples/alda/teachers.csv", sep=",", header=T)
ts <- survfit( Surv(t, 1-censor)~ 1, conf.type="none", data=teachers)
h<-ts$n.event/ts$n.risk
nlost<-ts$n.risk-ts$n.event- ts$n.risk[-1]
nlost[12] = ts$n.risk[12]-ts$n.event[12]
tab10.1<-cbind(time=ts$time, risk=ts$n.risk, left=ts$n.event, censored=nlost, hazard=h, survival=ts$surv)
tab10.1

      time risk left censored     hazard  survival
 [1,]    1 3941  456        0 0.11570667 0.8842933
 [2,]    2 3485  384        0 0.11018651 0.7868561
 [3,]    3 3101  359        0 0.11576911 0.6957625
 [4,]    4 2742  295        0 0.10758570 0.6209084
 [5,]    5 2447  218        0 0.08908868 0.5655925
 [6,]    6 2229  184        0 0.08254823 0.5189038
 [7,]    7 2045  123      280 0.06014670 0.4876935
 [8,]    8 1642   79      307 0.04811206 0.4642295
 [9,]    9 1256   53      255 0.04219745 0.4446402
[10,]   10  948   35      265 0.03691983 0.4282242
[11,]   11  648   16      241 0.02469136 0.4176508
[12,]   12  391    5      386 0.01278772 0.4123100

Figure 10.1, page 333.

plot(ts$time, h, type="l", ylab="Estimated hazard probability", xlab="years in teaching", 
     ylim=c(0, .15), xlim=c(0, 13))
text(8, .08, "All teachers", adj=0)

plot(ts$time, ts$surv, type="l", ylab="Estimated Survival Probability", xlab="years in teaching", 
     ylim=c(0, 1.), xlim=c(0, 13))
abline(h=c(.5), lty=2)
abline(v=c(6.6), lty=2)
text(8, .6, "All teachers(6.6 years)", adj=0, cex=.7)

Fig. 10.2, page 340

par(mfrow=c(1, 2), cex=0.7)

# Panel A

cocaine<-read.table("http://www.ats.ucla.edu/stat/examples/alda/cocaine_relapse.csv", sep=",", header=T)
ts <- survfit( Surv(week, 1-censor)~ 1, conf.type="none", data=cocaine)
h<-ts$n.event/ts$n.risk

plot(ts$time, h, type="l", ylab=" ", main="Estimated Hazard Function", xlim=c(0, 13))
plot(ts$time, ts$surv, type="l", ylab=" ", main="Estimated Survival Function", xlim=c(0, 13))
abline(h=c(.5), lty=2)
abline(v=c(7.6), lty=2)

# panel B

firstsex<-read.table("http://www.ats.ucla.edu/stat/examples/alda/firstsex.csv", sep=",", header=T)
ts <- survfit( Surv(time, 1-censor)~ 1, conf.type="none", data=firstsex)
h<-ts$n.event/ts$n.risk

plot(ts$time, h, type="l", ylab=" ", main="Estimated Hazard Function", xlim=c(6, 13))
plot(ts$time, ts$surv, type="l", ylab=" ", main="Estimated Survival Function",  xlim=c(6, 13))
abline(h=c(.5), lty=2)
abline(v=c(10.6), lty=2)

# panel C
rm(list=ls())
suicide<-read.table("http://www.ats.ucla.edu/stat/examples/alda/suicide.csv", sep=",", header=T)
ts <- survfit( Surv(time, 1-censor)~ 1, conf.type="none", data=suicide)
h<-ts$n.event/ts$n.risk

plot(ts$time, h, type="l", ylab=" ", main="Estimated Hazard Function", xlim=c(4, 24))
plot(ts$time, ts$surv, type="l", ylab=" ", main="Estimated Survival Function",  xlim=c(4, 24))
abline(h=c(.5), lty=2)
abline(v=c(14.8), lty=2)

# panel D
congress<-read.table("http://www.ats.ucla.edu/stat/examples/alda/congress.csv", sep=",", header=T)
ts <- survfit( Surv(time, 1-censor)~ 1, conf.type="none", data=congress)
h<-ts$n.event/ts$n.risk

plot(ts$time, h, type="l", ylab=" ", main="Estimated Hazard Function", xlim=c(0,8))
plot(ts$time, ts$surv, type="l", ylab=" ", main="Estimated Survival Function",  xlim=c(0,8))
abline(h=c(.5), lty=2)
abline(v=c(3.5), lty=2)

Table 10.2 on page 349. The middle column for the survival function is based on the formula (10.8) on page 350.

rm(list=ls())
teachers<-read.table("http://www.ats.ucla.edu/stat/examples/alda/teachers.csv", sep=",", header=T)
ts <- summary(survfit( Surv(t, 1-censor)~ 1, conf.type="none", error="greenwood", data=teachers))
h<-ts$n.event/ts$n.risk
se.h<-sqrt(h*(1-h)/ts$n.risk)
s.sqrt<-(ts$std.err)^2/(ts$surv)^2
tab10.2<-cbind(t=ts$time, n=ts$n.risk, hazard = h, se.h, survival=ts$surv, s.sqrt, stderr=ts$std.err)
tab10.2

       t    n     hazard        se.h  survival       s.sqrt      stderr
 [1,]  1 3941 0.11570667 0.005095352 0.8842933 3.320134e-05 0.005095352
 [2,]  2 3485 0.11018651 0.005304108 0.7868561 6.873392e-05 0.006523503
 [3,]  3 3101 0.11576911 0.005745506 0.6957625 1.109546e-04 0.007328813
 [4,]  4 2742 0.10758570 0.005917344 0.6209084 1.549210e-04 0.007728276
 [5,]  5 2447 0.08908868 0.005758804 0.5655925 1.948890e-04 0.007895820
 [6,]  6 2229 0.08254823 0.005828952 0.5189038 2.352549e-04 0.007958957
 [7,]  7 2045 0.06014670 0.005257621 0.4876935 2.665487e-04 0.007962239
 [8,]  8 1642 0.04811206 0.005281208 0.4642295 2.973305e-04 0.008004838
 [9,]  9 1256 0.04219745 0.005672654 0.4446402 3.324074e-04 0.008106700
[10,] 10  948 0.03691983 0.006124306 0.4282242 3.728453e-04 0.008268668
[11,] 11  648 0.02469136 0.006096155 0.4176508 4.119139e-04 0.008476499
[12,] 12  391 0.01278772 0.005682161 0.4123100 4.450427e-04 0.008698106

Figure 10.4, page 353

Demonstrating the difference between the person-oriented data set teachers and the person-period data set teachers_pp.

teachers<-read.table("http://www.ats.ucla.edu/stat/examples/alda/teachers.csv", sep=",", header=T)
teachers[id==20 | id == 126 | id ==129, ]
     id  t censor
19   20  3      0
125 126 12      0
128 129 12      1

rm(list=ls())
teachers.pp<-read.table("http://www.ats.ucla.edu/stat/examples/alda/teachers_pp.csv", sep=",", header=T)
t<-subset(teachers.pp, id %in% c(20, 126, 129))
t
     id period event
95   20      1     0
96   20      2     0
97   20      3     1
740 126      1     0
741 126      2     0
742 126      3     0
743 126      4     0
744 126      5     0
745 126      6     0
746 126      7     0
747 126      8     0
748 126      9     0
749 126     10     0
750 126     11     0
751 126     12     1
760 129      1     0
761 129      2     0
762 129      3     0
763 129      4     0
764 129      5     0
765 129      6     0
766 129      7     0
767 129      8     0
768 129      9     0
769 129     10     0
770 129     11     0
771 129     12     0

Table 10.3, page 355

Cross-tabulation of event indicator (event) and time-period indicator (period) in the person-period data set to yield components of the life table.

rm(list=ls())
teachers.pp<-read.table("http://www.ats.ucla.edu/stat/examples/alda/teachers_pp.csv", sep=",", header=T)

temp<-data.frame(table(period, event))
event0<-subset(temp, event==0)
event1<-subset(temp, event==1)
total<-event0$Freq + event1$Freq
proportion<-event1$Freq/total

tab10.3<-cbind(period=event0$period, event0=event0$Freq, event1=event1$Freq, total, proportion)
tab10.3

      period event0 event1 total proportion
 [1,]      1   3485    456  3941 0.11570667
 [2,]      2   3101    384  3485 0.11018651
 [3,]      3   2742    359  3101 0.11576911
 [4,]      4   2447    295  2742 0.10758570
 [5,]      5   2229    218  2447 0.08908868
 [6,]      6   2045    184  2229 0.08254823
 [7,]      7   1922    123  2045 0.06014670
 [8,]      8   1563     79  1642 0.04811206
 [9,]      9   1203     53  1256 0.04219745
[10,]     10    913     35   948 0.03691983
[11,]     11    632     16   648 0.02469136
[12,]     12    386      5   391 0.01278772