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

##  means within quintile groups of bilirubin level.
biliQuintileMean <- tapply(pbc3$bili, pbc3$biligroup, mean)

##  Fitting model, obtaining log of Hazard Ratios.
coxphModel <- coxph(Surv(followup, status != 0) ~ biligroup,
                    data = subset(pbc3, followup>0))

## Adding 0 to the estimates as reference level
logHazardRatio <- c(0, coxphModel$coefficients) 

##  The final plot
plot(biliQuintileMean, logHazardRatio, 
     xlab = "Bilirubin",
     ylab = "log(Hazard Ratio)")