Trend of record times in the mile run. See Chapter 3 in Regression and Other Stories.


Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
library("ggplot2")
theme_set(bayesplot::theme_default(base_family = "sans"))

Load data

mile <- read.csv(root("Mile/data","mile.csv"), header=TRUE)
head(mile)
    yr month min  sec     year seconds
1 1913     5   4 14.4 1913.417   254.4
2 1915     7   4 12.6 1915.583   252.6
3 1923     8   4 10.4 1923.667   250.4
4 1931    10   4  9.2 1931.833   249.2
5 1933     7   4  7.6 1933.583   247.6
6 1934     6   4  6.8 1934.500   246.8

Linear 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 <- stan_glm(seconds ~ year, data = mile, refresh = 0)
print(fit, digits=2)
stan_glm
 family:       gaussian [identity]
 formula:      seconds ~ year
 observations: 32
 predictors:   2
------
            Median  MAD_SD 
(Intercept) 1006.38   22.65
year          -0.39    0.01

Auxiliary parameter(s):
      Median MAD_SD
sigma 1.42   0.18  

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

Predictions for 1900 and 2000

print(1007 -.393*c(1900,2000))  # Approx
[1] 260.3 221.0
print(coef(fit)[1] + coef(fit)[2]*c(1900,2000), digits=4) # Exact
[1] 260.0 220.8

Example of increasing trend

a <- 0.15
b <- 0.4
par(mar=c(3,3,1,1), mgp=c(2,.5,0), tck=-.01)
plot(c(0,2.2), c(0,a+2.2*b), pch=20, cex=.5, main="y = a + bx (with b > 0)",
  bty="l", type="n", xlab="x", ylab="y", xaxt="n", yaxt="n", xaxs="i", yaxs="i")
axis(1, c(0,1,2))
axis(2, c(a,a+b,a+2*b), c("a","a+b","a+2b"))
abline(a, b)

Example of decreasing trend

a <- 0.95
b <- -0.4
par(mar=c(3,3,1,1), mgp=c(2,.5,0), tck=-.01)
plot(c(0,2.2), c(0,a+.2), pch=20, cex=.5, main="y = a + bx (with b < 0)",
  bty="l", type="n", xlab="x", ylab="y", xaxt="n", yaxt="n", xaxs="i", yaxs="i")
axis(1, c(0,1,2))
axis(2, c(a,a+b,a+2*b), c("a","a+b","a+2b"))
abline(a, b)

Approximate trend from the fit in range [0,2.1]

par(mar=c(3,3,1,1), mgp=c(2,.5,0), tck=-.01)
curve(1007 - 0.393*x, from=0, to=2.1, xlab="x", ylab="y", bty="l", xaxs="i",
  main="y = 1007 - 0.393x")

Approximate trend from the fit in range [0,100]

par(mar=c(3,3,1,1), mgp=c(2,.5,0), tck=-.01)
curve(1007 - 0.393*x, from=0, to=100, xlab="x", ylab="y", bty="l", xaxs="i",
  main="y = 1007 - 0.393x")

Approximate trend of record times in the mile run from 1900 to 2000

par(mar=c(3,3,1,1), mgp=c(2,.5,0), tck=-.01)
curve(1007 - 0.393*x, from=1900, to=2000,
      xlab="Year", ylab="Time (seconds)", bty="l", xaxs="i",
      main="Approx. trend of record times in the mile run",
      ylim=c(215, 265))
text(1960, 246, "y = 1007 - 0.393x")

ggplot version

ggplot(aes(x=year, y=seconds), data=mile) + geom_point(shape=1, size=2) +
    geom_abline(intercept=fit$coefficients[1], slope=fit$coefficients[2]) +
    labs(x="Year", y="Time (seconds)",
         title = "Approx. trend of record times in the mile run")

LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogTWlsZSIKYXV0aG9yOiAiQW5kcmV3IEdlbG1hbiwgSmVubmlmZXIgSGlsbCwgQWtpIFZlaHRhcmkiCmRhdGU6ICJgciBmb3JtYXQoU3lzLkRhdGUoKSlgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIHRoZW1lOiByZWFkYWJsZQogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDIKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQotLS0KVHJlbmQgb2YgcmVjb3JkIHRpbWVzIGluIHRoZSBtaWxlIHJ1bi4gU2VlIENoYXB0ZXIgMyBpbiBSZWdyZXNzaW9uCmFuZCBPdGhlciBTdG9yaWVzLgoKLS0tLS0tLS0tLS0tLQoKCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQobWVzc2FnZT1GQUxTRSwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGNvbW1lbnQ9TkEpCiMgc3dpdGNoIHRoaXMgdG8gVFJVRSB0byBzYXZlIGZpZ3VyZXMgaW4gc2VwYXJhdGUgZmlsZXMKc2F2ZWZpZ3MgPC0gRkFMU0UKYGBgCgojIyMjIExvYWQgcGFja2FnZXMKCmBgYHtyIH0KbGlicmFyeSgicnByb2pyb290IikKcm9vdDwtaGFzX2ZpbGUoIi5ST1MtRXhhbXBsZXMtcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpsaWJyYXJ5KCJyc3RhbmFybSIpCmxpYnJhcnkoImdncGxvdDIiKQp0aGVtZV9zZXQoYmF5ZXNwbG90Ojp0aGVtZV9kZWZhdWx0KGJhc2VfZmFtaWx5ID0gInNhbnMiKSkKYGBgCgojIyMjIExvYWQgZGF0YQoKYGBge3IgfQptaWxlIDwtIHJlYWQuY3N2KHJvb3QoIk1pbGUvZGF0YSIsIm1pbGUuY3N2IiksIGhlYWRlcj1UUlVFKQpoZWFkKG1pbGUpCmBgYAoKIyMjIyBMaW5lYXIgbW9kZWwKClRoZSBvcHRpb24gYHJlZnJlc2ggPSAwYCBzdXByZXNzZXMgdGhlIGRlZmF1bHQgU3RhbiBzYW1wbGluZwpwcm9ncmVzcyBvdXRwdXQuIFRoaXMgaXMgdXNlZnVsIGZvciBzbWFsbCBkYXRhIHdpdGggZmFzdApjb21wdXRhdGlvbi4gRm9yIG1vcmUgY29tcGxleCBtb2RlbHMgYW5kIGJpZ2dlciBkYXRhLCBpdCBjYW4gYmUKdXNlZnVsIHRvIHNlZSB0aGUgcHJvZ3Jlc3MuCgpgYGB7ciB9CmZpdCA8LSBzdGFuX2dsbShzZWNvbmRzIH4geWVhciwgZGF0YSA9IG1pbGUsIHJlZnJlc2ggPSAwKQpwcmludChmaXQsIGRpZ2l0cz0yKQpgYGAKCiMjIyMgUHJlZGljdGlvbnMgZm9yIDE5MDAgYW5kIDIwMDAKCmBgYHtyIH0KcHJpbnQoMTAwNyAtLjM5MypjKDE5MDAsMjAwMCkpICAjIEFwcHJveApwcmludChjb2VmKGZpdClbMV0gKyBjb2VmKGZpdClbMl0qYygxOTAwLDIwMDApLCBkaWdpdHM9NCkgIyBFeGFjdApgYGAKCiMjIyMgRXhhbXBsZSBvZiBpbmNyZWFzaW5nIHRyZW5kCgpgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJNaWxlL2ZpZ3MiLCJhcGx1c2J4MWEucGRmIiksIGhlaWdodD0zLjUsIHdpZHRoPTUpCmBgYApgYGB7ciB9CmEgPC0gMC4xNQpiIDwtIDAuNApwYXIobWFyPWMoMywzLDEsMSksIG1ncD1jKDIsLjUsMCksIHRjaz0tLjAxKQpwbG90KGMoMCwyLjIpLCBjKDAsYSsyLjIqYiksIHBjaD0yMCwgY2V4PS41LCBtYWluPSJ5ID0gYSArIGJ4ICh3aXRoIGIgPiAwKSIsCiAgYnR5PSJsIiwgdHlwZT0ibiIsIHhsYWI9IngiLCB5bGFiPSJ5IiwgeGF4dD0ibiIsIHlheHQ9Im4iLCB4YXhzPSJpIiwgeWF4cz0iaSIpCmF4aXMoMSwgYygwLDEsMikpCmF4aXMoMiwgYyhhLGErYixhKzIqYiksIGMoImEiLCJhK2IiLCJhKzJiIikpCmFibGluZShhLCBiKQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKYGBgCgojIyMjIEV4YW1wbGUgb2YgZGVjcmVhc2luZyB0cmVuZAoKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBwZGYocm9vdCgiTWlsZS9maWdzIiwiYXBsdXNieDFiLnBkZiIpLCBoZWlnaHQ9My41LCB3aWR0aD01KQpgYGAKYGBge3IgfQphIDwtIDAuOTUKYiA8LSAtMC40CnBhcihtYXI9YygzLDMsMSwxKSwgbWdwPWMoMiwuNSwwKSwgdGNrPS0uMDEpCnBsb3QoYygwLDIuMiksIGMoMCxhKy4yKSwgcGNoPTIwLCBjZXg9LjUsIG1haW49InkgPSBhICsgYnggKHdpdGggYiA8IDApIiwKICBidHk9ImwiLCB0eXBlPSJuIiwgeGxhYj0ieCIsIHlsYWI9InkiLCB4YXh0PSJuIiwgeWF4dD0ibiIsIHhheHM9ImkiLCB5YXhzPSJpIikKYXhpcygxLCBjKDAsMSwyKSkKYXhpcygyLCBjKGEsYStiLGErMipiKSwgYygiYSIsImErYiIsImErMmIiKSkKYWJsaW5lKGEsIGIpCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIGRldi5vZmYoKQpgYGAKCiMjIyMgQXBwcm94aW1hdGUgdHJlbmQgZnJvbSB0aGUgZml0IGluIHJhbmdlIFswLDIuMV0KCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgcGRmKHJvb3QoIk1pbGUvZmlncyIsImFwbHVzYngyYS5wZGYiKSwgaGVpZ2h0PTMuNSwgd2lkdGg9NSkKYGBgCmBgYHtyIH0KcGFyKG1hcj1jKDMsMywxLDEpLCBtZ3A9YygyLC41LDApLCB0Y2s9LS4wMSkKY3VydmUoMTAwNyAtIDAuMzkzKngsIGZyb209MCwgdG89Mi4xLCB4bGFiPSJ4IiwgeWxhYj0ieSIsIGJ0eT0ibCIsIHhheHM9ImkiLAogIG1haW49InkgPSAxMDA3IC0gMC4zOTN4IikKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgZGV2Lm9mZigpCmBgYAoKIyMjIyBBcHByb3hpbWF0ZSB0cmVuZCBmcm9tIHRoZSBmaXQgaW4gcmFuZ2UgWzAsMTAwXQoKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBwZGYocm9vdCgiTWlsZS9maWdzIiwiYXBsdXNieDJiLnBkZiIpLCBoZWlnaHQ9My41LCB3aWR0aD01KQpgYGAKYGBge3IgfQpwYXIobWFyPWMoMywzLDEsMSksIG1ncD1jKDIsLjUsMCksIHRjaz0tLjAxKQpjdXJ2ZSgxMDA3IC0gMC4zOTMqeCwgZnJvbT0wLCB0bz0xMDAsIHhsYWI9IngiLCB5bGFiPSJ5IiwgYnR5PSJsIiwgeGF4cz0iaSIsCiAgbWFpbj0ieSA9IDEwMDcgLSAwLjM5M3giKQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkgICAgIApgYGAKCiMjIyMgQXBwcm94aW1hdGUgdHJlbmQgb2YgcmVjb3JkIHRpbWVzIGluIHRoZSBtaWxlIHJ1biBmcm9tIDE5MDAgdG8gMjAwMAoKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBwZGYocm9vdCgiTWlsZS9maWdzIiwiYXBsdXNieDMucGRmIiksIGhlaWdodD0zLjUsIHdpZHRoPTUpCmBgYApgYGB7ciB9CnBhcihtYXI9YygzLDMsMSwxKSwgbWdwPWMoMiwuNSwwKSwgdGNrPS0uMDEpCmN1cnZlKDEwMDcgLSAwLjM5Myp4LCBmcm9tPTE5MDAsIHRvPTIwMDAsCiAgICAgIHhsYWI9IlllYXIiLCB5bGFiPSJUaW1lIChzZWNvbmRzKSIsIGJ0eT0ibCIsIHhheHM9ImkiLAogICAgICBtYWluPSJBcHByb3guIHRyZW5kIG9mIHJlY29yZCB0aW1lcyBpbiB0aGUgbWlsZSBydW4iLAogICAgICB5bGltPWMoMjE1LCAyNjUpKQp0ZXh0KDE5NjAsIDI0NiwgInkgPSAxMDA3IC0gMC4zOTN4IikKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgZGV2Lm9mZigpCmBgYAoKIyMjIyBnZ3Bsb3QgdmVyc2lvbgoKYGBge3IgfQpnZ3Bsb3QoYWVzKHg9eWVhciwgeT1zZWNvbmRzKSwgZGF0YT1taWxlKSArIGdlb21fcG9pbnQoc2hhcGU9MSwgc2l6ZT0yKSArCiAgICBnZW9tX2FibGluZShpbnRlcmNlcHQ9Zml0JGNvZWZmaWNpZW50c1sxXSwgc2xvcGU9Zml0JGNvZWZmaWNpZW50c1syXSkgKwogICAgbGFicyh4PSJZZWFyIiwgeT0iVGltZSAoc2Vjb25kcykiLAogICAgICAgICB0aXRsZSA9ICJBcHByb3guIHRyZW5kIG9mIHJlY29yZCB0aW1lcyBpbiB0aGUgbWlsZSBydW4iKQoKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CiMgQSBzaW1wbGUgZ3JhcGgKIyBXb3JsZCByZWNvcmQgdGltZXMgaW4gdGhlIG1pbGUgcnVuIGZyb20gMTkwMCB0byAyMDAwCmlmIChzYXZlZmlncykgcGRmKHJvb3QoIk1pbGUvZmlncyIsIm1pbGUxYS5wZGYiKSwgaGVpZ2h0PTMuNSwgd2lkdGg9NSkKcGxvdChtaWxlJHllYXIsIG1pbGUkc2Vjb25kcywKICAgICBtYWluPSJXb3JsZCByZWNvcmQgdGltZXMgaW4gdGhlIG1pbGUgcnVuIikKaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQojIFdvcmxkIHJlY29yZCB0aW1lcyBpbiB0aGUgbWlsZSBydW4gZnJvbSAxOTAwIHRvIDIwMDAKIyBJbXByb3ZlZCBncmFwaAppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJNaWxlL2ZpZ3MiLCJtaWxlMWIucGRmIiksIGhlaWdodD0zLjUsIHdpZHRoPTUpCnBhcihtYXI9YygzLDMsMywxKSwgbWdwPWMoMiwuNSwwKSwgdGNrPS0uMDEpCnBsb3QobWlsZSR5ZWFyLCBtaWxlJHNlY29uZHMsIGJ0eT0ibCIsCiAgICAgbWFpbj0iV29ybGQgcmVjb3JkIHRpbWVzIGluIHRoZSBtaWxlIHJ1biIsIHhsYWI9IlllYXIiLCB5bGFiPSJTZWNvbmRzIikKaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQojIFdvcmxkIHJlY29yZCB0aW1lcyBpbiB0aGUgbWlsZSBydW4gZnJvbSAxOTAwIHRvIDIwMDAKIyBGaXR0ZWQgbGluZQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJNaWxlL2ZpZ3MiLCJtaWxlMi5wZGYiKSwgaGVpZ2h0PTMuNSwgd2lkdGg9NSkKcGFyKG1hcj1jKDMsMywzLDEpLCBtZ3A9YygyLC41LDApLCB0Y2s9LS4wMSkKcGxvdChtaWxlJHllYXIsIG1pbGUkc2Vjb25kcywgYnR5PSJsIiwKICAgICBtYWluPSJXb3JsZCByZWNvcmQgdGltZXMgaW4gdGhlIG1pbGUgcnVuIiwgeGxhYj0iWWVhciIsIHlsYWI9IlNlY29uZHMiKQpjdXJ2ZShjb2VmKGZpdClbMV0gKyBjb2VmKGZpdClbMl0qeCwgYWRkPVRSVUUpCmlmIChzYXZlZmlncykgZGV2Lm9mZigpCmBgYAoK