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

coxphModel <- coxph(Surv(followup,status!=0) ~
                    I((log(bili)-log(17.1))*ifelse(bili >= 17.1, 1, 0)) + factor(ifelse(bili  >= 17.1, 1, 0)),
                    model = TRUE,
                    data = subset(pbc3,followup>0))

## we want the linear predictor translated such that it is 0 for
## bilirubin values under 17.1
lp <- predict(coxphModel, data.frame(bili = min(pbc3$bili)))

lowbili <- seq(min(pbc3$bili), 17.09, length.out = 200)
highbili <- seq(17.1, max(pbc3$bili), length.out = 200)

plot(0,0, type = "n",
     xlim = c(log(min(pbc3$bili)), log(max(pbc3$bili))),
     xlab = "log(Bilirubin)",
     ylab = "Linear Predictor",
     ylim = c(-1.5,6)
     )
lines(log(highbili), predict(coxphModel, data.frame(bili = highbili))-lp)
lines(log(lowbili), predict(coxphModel, data.frame(bili = lowbili))-lp)

curve(coxphModel$coef[2] + coxphModel$coef[1]*(x-log(17.1)),
      log(min(pbc3$bili)), log(max(pbc3$bili)),
      add=TRUE,
      lty = "43")

## marginal histogram over the bilirubin values
hist <- hist(log(pbc3$bili), breaks=ceiling(max(log(pbc3$bili))*10), plot=FALSE)
translation <- 1.5
biliWeight <- (hist$counts)/50
p <- hist$mids[biliWeight != 0]
segments(p,-translation, p, biliWeight[biliWeight != 0]-translation)