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=