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

## Define log-variables, for convenience
pbc3$logbili <- log(pbc3$bili)
pbc3$logasptr <- log(pbc3$asptr)
pbc3$logalkph <- log(pbc3$alkph)

## Recoding the status variable for use with package "pseudo"
pbc3$event <- ifelse(pbc3$status != 0, 1, 0)

## Model of interest
cox.model <- coxph(Surv(followup, event) ~ tment + logbili + alb + logalkph + logasptr + age,
                   data = pbc3)

## The dataset actually used
used.data <- na.omit(get_all_vars(cox.model, pbc3))

fit <- survfit(cox.model,
               newdata = used.data,
               type = "kaplan-meier")

## function for calulating pseudo residuals for the model at time t

pseudoresiduals <- function(t, fit){
  
  ## pseudobservations
  p <- as.vector(pseudosurv(used.data$followup,
                            used.data$event, tmax = t)[,3])
  ## predictions
  s <- as.vector(summary(fit, times = t)$surv)
  ## pseudo-residuals
  (p-s)/sqrt(s*(1-s))
}
  
## we select 4 time points corresponding to quintiles of observed event times
times <- quantile(pbc3$followup[pbc3$event == 1], c(.2,.4,.6,.8))  

## pseudo residuals vs log of bilirubin

par(mfrow = c(2,2))
range <- range(mapply(function(t){pseudoresiduals(t, fit)}, times))

for (t in times){
  used.data$pseudores <- NULL
  used.data$pseudores <- pseudoresiduals(t,fit)
  plot(used.data$pseudores ~ used.data$alb,
       ylim = range, 
       xlab = "Albumin",
       ylab = paste("Pseudo-residual, time=",round(t,2)),
       col = "grey")
  
    ## New alb values in which to evaluate the smoother
  newalb <- seq(min(used.data$alb),max(used.data$alb),length.out = 500)
  smoothed <- predict(loess(pseudores ~ alb,
                            data = used.data, degree = 1), 
                      data.frame(alb = newalb))
  lines(smoothed ~ newalb)
  abline(h = 0, lty = "21")
}