source("http://192.38.117.59/~linearpredictors/datafiles/readFever.R")

model <- glm(death ~ alco, 
             family = binomial,
             data = fever)

p.estimates <- predict(model, type="response")
modeldata <- model$model ## the data.frame actually used in the model

## residuals: Difference in observed death and fitted probability of death
modeldata$residuals <- ((unclass(modeldata$death)-1) - p.estimates)/sqrt(p.estimates*(1-p.estimates))

plot(residuals ~ alco,
     data = modeldata,
     xlab = "Drinks per week",
     ylab = "Residual",
     pch = ifelse(modeldata$death==1,3,1)
     )

smooth.res <- loess(residuals ~ alco, degree=1, span = 0.8, data = modeldata)
smoothed.res <- predict(smooth.res, data.frame(alco = seq(0, 15, .1)))
lines(smoothed.res~seq(0,15,.1))