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

## dataset with n points evenly spaced over the range of observed bili
n <- 10000
newbili <- seq(min(pbc3$bili), max(pbc3$bili), length.out = n)

## Cox Proportional Hazards regression model
model <- coxph(Surv(followup, status!=0) ~ I(bili^(-2)) + I(bili^(-1/2)) + tment,
               data = pbc3)

## Scatterplot and the two predicted curves
linpredPlacebo <- predict(model,
                          data.frame(bili = newbili,
                                     tment = factor(rep(0, n))),
                          type = "lp"
                          )

linpredTment <- predict(model,
                        data.frame(bili = newbili,
                                   tment = factor(rep(1, n))),
                        type = "lp"
                        )

plot(linpredPlacebo ~ log(newbili),
     type = "l",
     lty = "32",
     ylim = c(-2,3),
     xlab = "log(bilirubin)",
     ylab = "Linear predictor")
lines(linpredTment ~ log(newbili))