Regression and Other Stories: Mile

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")