Survival Analysis Part 2: Cox Model and Applications (Based on the book by Hosmer, Lemeshow, May)

rm(list = ls())
setwd("C:\\ASA")
# import data
actg320 = read.table("C:\\ASA\\actg320.DAT", sep = "", header = F)
names(actg320) = c("id",
                   "time",
                   "censor",
                   "time_d",
                   "censor_d",
                   "tx",
                   "txgrp",
                   "strat2",
                   "sex",
                   "raceth",
                   "ivdrug",
                   "hemophil",
                   "karnof",
                   "cd4",
                   "priorzdv",
                   "age")
library(survival)
tx.coxph = coxph( Surv(time, censor)~tx, data = actg320, method = "breslow")
# Table 3.1
summary(tx.coxph)
v2.coxph = coxph( Surv(time, censor)~tx + age + sex + cd4 + priorzdv,
                  data = actg320, method = "breslow")
# Table 3.2
summary(v2.coxph)
v3.coxph = coxph( Surv(time, censor)~tx + age + cd4,
                  data = actg320, method = "breslow")
# Table 3.3
summary(v3.coxph)
v4.coxph = coxph( Surv(time, censor)~tx + age + cd4,
                  data = actg320, method = "efron")
summary(v4.coxph)
# baseline cumulative hazard function
H = basehaz(v4.coxph, centered = TRUE)
plot(H[,2], H[,1], xlab = "Time", ylab = "cumulative hazard function", type = "l")
###################################################################
###################################################################
gbcs = read.table("C:\\ASA\\gbcs.DAT", sep = "", header = F)
names(gbcs) = c("id",
                "diagdate",
                "recdate",
                "deathdate",
                "age",
                "menopause",
                "hormone",
                "size",
                "grade",
                "nodes",
                "prog_recp",
                "estrg_recp",
                "rectime",
                "censrec",
                "survtime",
                "censdead")
crude = coxph( Surv(rectime, censrec)~hormone, data = gbcs, method = "breslow")
summary(crude)
adjusted = coxph( Surv(rectime, censrec)~hormone + size, data = gbcs, method = "breslow")
summary(adjusted)
interaction = coxph( Surv(rectime, censrec)~hormone + size + hormone*size,
                     data = gbcs, method = "breslow")
summary(interaction)
hormone.new = data.frame(hormone=c(0,1))
model = coxph( Surv(rectime, censrec)~hormone, data = gbcs, method = "breslow")
plot(survfit(model, newdata=hormone.new), xlab="Time", ylab="Survival Probability")
points(survfit(model, newdata=hormone.new), pch=c(1,2))
legend(200, 0.5, c("Hormone Absent", "Hormone Present"), pch=c(1,2), bg="transparent")