library(survival)
?coxph

# Create the simplest test data set 
test1 <- list(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              x=c(0,2,1,1,1,0,0))
test1

# Fit a coxph model
coxph(Surv(time, status) ~ x, data=test1) 
# interpetation: increasing x by 1 increases the log hazard by .461
#                increasing x by 1 multiplies the hazard by exp(.461)=1.59


# Real data on recidivism
# (see http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-cox-regression.pdf 
#  for an explanation)
# Rossi <- read.table("Rossi.txt", header=T)
Rossi <- read.table("http://cran.r-project.org/doc/contrib/Fox-Companion/Rossi.txt", header=T)
Rossi[1:5, 1:10]

fit1 <- coxph(Surv(week, arrest) ~ fin + age + race + wexp + mar + paro + prio,
               data=Rossi)
fit1
summary(fit1)
# the last three lines describe the global test comparing
# H0: beta_j = 0 for all j
# Ha: beta_j != 0 for some j

# estimated survival curve at the mean value of the covariates,
#   including 95% pointwise confidence intervals:
plot(survfit(fit1), ylim=c(0.6,1), xlab="Weeks",
     ylab="Proportion Not Rearrested")


# look at effect of financial aid:
attach(Rossi)
# construct new data-frame
Rossi.fin <- data.frame(fin=c(0,1), age=rep(mean(age),2), race=rep(mean(race),2),
                        wexp=rep(mean(wexp),2), mar=rep(mean(mar),2), paro=rep(mean(paro),2),
                        prio=rep(mean(prio),2))
Rossi.fin

# plot survival curves for fin=0 and fin=1:
plot(survfit(fit1, newdata=Rossi.fin),
     lty=c(1,2), col=c("red","blue"), ylim=c(.6, 1),
     ylab="proportion not rearrested",
     xlab="time (weeks)")
legend("topright", legend=c("fin = 0", "fin = 1"), 
       lty=c(1,2), col=c("red", "blue"))

# include pointwise 95% confidence intervals:
plot(survfit(fit1, newdata=Rossi.fin), conf.int=T,
     lty=c(1,2), col=c("red","blue"), ylim=c(.6, 1),
     ylab="proportion not rearrested",
     xlab="time (weeks)")
legend("topright", legend=c("fin = 0", "fin = 1"), 
       lty=c(1,2), col=c("red", "blue"))