Regression and Other Stories: Elections Economy -- model checking

Elections Economy -- model checking. Checking the model-fitting procedure using fake-data simulation. See Chapter 7 in Regression and Other Stories.

---

Load packages

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

Load data

hibbs <- read.table(root("ElectionsEconomy/data","hibbs.dat"), header=TRUE)
head(hibbs)

  year growth  vote inc_party_candidate other_candidate
1 1952   2.40 44.60           Stevenson      Eisenhower
2 1956   2.89 57.76          Eisenhower       Stevenson
3 1960   0.85 49.91               Nixon         Kennedy
4 1964   4.21 61.34             Johnson       Goldwater
5 1968   3.02 49.60            Humphrey           Nixon
6 1972   3.62 61.79               Nixon        McGovern

Step 1: Creating the pretend world

a <- 46.2
b <- 3.1
sigma <- 3.8
x <- hibbs$growth
n <- length(x)

Step 2: Simulating fake data

set.seed(1)
y <- a + b*x + rnorm(n, 0, sigma)
fake <- data.frame(x, y)

Step 3: Fitting the model and comparing fitted to pretend values

fit <- stan_glm(y ~ x, data = fake, refresh = 0)

print(fit)

stan_glm
 family:       gaussian [identity]
 formula:      y ~ x
 observations: 16
 predictors:   2
------
            Median MAD_SD
(Intercept) 47.1    1.7  
x            2.8    0.8  

Auxiliary parameter(s):
      Median MAD_SD
sigma 4.0    0.7   

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

pi50 <- posterior_interval(fit, prob = 0.5)['x',]
pi90 <- posterior_interval(fit, prob = 0.9)['x',]
cover_50 <- as.numeric(b > pi50[1] & b < pi50[2])
cover_90 <- as.numeric(b > pi90[1] & b < pi90[2])
cat(paste("50% coverage: ", cover_50, "\n"))

50% coverage:  1 

cat(paste("90% coverage: ", cover_90, "\n"))

90% coverage:  1 

Step 4: Embedding the simulation in a loop

n_fake <- 1000
cover_50 <- rep(NA, n_fake)
cover_68 <- rep(NA, n_fake)
cover_68b <- rep(NA, n_fake)
cover_90 <- rep(NA, n_fake)
cover_95 <- rep(NA, n_fake)
cover_95b <- rep(NA, n_fake)
pb <- txtProgressBar(min=0, max=n_fake, initial=0, style=3)
for (s in 1:n_fake){
  setTxtProgressBar(pb, s)
  set.seed(s)
  y <- a + b*x + rnorm(n, 0, sigma)
  fake <- data.frame(x, y)
  fit <- stan_glm(y ~ x, data = fake, warmup = 500, iter = 1500, refresh = 0,
                      save_warmup = FALSE, cores = 1, open_progress = FALSE,
                      seed = s)
  pi50 <- posterior_interval(fit, prob = 0.5)['x',]
  pi68 <- posterior_interval(fit, prob = 0.68)['x',]
  pi90 <- posterior_interval(fit, prob = 0.9)['x',]
  pi95 <- posterior_interval(fit, prob = 0.95)['x',]
  cover_50[s] <- as.numeric(b > pi50[1] & b < pi50[2])
  cover_68[s] <- as.numeric(b > pi68[1] & b < pi68[2])
  cover_90[s] <- as.numeric(b > pi90[1] & b < pi90[2])
  cover_95[s] <- as.numeric(b > pi95[1] & b < pi95[2])
  cover_68b[s] <- as.numeric(abs(b-coef(fit)['x']) < se(fit)['x'])
  cover_95b[s] <- as.numeric(abs(b-coef(fit)['x']) < 2*se(fit)['x'])
}
close(pb)

cat(paste("50% coverage: ", mean(cover_50), "\n"))

50% coverage:  0.516 

cat(paste("68% coverage: ", mean(cover_68), "\n"))

68% coverage:  0.714 

cat(paste("90% coverage: ", mean(cover_90), "\n"))

90% coverage:  0.911 

cat(paste("95% coverage: ", mean(cover_95), "\n"))

95% coverage:  0.949 

cat(paste("68% coverage: ", mean(cover_68b), "\n"))

68% coverage:  0.71 

cat(paste("95% coverage: ", mean(cover_95b), "\n"))

95% coverage:  0.943