Regression and Other Stories: Introclass

Plot residuals vs. predicted values, or residuals vs. observed values? See Chapter 11 in Regression and Other Stories.

---

Load packages

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

Load data

introclass <- read.table(root("Introclass/data","gradesW4315.dat"), header=TRUE)
head(introclass)

  hw1 hw2 hw3 hw4 midterm hw5 hw6 hw7 final
1  95  88 100  95      80  96  99   0   103
2   0  74  74   0      53  83  97   0    79
3 100   0 105 100      91  96 100  96   122
4   0  90  76 100      63  91  95   0    78
5 100  96  99 100      91  93 100  92   135
6  90  83  95 100      73  89 100  90   117

Linear regression model

The option refresh = 0 supresses the default Stan sampling progress output. This is useful for small data with fast computation. For more complex models and bigger data, it can be useful to see the progress.

fit_1 <- stan_glm(final ~ midterm, data = introclass, refresh = 0)
print(fit_1)

stan_glm
 family:       gaussian [identity]
 formula:      final ~ midterm
 observations: 52
 predictors:   2
------
            Median MAD_SD
(Intercept) 65.4   17.6  
midterm      0.7    0.2  

Auxiliary parameter(s):
      Median MAD_SD
sigma 14.9    1.5  

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Compute residuals

compute predictions from simulations

sims <- as.matrix(fit_1)
predicted <- colMeans(sims[,1] + sims[,2] %*% t(introclass$midterm))

or with built-in function

predicted <- predict(fit_1)
resid <- introclass$final - predicted

Plot residuals vs predicted

plot(predicted, resid, xlab="predicted value", ylab="residual",
     main="Residuals vs.\ predicted values", mgp=c(1.5,.5,0), pch=20, yaxt="n")
axis(2, seq(-40,40,20), mgp=c(1.5,.5,0))
abline(0, 0, col="gray", lwd=.5)

Plot residuals vs observed

plot(introclass$final, resid, xlab="observed value", ylab="residual", main="Residuals vs.\ observed values", mgp=c(1.5,.5,0), pch=20, yaxt="n")
axis(2, seq(-40,40,20), mgp=c(1.5,.5,0))
abline(0, 0, col="gray", lwd=.5)

Simulated fake data

a <- 65
b <- 0.7
sigma <- 15
n <- nrow(introclass)
introclass$final_fake <- a + b*introclass$midterm + rnorm(n, 0, 15)
fit_fake <- stan_glm(final_fake ~ midterm, data = introclass, refresh = 0)

Compute residuals

compute predictions from simulations

sims <- as.matrix(fit_fake)
predicted_fake <- colMeans(sims[,1] + sims[,2] %*% t(introclass$midterm))

or with built-in function

predicted_fake <- predict(fit_fake)
resid_fake <- introclass$final_fake - predicted_fake

Plot residuals vs predicted

plot(predicted_fake, resid_fake, xlab="predicted value", ylab="residual", main="Fake data:  resids vs.\ predicted", mgp=c(1.5,.5,0), pch=20, yaxt="n")
axis(2, seq(-40,40,20), mgp=c(1.5,.5,0))
abline(0, 0, col="gray", lwd=.5)

Plot residuals vs observed

plot(introclass$final_fake, resid_fake, xlab="observed value", ylab="residual", main="Fake data:  resids vs.\ observed", mgp=c(1.5,.5,0), pch=20, yaxt="n")
axis(2, seq(-40,40,20), mgp=c(1.5,.5,0))
abline(0, 0, col="gray", lwd=.5)