source("http://192.38.117.59/~linearpredictors/datafiles/readPbc3.R")
library(survival)

model <- coxph(Surv(followup, status != 0) ~ strata(tment)
               + log(bili) + alb + log(alkph) + log(asptr) + age
               )

## Get event times and estimate survival curve
fit <- survfit(model)
event.times <- sort(fit$time)

## survival curve estimated at all event times.
fit.all <- summary(fit,
                   times = event.times,
                   extend = TRUE)
## last argument as we need to extend the survival function beyond
## the last event time for one group.

## integrated intensities plotted against
## integrated intensities in reference level.

plot(-log(fit.all$surv[fit.all$strata=="tment=0"]),
     -log(fit.all$surv[fit.all$strata=="tment=1"]),
     type = "s", ## plot as step function
     xlab = "Cumulative hazard of placebo",
     ylab = "Cumulative hazard of CyA"
     )