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


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

range <- range(mapply(function(t){pseudoresiduals(t, fit)}, times))

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