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 8. ####################################################################library("survival")################################# NCCTG lung cancer study ################################### read in the NCCTG lung cancer study## (clustered data by institution)data <-read.table("Data//NCCTG//lung.txt")head(data)## Follow up plotlibrary(tidyverse)library(patchwork)# function to plot follow-up by# institution and sexinst_by_sex_fu_plot <-function(df){ df |>ggplot(aes(y =reorder(id, time), x = time, color =factor(2- sex))) +geom_linerange(aes(xmin =0, xmax = time)) +geom_point(aes(shape =factor(status)), size =2, fill ="white") +geom_vline(xintercept =0, linewidth =1) +facet_grid(inst ~ ., scales ="free", space ="free", switch ="y") +theme_minimal() +scale_x_continuous("Time (months)", limits =c(0, 36), breaks =seq(0, 36, by =12),expand =c(0, 0.25)) +scale_y_discrete("Patients (by institution)") +scale_shape_manual(values =c(23, 19), labels =c("Censoring", "Death")) +scale_color_brewer(palette ="Set1", labels =c("Female", "Male"))+theme(strip.background =element_rect(fill ="gray90", color ="gray90"),axis.text.y =element_blank(),axis.ticks.y =element_blank(),axis.line.y =element_blank(),panel.grid.major.y =element_blank(),legend.title =element_blank() )}p1 <-inst_by_sex_fu_plot(data |>filter(inst <=11))p2 <-inst_by_sex_fu_plot(data |>filter(inst >11))mul_lung_fu <- p1 + p2 +plot_layout(ncol =2, guides ="collect") &theme(legend.position ="top")# ggsave("mul_lung_fu.pdf", mul_lung_fu, width = 8, height = 10)# ggsave("mul_lung_fu.eps", mul_lung_fu, width = 8, height = 10)# Fit a Cox model with institution-specific frailty# to account for correlation within institutionobj <-coxph(Surv(time, status) ~ age+factor(sex) + phec + phkn + ptkn + wl +frailty(inst, distribution="gamma"), data = data)summary(obj)# fit a naive Cox model without institution-specific frailtyobj.naive <-coxph(Surv(time,status)~age+factor(sex)+phec+phkn+ptkn + wl,data=data)summary(obj.naive)##################################################### Prediction of subject-specific survival curves# ################################################# Median agemed_age <-median(data$age)# Median ph.karnomed_phkn <-median(data$phkn,na.rm=T)# Median pat.karnomed_ptkn <-median(data$ptkn,na.rm=T)# Median wt.lossmed_wl <-median(data$wl,na.rm=T)# Extract the regression coefficientsbeta <- obj$coefficients# Extract the (only) baseline functionbase_obj <-basehaz(obj,centered=F)eta <- base_obj$hazardt <- base_obj$time# Figure 8.2 Prediction of survival probabilities for a typical patient # of median age (63 years), with median physician-# and patient-rated Karnofsky scores (each 80), and with median# weighted loss (7 pounds) by sex and ECOG score.# Obtain the covariate profiles.## Femalezf0 <-c(med_age,1,0,med_phkn,med_ptkn,med_wl)zf1 <-c(med_age,1,1,med_phkn,med_ptkn,med_wl) zf2 <-c(med_age,1,2,med_phkn,med_ptkn,med_wl) zf3 <-c(med_age,1,3,med_phkn,med_ptkn,med_wl) ## Malezm0 <-c(med_age,0,0,med_phkn,med_ptkn,med_wl)zm1 <-c(med_age,0,1,med_phkn,med_ptkn,med_wl) zm2 <-c(med_age,0,2,med_phkn,med_ptkn,med_wl) zm3 <-c(med_age,0,3,med_phkn,med_ptkn,med_wl) # Plot the preducted survival curvespar(mfrow=c(1,2))plot(t,exp(-exp(sum(beta*zf0))*eta),type="s",xlim=c(0,35),ylim=c(0,1),frame=F,lty=1,main="Female",xlab="Time (months)",ylab="Survival probabilities",lwd=2,cex.lab=1.3,cex.axis=1.3,cex.main=1.3)lines(t,exp(-exp(sum(beta*zf1))*eta),lty=2,lwd=2)lines(t,exp(-exp(sum(beta*zf2))*eta),lty=3,lwd=2)lines(t,exp(-exp(sum(beta*zf3))*eta),lty=4,lwd=2)legend("topright",lty=1:4,lwd=2,cex=1.2,paste("ECOG",0:3))plot(t,exp(-exp(sum(beta*zm0))*eta),type="s",xlim=c(0,35),ylim=c(0,1),frame=F,lty=1,main="Male",xlab="Time (months)",ylab="Survival probabilities",lwd=2,cex.lab=1.3,cex.axis=1.3,cex.main=1.3)lines(t,exp(-exp(sum(beta*zm1))*eta),lty=2,lwd=2)lines(t,exp(-exp(sum(beta*zm2))*eta),lty=3,lwd=2)lines(t,exp(-exp(sum(beta*zm3))*eta),lty=4,lwd=2)legend("topright",lty=1:4,lwd=2,cex=1.2,paste("ECOG",0:3))################################# Diabetic retinopathy study ################################## read in the datadata <-read.table("Data//Diabetic Retinopathy Study//drs.txt")head(data)# fit a bivariate marginal Cox model# with treatment, diabetic type# risk score, and treatment*type interaction# as covariatesobj <-coxph(Surv(time, status) ~ trt + type + trt * type + risk+cluster(id), data = data)summary(obj)# Table 8.1 Marginal Cox model analysis of the Diabetic Retinopathy Study# output tablecoeff <-summary(obj)$coeff# beta estimatec1 <- coeff[,1]# robust se and p-valuec2 <- coeff[,4]c3 <- coeff[,6]# naive se and p-valuec4 <- coeff[,3]c5 <-1-pchisq((c1/c4)^2,1)#output the tablenoquote(round(cbind(c1,c2,c3,c4,c5),3))# Fig. 8.4 Prediction of vision-retention probabilities # for patients with a median risk# score (10) by treatment for each diabetic type.# Lambda_0(t) and tLt <-basehaz(obj,centered = F)t <- Lt$timeL <- Lt$hazard# betabeta <- coeff[,1]# plot the predicted survival functionspar(mfrow=c(1,2))# Compute the survival function for # adult and juvenile patients in control and treatmentadult.contr <-exp(-exp(sum(beta*c(0,0,10,0)))*L)adult.trt <-exp(-exp(sum(beta*c(1,0,10,0)))*L)juv.contr <-exp(-exp(sum(beta*c(0,1,10,0)))*L)juv.trt <-exp(-exp(sum(beta*c(1,1,10,1)))*L)# Plot the predicted survival curvesplot(t,adult.contr,type="s",xlim=c(0,80),ylim=c(0,1),frame.plot=F,lty=3,main="Adult",xlab="Time (months)",ylab="Vision-retention probabilities",lwd=2, cex.lab=1.2,cex.axis=1.2,cex.main=1.2)lines(t,adult.trt,lty=1,lwd=2)plot(t,juv.contr,type="s",xlim=c(0,80),ylim=c(0,1),frame.plot=F,lty=3,main="Juvenile",xlab="Time (months)",ylab="Vision-retention probabilities",lwd=2,cex.lab=1.2,cex.axis=1.2,cex.main=1.2)lines(t,juv.trt,lty=1,lwd=2)