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 log of bilirubin

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

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