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

## means of log(bilirubin) within quintile groups of bilirubin

logbiliQuintileMean <- tapply(log(pbc3$bili), pbc3$biligroup, mean)

names(logbiliQuintileMean) <- NULL

## factor with means of log(bilirubin) within quintile groups of bilirubin
## as levels

pbc3$logbgm <- as.factor(logbiliQuintileMean[as.numeric(pbc3$biligroup)])

## the cox proportional hazards model 

coxphModel <- coxph(Surv(followup,status!=0)~logbgm,
                    data=subset(pbc3,followup>0))

## estimated log(hazard ratio)'s
## added 0 as baseline level

loghr <- c(0,coxphModel$coefficients) 

plot(logbiliQuintileMean, loghr, 
     ylab = "log(hazard ratio)",
     xlab = "log(bilirubin)"
     )