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

pbc3$logbili <- log(pbc3$bili)
pbc3$event <- ifelse(pbc3$status != 0, 1, 0)

## we select 4 time points corresponding to quintiles of observed event times
times <- quantile(pbc3$followup[pbc3$event == 1], c(.2,.4,.6,.8))  

## bili values in which to evaluate the smoother
newbili <- seq(min(pbc3$bili), max(pbc3$bili),length.out = 500)
newlogbili <- seq(min(pbc3$logbili), max(pbc3$logbili),length.out = 500)

## scatter plot: pseudo observations vs bilirubin

par(mfrow=c(2,2))
## common range for all 4 plots
range <- range(sapply(times, function(t) pseudosurv(pbc3$followup, pbc3$event, tmax = t)[,3]))

for (t in times){
  pbc3$pseudo <- pseudosurv(pbc3$followup, pbc3$event, tmax = t)[,3]
  plot(pbc3$bili, pbc3$pseudo,
       ylim = range,
       xlab = "Bilirubin",
       ylab = paste("Pseudo-observation, time=", round(t,2)),
       col = "grey"
       )
  smoothed <- predict(loess(pseudo ~ bili, data = pbc3, degree = 1),
                      data.frame(bili = newbili))
  lines(smoothed ~ newbili)
}