Chapter 9 - Recurrent Events

Slides

Lecture slides here. (To convert html to pdf, press E \(\to\) Print \(\to\) Destination: Save to pdf)

Base R Code

Show the code

##################################################################
# This code generates all numerical results in chapter 9.       ##
##################################################################



library(survival)

# read in the cgd dataset in counting process format
cgd <- read.table("Data\\Chronic Granulomatous Disease Study\\cgd_counting.txt")
head(cgd)


# Andersen-Gill model
obj.AG <- coxph(Surv(tstart, tstop, status) ~ treat + sex + age + inherit + steroids +
                 propylac, data = cgd)
summary(obj.AG)

# Frailty model
obj.frail <- coxph(Surv(tstart, tstop, status) ~ treat + sex + age + inherit + steroids +
                    propylac + frailty(id, distribution = "gamma"), data = cgd)

summary(obj.frail)


# proportional mean model (LWYY)
obj.pm <- coxph(Surv(tstart, tstop, status) ~ treat + sex + age + inherit + steroids +
                 propylac + cluster(id), data = cgd)

summary(obj.pm)


# extract the beta's from the three models
coeff.AG <- summary(obj.AG)$coeff
coeff.frail <- summary(obj.frail)$coeff
coeff.pm <- summary(obj.pm)$coeff


#########################################
# Table 9.1. beta, se(beta), and p-vaues
# from the three models
#########################################

## Andersen-Gill model
c1 <- coeff.AG[,1]
c2 <- coeff.AG[,3]
c3 <- coeff.AG[,5]

# Frailty model
c4 <- coeff.frail[1:6,1]
c5 <- coeff.frail[1:6,2]
c6 <- coeff.frail[1:6,6]

# Proportional means model
c7 <- coeff.pm[1:6,1]
c8 <- coeff.pm[1:6,4]
c9 <- coeff.pm[1:6,6]

#print out Table 9.1
noquote(round(cbind(c1,c2,c3,c4,c5,c6,c7,c8,c9),3))


### Figure 9.2 #############################################
# predicted  mean functions by treatment
# for a female/male patient of 12 years old with X-linked 
# inheritance pattern and use of both steroids and 
#  prophylactic antibiotics
############################################################

# get beta
beta <- obj.pm$coeff

# get baseline mean function mu_0(t)
# and t
Lt <- basehaz(obj.pm,centered = F)
t <- Lt$time
mu0 <- Lt$hazard

# covariate vector (besides treatement) for female patient
zf <- c(0,12,1,1,1)
# covariate vector (besides treatement) for male patient
zm <- c(1,12,1,1,1)

# female in treatment and control
mu.f.trt <- exp(sum(c(1,zf)*beta))*mu0
mu.f.contr <- exp(sum(c(0,zf)*beta))*mu0

# male in treatment and control
mu.m.trt <- exp(sum(c(1,zm)*beta))*mu0
mu.m.contr <- exp(sum(c(0,zm)*beta))*mu0

# Plot the figure
par(mfrow=c(1,2))

# for female (left panel)
plot(t/30.5, mu.f.trt, type="s",xlim=c(0, 12), ylim=c(0,6),frame.plot =F,lty=1, main="Female",
     xlab="Time (months)",ylab = "Mean number of infections", lwd=2)
lines(t/30.5,mu.f.contr,lty=3,lwd=2)

#for male (right panel)
plot(t/30.5, mu.m.trt, type="s", xlim=c(0, 12), ylim=c(0,6),frame.plot =F,lty=1,main="Male",
     xlab="Time (months)",ylab = "Mean number of infections",lwd=2)
lines(t/30.5,mu.m.contr,lty=3,lwd=2)