Fake dataset of 1000 students’ scores on a midterm and final exam. See Chapter 6 in Regression and Other Stories.


Load packages

library("rstanarm")

Simulate fake data

n <- 1000
# set the random seed to get reproducible results
# change the seed to experiment with variation due to random noise
set.seed(2243)
true_ability <- rnorm(n, 50, 10)
noise_1 <- rnorm(n, 0, 10)
noise_2 <- rnorm(n, 0, 10)
midterm <- true_ability + noise_1
final <- true_ability + noise_2
exams <- data.frame(midterm, final)

Linear regression

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=exams, refresh = 0)
print(fit_1, digits=2)
stan_glm
 family:       gaussian [identity]
 formula:      final ~ midterm
 observations: 1000
 predictors:   2
------
            Median MAD_SD
(Intercept) 24.82   1.40 
midterm      0.51   0.03 

Auxiliary parameter(s):
      Median MAD_SD
sigma 11.59   0.27 

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

Plot midterm and final exam scores

par(mar=c(3, 3, 2, 1), mgp=c(1.7, .5, 0), tck=-.01)
par(pty="s")
plot(midterm, final, xlab="Midterm exam score", ylab="Final exam score", xlim=c(0,100), ylim=c(0,100), xaxs="i", yaxs="i", xaxt="n", yaxt="n", pch=20, cex=.5)
x <- seq(0,100,20)
axis(1, x)
axis(2, x)
for (i in x){
  abline(h=i, col="gray70", lty=2)
  abline(v=i, col="gray70", lty=2)
}
abline(coef(fit_1))

LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogRmFrZSBtaWR0ZXJtIGFuZCBmaW5hbCBleGFtIgphdXRob3I6ICJBbmRyZXcgR2VsbWFuLCBKZW5uaWZlciBIaWxsLCBBa2kgVmVodGFyaSIKZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IHJlYWRhYmxlCiAgICB0b2M6IHRydWUKICAgIHRvY19kZXB0aDogMgogICAgdG9jX2Zsb2F0OiB0cnVlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCi0tLQpGYWtlIGRhdGFzZXQgb2YgMTAwMCBzdHVkZW50cycgc2NvcmVzIG9uIGEgbWlkdGVybSBhbmQgZmluYWwKZXhhbS4gU2VlIENoYXB0ZXIgNiBpbiBSZWdyZXNzaW9uIGFuZCBPdGhlciBTdG9yaWVzLgoKLS0tLS0tLS0tLS0tLQoKCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQobWVzc2FnZT1GQUxTRSwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGNvbW1lbnQ9TkEpCiMgc3dpdGNoIHRoaXMgdG8gVFJVRSB0byBzYXZlIGZpZ3VyZXMgaW4gc2VwYXJhdGUgZmlsZXMKc2F2ZWZpZ3MgPC0gRkFMU0UKYGBgCgojIyMjIExvYWQgcGFja2FnZXMKCmBgYHtyIH0KbGlicmFyeSgicnN0YW5hcm0iKQpgYGAKCiMjIyMgU2ltdWxhdGUgZmFrZSBkYXRhCgpgYGB7ciB9Cm4gPC0gMTAwMAojIHNldCB0aGUgcmFuZG9tIHNlZWQgdG8gZ2V0IHJlcHJvZHVjaWJsZSByZXN1bHRzCiMgY2hhbmdlIHRoZSBzZWVkIHRvIGV4cGVyaW1lbnQgd2l0aCB2YXJpYXRpb24gZHVlIHRvIHJhbmRvbSBub2lzZQpzZXQuc2VlZCgyMjQzKQp0cnVlX2FiaWxpdHkgPC0gcm5vcm0obiwgNTAsIDEwKQpub2lzZV8xIDwtIHJub3JtKG4sIDAsIDEwKQpub2lzZV8yIDwtIHJub3JtKG4sIDAsIDEwKQptaWR0ZXJtIDwtIHRydWVfYWJpbGl0eSArIG5vaXNlXzEKZmluYWwgPC0gdHJ1ZV9hYmlsaXR5ICsgbm9pc2VfMgpleGFtcyA8LSBkYXRhLmZyYW1lKG1pZHRlcm0sIGZpbmFsKQpgYGAKCiMjIyMgTGluZWFyIHJlZ3Jlc3Npb24KClRoZSBvcHRpb24gYHJlZnJlc2ggPSAwYCBzdXByZXNzZXMgdGhlIGRlZmF1bHQgU3RhbiBzYW1wbGluZwpwcm9ncmVzcyBvdXRwdXQuIFRoaXMgaXMgdXNlZnVsIGZvciBzbWFsbCBkYXRhIHdpdGggZmFzdApjb21wdXRhdGlvbi4gRm9yIG1vcmUgY29tcGxleCBtb2RlbHMgYW5kIGJpZ2dlciBkYXRhLCBpdCBjYW4gYmUKdXNlZnVsIHRvIHNlZSB0aGUgcHJvZ3Jlc3MuCgpgYGB7ciB9CmZpdF8xIDwtIHN0YW5fZ2xtKGZpbmFsIH4gbWlkdGVybSwgZGF0YT1leGFtcywgcmVmcmVzaCA9IDApCnByaW50KGZpdF8xLCBkaWdpdHM9MikKYGBgCgojIyMjIFBsb3QgbWlkdGVybSBhbmQgZmluYWwgZXhhbSBzY29yZXMKCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgcGRmKGhlcmUoIkZha2VNaWR0ZXJtRmluYWwvZmlncyIsIkZha2VNaWR0ZXJtRmluYWwxLnBkZiIpLCBoZWlnaHQ9NCwgd2lkdGg9NCkKYGBgCmBgYHtyIH0KcGFyKG1hcj1jKDMsIDMsIDIsIDEpLCBtZ3A9YygxLjcsIC41LCAwKSwgdGNrPS0uMDEpCnBhcihwdHk9InMiKQpwbG90KG1pZHRlcm0sIGZpbmFsLCB4bGFiPSJNaWR0ZXJtIGV4YW0gc2NvcmUiLCB5bGFiPSJGaW5hbCBleGFtIHNjb3JlIiwgeGxpbT1jKDAsMTAwKSwgeWxpbT1jKDAsMTAwKSwgeGF4cz0iaSIsIHlheHM9ImkiLCB4YXh0PSJuIiwgeWF4dD0ibiIsIHBjaD0yMCwgY2V4PS41KQp4IDwtIHNlcSgwLDEwMCwyMCkKYXhpcygxLCB4KQpheGlzKDIsIHgpCmZvciAoaSBpbiB4KXsKICBhYmxpbmUoaD1pLCBjb2w9ImdyYXk3MCIsIGx0eT0yKQogIGFibGluZSh2PWksIGNvbD0iZ3JheTcwIiwgbHR5PTIpCn0KYWJsaW5lKGNvZWYoZml0XzEpKQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKYGBgCgo=