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

library(survival)

biliQuintileMean <- tapply(pbc3$bili,  pbc3$biligroup, mean)

model1 <- coxph(Surv(followup,status!=0)~ log(bili) + bili,
                model = TRUE,
                data = subset(pbc3,followup>0))
model2 <- coxph(Surv(followup,status!=0)~ I(bili^(-2)) + I(bili^(-0.5)) + bili,
                model = TRUE,
                data=subset(pbc3,followup>0))

## predicted value in the mean value of the first Quintile group
lp1 <- predict(model1, data.frame(bili = biliQuintileMean[1]))
lp2 <- predict(model2, data.frame(bili = biliQuintileMean[1]))

## evenly distributed values in the observed range of bilirubin
newbili <- seq(min(pbc3$bili), max(pbc3$bili), length.out = 100)

## Linear predictors, set to be zero in the mean value of the first Quintile group
plot(newbili, predict(model1, data.frame(bili = newbili)) - lp1,
     type = "l",
     ylim=c(-1.5,4.5),
     yaxp = c(-1,5,6),
     lty = "dashed",
     xlab = "Bilirubin",
     ylab = "Linear predictor")

lines(newbili, predict(model2, data.frame(bili = newbili)) - lp2)

hist <- hist(pbc3$bili, breaks = ceiling(max(pbc3$bili)), plot = FALSE)
  
translation <- 1.5
biliWeight <- (hist$counts)/50
p <- hist$mids[biliWeight != 0]
segments(p,-translation, p, biliWeight[biliWeight != 0]-translation)