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)

## function for calulating pseudo residuals for a model using pbc3 data
pseudoresiduals <- function(t, fit){
  ## pseudobservations
  p <- as.vector(pseudosurv(pbc3$followup, pbc3$event, tmax = t)[,3])
  ## predictions
  s <- as.vector(summary(fit, times = t)$surv)
  ## pseudo-residuals
  (p-s)/sqrt(s*(1-s))
}

## pseudo residuals vs bilirubin

fit <- survfit(coxph(Surv(followup, event) ~ bili, data = pbc3),
               newdata = data.frame(bili = pbc3$bili),
               type = "kaplan-meier")

par(mfrow = c(2,2))
for (t in times){
  pbc3$pseudores <- NULL
  pbc3$pseudores <- pseudoresiduals(t,fit)
  
  plot(pbc3$pseudores ~ pbc3$bili,
       ylim = c(-6.5,5),
       xlab = "Bilirubin",
       ylab = paste("Pseudo-residual, time=", round(t,2)),
       col = "grey")
  smoothed <- predict(loess(pseudores ~ bili, degree = 1,
                            data = pbc3),
                      data.frame(bili = newbili))
  lines(smoothed ~ newbili)
  abline(h = 0,lty = "21")
}