Example where an informative prior makes a difference. See Chapter 9 in Regression and Other Stories.


Load packages

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

Data

x <- seq(-2,2,1)
y <- c(50, 44, 50, 47, 56)
sexratio <- data.frame(x, y)

Informative priors

theta_hat_prior <- 0
se_prior <- 0.25
theta_hat_data <- 8
se_data <- 3
theta_hat_bayes <- (theta_hat_prior/se_prior^2 + theta_hat_data/se_data^2)/(1/se_prior^2 + 1/se_data^2)
se_bayes <- sqrt(1/(1/se_prior^2 + 1/se_data^2))

Least-squares regression

fit <- lm(y ~ x, data = sexratio)
display(fit)
lm(formula = y ~ x, data = sexratio)
            coef.est coef.se
(Intercept) 49.40     1.94  
x            1.50     1.37  
---
n = 5, k = 2
residual sd = 4.35, R-Squared = 0.28

Plot data and least-squares regression line

par(mfrow=c(1,2), mar=c(3,3,3,2), mgp=c(1.7,.5,0), tck=-.01)
plot(x, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="Data on beauty and sex ratio",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
plot(x, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="Data and least-squares regression line",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
abline(coef(fit)[1], coef(fit)[2])
text(1, 52.2, paste("y = ", fround(coef(fit)[1], 1), " + ", fround(coef(fit)[2], 1), " x\n(Std err of slope is ", fround(se.coef(fit)[2], 1), ")", sep=""))

Bayesian regression with weakly informative prior

fit_default <- stan_glm(y ~ x, data = sexratio, refresh = 0)
prior_summary(fit_default)
Priors for model 'fit_default' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 49, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 49, scale = 11)

Coefficients
  Specified prior:
    ~ normal(location = 0, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 0, scale = 7)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.22)
------
See help('prior_summary.stanreg') for more details
print(fit_default)
stan_glm
 family:       gaussian [identity]
 formula:      y ~ x
 observations: 5
 predictors:   2
------
            Median MAD_SD
(Intercept) 49.4    2.0  
x            1.4    1.4  

Auxiliary parameter(s):
      Median MAD_SD
sigma 4.6    1.7   

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

Bayesian regression with informative prior

fit_post <- stan_glm(y ~ x, data = sexratio,
                     prior = normal(0, 0.2),
                     prior_intercept = normal(48.8, 0.5),
                     refresh = 0)
prior_summary(fit_post)
Priors for model 'fit_post' 
------
Intercept (after predictors centered)
 ~ normal(location = 49, scale = 0.5)

Coefficients
 ~ normal(location = 0, scale = 0.2)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.22)
------
See help('prior_summary.stanreg') for more details
print(fit_post)
stan_glm
 family:       gaussian [identity]
 formula:      y ~ x
 observations: 5
 predictors:   2
------
            Median MAD_SD
(Intercept) 48.8    0.5  
x            0.0    0.2  

Auxiliary parameter(s):
      Median MAD_SD
sigma 4.3    1.3   

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

Plot Posterior simulations under weakly informative and informative prior

par(mfrow=c(2,2), mar=c(5,3,3,2), mgp=c(1.7,.5,0), tck=-.01)
fit_bayes <- list(as.data.frame(fit_default), as.data.frame(fit_post))
for (k in 1:2){
  sims <- fit_bayes[[k]]
  coef_est <- apply(sims, 2, median)
  b_se <- 1.483*median(abs(sims[,2] - median(sims[,2])))
  plot(sims[,1], sims[,2], xlim=range(fit_bayes[[1]][,1], fit_bayes[[2]][,1]), ylim=range(fit_bayes[[1]][,2], fit_bayes[[2]][,2]), xlab="Intercept, a", ylab="Slope, b", main=paste("Posterior simulations under", if (k==1) "default prior" else "informative prior"), pch=19, cex=.2, bty="l")
  plot(x, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main=if (k==1) "Least-squares regression line and\nposterior uncertainty given default prior" else "Bayes estimated regression line and\nposterior uncertainty given informative prior", pch=19, cex=1)
  axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
  for (i in 1:100){
    abline(sims[i,1], sims[i,2], lwd=.5, col="gray50")
  }
  abline(coef_est[1], coef_est[2], lwd=2)
}

LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogQmVhdXR5IGFuZCBzZXggcmF0aW8iCmF1dGhvcjogIkFuZHJldyBHZWxtYW4sIEplbm5pZmVyIEhpbGwsIEFraSBWZWh0YXJpIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiAyCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCkV4YW1wbGUgd2hlcmUgYW4gaW5mb3JtYXRpdmUgcHJpb3IgbWFrZXMgYSBkaWZmZXJlbmNlLiBTZWUgQ2hhcHRlciA5CmluIFJlZ3Jlc3Npb24gYW5kIE90aGVyIFN0b3JpZXMuCgotLS0tLS0tLS0tLS0tCgoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRSwgY29tbWVudD1OQSkKIyBzd2l0Y2ggdGhpcyB0byBUUlVFIHRvIHNhdmUgZmlndXJlcyBpbiBzZXBhcmF0ZSBmaWxlcwpzYXZlZmlncyA8LSBGQUxTRQpgYGAKCiMjIyMgTG9hZCBwYWNrYWdlcwoKYGBge3IgfQpsaWJyYXJ5KCJycHJvanJvb3QiKQpyb290PC1oYXNfZmlsZSgiLlJPUy1FeGFtcGxlcy1yb290IikkbWFrZV9maXhfZmlsZSgpCmxpYnJhcnkoImFybSIpCmxpYnJhcnkoInJzdGFuYXJtIikKYGBgCgojIyBEYXRhCgpgYGB7ciB9CnggPC0gc2VxKC0yLDIsMSkKeSA8LSBjKDUwLCA0NCwgNTAsIDQ3LCA1NikKc2V4cmF0aW8gPC0gZGF0YS5mcmFtZSh4LCB5KQpgYGAKCiMjIEluZm9ybWF0aXZlIHByaW9ycwoKYGBge3IgfQp0aGV0YV9oYXRfcHJpb3IgPC0gMApzZV9wcmlvciA8LSAwLjI1CnRoZXRhX2hhdF9kYXRhIDwtIDgKc2VfZGF0YSA8LSAzCnRoZXRhX2hhdF9iYXllcyA8LSAodGhldGFfaGF0X3ByaW9yL3NlX3ByaW9yXjIgKyB0aGV0YV9oYXRfZGF0YS9zZV9kYXRhXjIpLygxL3NlX3ByaW9yXjIgKyAxL3NlX2RhdGFeMikKc2VfYmF5ZXMgPC0gc3FydCgxLygxL3NlX3ByaW9yXjIgKyAxL3NlX2RhdGFeMikpCmBgYAoKIyMgTGVhc3Qtc3F1YXJlcyByZWdyZXNzaW9uCgpgYGB7ciB9CmZpdCA8LSBsbSh5IH4geCwgZGF0YSA9IHNleHJhdGlvKQpkaXNwbGF5KGZpdCkKYGBgCgojIyMjIFBsb3QgZGF0YSBhbmQgbGVhc3Qtc3F1YXJlcyByZWdyZXNzaW9uIGxpbmUKCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgcGRmKHJvb3QoIlNleFJhdGlvL2ZpZ3MiLCJzZXhyYXRpb19iYXllc18xLnBkZiIpLCBoZWlnaHQ9NCwgd2lkdGg9MTApCmBgYApgYGB7ciB9CnBhcihtZnJvdz1jKDEsMiksIG1hcj1jKDMsMywzLDIpLCBtZ3A9YygxLjcsLjUsMCksIHRjaz0tLjAxKQpwbG90KHgsIHksIHlsaW09Yyg0MywgNTcpLCB4bGFiPSJBdHRyYWN0aXZlbmVzcyBvZiBwYXJlbnQiLCB5bGFiPSJQZXJjZW50YWdlIG9mIGdpcmwgYmFiaWVzIiwgYnR5PSJsIiwgeWF4dD0ibiIsIG1haW49IkRhdGEgb24gYmVhdXR5IGFuZCBzZXggcmF0aW8iLCAgcGNoPTE5LCBjZXg9MSkKYXhpcygyLCBjKDQ1LDUwLDU1KSwgcGFzdGUoYyg0NSw1MCw1NSksICIlIiwgc2VwPSIiKSkKcGxvdCh4LCB5LCB5bGltPWMoNDMsIDU3KSwgeGxhYj0iQXR0cmFjdGl2ZW5lc3Mgb2YgcGFyZW50IiwgeWxhYj0iUGVyY2VudGFnZSBvZiBnaXJsIGJhYmllcyIsIGJ0eT0ibCIsIHlheHQ9Im4iLCBtYWluPSJEYXRhIGFuZCBsZWFzdC1zcXVhcmVzIHJlZ3Jlc3Npb24gbGluZSIsICBwY2g9MTksIGNleD0xKQpheGlzKDIsIGMoNDUsNTAsNTUpLCBwYXN0ZShjKDQ1LDUwLDU1KSwgIiUiLCBzZXA9IiIpKQphYmxpbmUoY29lZihmaXQpWzFdLCBjb2VmKGZpdClbMl0pCnRleHQoMSwgNTIuMiwgcGFzdGUoInkgPSAiLCBmcm91bmQoY29lZihmaXQpWzFdLCAxKSwgIiArICIsIGZyb3VuZChjb2VmKGZpdClbMl0sIDEpLCAiIHhcbihTdGQgZXJyIG9mIHNsb3BlIGlzICIsIGZyb3VuZChzZS5jb2VmKGZpdClbMl0sIDEpLCAiKSIsIHNlcD0iIikpCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIGRldi5vZmYoKQpgYGAKCiMjIEJheWVzaWFuIHJlZ3Jlc3Npb24gd2l0aCB3ZWFrbHkgaW5mb3JtYXRpdmUgcHJpb3IKCmBgYHtyIH0KZml0X2RlZmF1bHQgPC0gc3Rhbl9nbG0oeSB+IHgsIGRhdGEgPSBzZXhyYXRpbywgcmVmcmVzaCA9IDApCnByaW9yX3N1bW1hcnkoZml0X2RlZmF1bHQpCnByaW50KGZpdF9kZWZhdWx0KQpgYGAKCiMjIEJheWVzaWFuIHJlZ3Jlc3Npb24gd2l0aCBpbmZvcm1hdGl2ZSBwcmlvcgoKYGBge3IgfQpmaXRfcG9zdCA8LSBzdGFuX2dsbSh5IH4geCwgZGF0YSA9IHNleHJhdGlvLAogICAgICAgICAgICAgICAgICAgICBwcmlvciA9IG5vcm1hbCgwLCAwLjIpLAogICAgICAgICAgICAgICAgICAgICBwcmlvcl9pbnRlcmNlcHQgPSBub3JtYWwoNDguOCwgMC41KSwKICAgICAgICAgICAgICAgICAgICAgcmVmcmVzaCA9IDApCnByaW9yX3N1bW1hcnkoZml0X3Bvc3QpCnByaW50KGZpdF9wb3N0KQpgYGAKCiMjIyMgUGxvdCBQb3N0ZXJpb3Igc2ltdWxhdGlvbnMgdW5kZXIgd2Vha2x5IGluZm9ybWF0aXZlIGFuZCBpbmZvcm1hdGl2ZSBwcmlvcgoKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBwZGYocm9vdCgiU2V4UmF0aW8vZmlncyIsInNleHJhdGlvX2JheWVzXzIucGRmIiksIGhlaWdodD04LCB3aWR0aD0xMCkKYGBgCmBgYHtyIH0KcGFyKG1mcm93PWMoMiwyKSwgbWFyPWMoNSwzLDMsMiksIG1ncD1jKDEuNywuNSwwKSwgdGNrPS0uMDEpCmZpdF9iYXllcyA8LSBsaXN0KGFzLmRhdGEuZnJhbWUoZml0X2RlZmF1bHQpLCBhcy5kYXRhLmZyYW1lKGZpdF9wb3N0KSkKZm9yIChrIGluIDE6Mil7CiAgc2ltcyA8LSBmaXRfYmF5ZXNbW2tdXQogIGNvZWZfZXN0IDwtIGFwcGx5KHNpbXMsIDIsIG1lZGlhbikKICBiX3NlIDwtIDEuNDgzKm1lZGlhbihhYnMoc2ltc1ssMl0gLSBtZWRpYW4oc2ltc1ssMl0pKSkKICBwbG90KHNpbXNbLDFdLCBzaW1zWywyXSwgeGxpbT1yYW5nZShmaXRfYmF5ZXNbWzFdXVssMV0sIGZpdF9iYXllc1tbMl1dWywxXSksIHlsaW09cmFuZ2UoZml0X2JheWVzW1sxXV1bLDJdLCBmaXRfYmF5ZXNbWzJdXVssMl0pLCB4bGFiPSJJbnRlcmNlcHQsIGEiLCB5bGFiPSJTbG9wZSwgYiIsIG1haW49cGFzdGUoIlBvc3RlcmlvciBzaW11bGF0aW9ucyB1bmRlciIsIGlmIChrPT0xKSAiZGVmYXVsdCBwcmlvciIgZWxzZSAiaW5mb3JtYXRpdmUgcHJpb3IiKSwgcGNoPTE5LCBjZXg9LjIsIGJ0eT0ibCIpCiAgcGxvdCh4LCB5LCB5bGltPWMoNDMsIDU3KSwgeGxhYj0iQXR0cmFjdGl2ZW5lc3Mgb2YgcGFyZW50IiwgeWxhYj0iUGVyY2VudGFnZSBvZiBnaXJsIGJhYmllcyIsIGJ0eT0ibCIsIHlheHQ9Im4iLCBtYWluPWlmIChrPT0xKSAiTGVhc3Qtc3F1YXJlcyByZWdyZXNzaW9uIGxpbmUgYW5kXG5wb3N0ZXJpb3IgdW5jZXJ0YWludHkgZ2l2ZW4gZGVmYXVsdCBwcmlvciIgZWxzZSAiQmF5ZXMgZXN0aW1hdGVkIHJlZ3Jlc3Npb24gbGluZSBhbmRcbnBvc3RlcmlvciB1bmNlcnRhaW50eSBnaXZlbiBpbmZvcm1hdGl2ZSBwcmlvciIsIHBjaD0xOSwgY2V4PTEpCiAgYXhpcygyLCBjKDQ1LDUwLDU1KSwgcGFzdGUoYyg0NSw1MCw1NSksICIlIiwgc2VwPSIiKSkKICBmb3IgKGkgaW4gMToxMDApewogICAgYWJsaW5lKHNpbXNbaSwxXSwgc2ltc1tpLDJdLCBsd2Q9LjUsIGNvbD0iZ3JheTUwIikKICB9CiAgYWJsaW5lKGNvZWZfZXN0WzFdLCBjb2VmX2VzdFsyXSwgbHdkPTIpCn0KYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgZGV2Lm9mZigpCmBgYAoK