Regression and Other Stories: Arsenic model residuals

Binned residual plots for a logistic regression model: wells in Bangladesh. See Chapter 14 in Regression and Other Stories.

---

Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
library("loo")
invlogit <- plogis

Load data

wells <- read.csv(root("Arsenic/data","wells.csv"))
head(wells)

  switch arsenic   dist dist100 assoc educ educ4
1      1    2.36 16.826 0.16826     0    0  0.00
2      1    0.71 47.322 0.47322     0    0  0.00
3      0    2.07 20.967 0.20967     0   10  2.50
4      1    1.15 21.486 0.21486     0   12  3.00
5      1    1.10 40.874 0.40874     1   14  3.50
6      1    3.90 69.518 0.69518     1    9  2.25

n <- nrow(wells)

Predict switching

with distance, arsenic, education and intercations

wells$c_dist100 <- wells$dist100 - mean(wells$dist100)
wells$c_arsenic <- wells$arsenic - mean(wells$arsenic)
wells$c_educ4 <- wells$educ4 - mean(wells$educ4)

fit_8 <- stan_glm(switch ~ c_dist100 + c_arsenic + c_educ4 +
                      c_dist100:c_educ4 + c_arsenic:c_educ4,
                  family = binomial(link="logit"), data = wells)
pred8 <- fitted(fit_8)

Error rates

error_rate_null <- mean(round(abs(wells$switch-mean(pred8))))
round(error_rate_null, 2)

[1] 0.42

error_rate <- mean(round(abs(wells$switch-pred8)))
round(error_rate, 2)

[1] 0.38

Residual plot

plot(c(0,1), c(-1,1), xlab="Estimated  Pr (switching)", ylab="Observed - estimated",
     type="n", main="Residual plot", mgp=c(2,.5,0))
abline(0,0, col="gray", lwd=.5)
points(pred8, wells$switch-pred8, pch=20, cex=.2)

Binned residual plots

Function for binning residuals

binned_resids <- function (x, y, nclass=sqrt(length(x))){
  breaks.index <- floor(length(x)*(1:(nclass-1))/nclass)
  breaks <- c (-Inf, sort(x)[breaks.index], Inf)
  output <- NULL
  xbreaks <- NULL
  x.binned <- as.numeric (cut (x, breaks))
  for (i in 1:nclass){
    items <- (1:length(x))[x.binned==i]
    x.range <- range(x[items])
    xbar <- mean(x[items])
    ybar <- mean(y[items])
    n <- length(items)
    sdev <- sd(y[items])
    output <- rbind (output, c(xbar, ybar, n, x.range, 2*sdev/sqrt(n)))
  }
  colnames (output) <- c ("xbar", "ybar", "n", "x.lo", "x.hi", "2se")
  return (list (binned=output, xbreaks=xbreaks))
}

Binned residual plot with respect to predicted probability

br8 <- binned_resids(pred8, wells$switch-pred8, nclass=40)$binned
plot(range(br8[,1]), range(br8[,2],br8[,6],-br8[,6]),
     xlab="Estimated  Pr (switching)", ylab="Average residual",
     type="n", main="Binned residual plot", mgp=c(2,.5,0))
abline(0,0, col="gray", lwd=.5)
lines(br8[,1], br8[,6], col="gray", lwd=.5)
lines(br8[,1], -br8[,6], col="gray", lwd=.5)
points(br8[,1], br8[,2], pch=20, cex=.5)

Binned residual plots with respect to predictors

br <- binned_resids(wells$dist, wells$switch-pred8, nclass=40)$binned
plot(range(br[,1]), range(br[,2],br[,6],-br[,6]),
     xlab="Distance to nearest safe well", ylab="Average residual",
     type="n", main="Binned residual plot", mgp=c(2,.5,0))
abline(0,0, col="gray", lwd=.5)
n_within_bin <- length(wells$switch)/nrow(br)
lines(br[,1], br[,6], col="gray", lwd=.5)
lines(br[,1], -br[,6], col="gray", lwd=.5)
points(br[,1], br[,2], pch=20, cex=.5)

br <- binned_resids(wells$arsenic, wells$switch-pred8, nclass=40)$binned
plot(range(0,br[,1]), range(br[,2],br[,6],-br[,6]),
     xlab="Arsenic level", ylab="Average residual",
     type="n", main="Binned residual plot", mgp=c(2,.5,0))
abline (0,0, col="gray", lwd=.5)
lines (br[,1], br[,6], col="gray", lwd=.5)
lines (br[,1], -br[,6], col="gray", lwd=.5)
points (br[,1], br[,2], pch=20, cex=.5)

Predict switching with distance, log(arsenic), education and intercations

Use non-centered predictors for easier plotting

wells$log_arsenic <- log(wells$arsenic)
fit_8b <- stan_glm(switch ~ dist100 + log_arsenic + educ4 +
                      dist100:educ4 + log_arsenic:educ4,
                   family = binomial(link="logit"), data = wells)

Predict switching with distance, log(arsenic), education and intercations

Use non-centered predictors for easier plotting

wells$log_arsenic <- log(wells$arsenic)
fit_8b <- stan_glm(switch ~ dist100 + log_arsenic + educ4 +
                      dist100:educ4 + log_arsenic:educ4,
                   family = binomial(link="logit"), data = wells)
pred8b <- fitted(fit_8b)

Error rate

error_rate <- mean(round(abs(wells$switch-pred8b)))
round(error_rate, 2)

[1] 0.36

Plots for log model

jitter_binary <- function(a, jitt=.05){
  a + (1-2*a)*runif(length(a),0,jitt)
}
plot(c(0,max(wells$arsenic,na.rm=T)*1.02), c(0,1),
     xlab="Arsenic concentration in well water", ylab="Pr (switching)",
     type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0))
points(wells$arsenic, jitter_binary(wells$switch), pch=20, cex=.1)
curve(invlogit(coef(fit_8b)[1]+coef(fit_8b)[2]*0+coef(fit_8b)[3]*log(x)+coef(fit_8b)[4]*mean(wells$educ4)+coef(fit_8b)[5]*0*mean(wells$educ4)+coef(fit_8b)[6]*log(x)*mean(wells$educ4)), from=.5, lwd=.5, add=T)
curve(invlogit(coef(fit_8b)[1]+coef(fit_8b)[2]*.5+coef(fit_8b)[3]*log(x)+coef(fit_8b)[4]*mean(wells$educ4)+coef(fit_8b)[5]*.5*mean(wells$educ4)+coef(fit_8b)[6]*log(x)*mean(wells$educ4)), from=.5, lwd=.5, add=T)
text(.25, .80, "if dist = 0", adj=0, cex=.8)
text(2, .63, "if dist = 50", adj=0, cex=.8)

br <- binned_resids(wells$arsenic, wells$switch-pred8b, nclass=40)$binned
plot(range(0,br[,1]), range(br[,2],br[,6],-br[,6]),
     xlab="Arsenic level", ylab="Average residual", type="n",
     main="Binned residual plot\nfor model with log (arsenic)", mgp=c(2,.5,0))
abline(0,0, col="gray", lwd=.5)
n.within.bin <- length(wells$switch)/nrow(br)
lines(br[,1], br[,6], col="gray", lwd=.5)
lines(br[,1], -br[,6], col="gray", lwd=.5)
points(br[,1], br[,2], pch=20, cex=.5)