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 
LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogRWxlY3Rpb25zIEVjb25vbXkgLS0gbW9kZWwgY2hlY2tpbmciCmF1dGhvcjogIkFuZHJldyBHZWxtYW4sIEplbm5pZmVyIEhpbGwsIEFraSBWZWh0YXJpIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiAyCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCkVsZWN0aW9ucyBFY29ub215IC0tIG1vZGVsIGNoZWNraW5nLiBDaGVja2luZyB0aGUgbW9kZWwtZml0dGluZwpwcm9jZWR1cmUgdXNpbmcgZmFrZS1kYXRhIHNpbXVsYXRpb24uIFNlZSBDaGFwdGVyIDcgaW4gUmVncmVzc2lvbgphbmQgT3RoZXIgU3Rvcmllcy4KCi0tLS0tLS0tLS0tLS0KCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBjb21tZW50PU5BKQpgYGAKCiMjIyMgTG9hZCBwYWNrYWdlcwoKYGBge3IgfQpsaWJyYXJ5KCJycHJvanJvb3QiKQpyb290PC1oYXNfZmlsZSgiLlJPUy1FeGFtcGxlcy1yb290IikkbWFrZV9maXhfZmlsZSgpCmxpYnJhcnkoInJzdGFuYXJtIikKYGBgCgojIyMjIExvYWQgZGF0YQoKYGBge3IgfQpoaWJicyA8LSByZWFkLnRhYmxlKHJvb3QoIkVsZWN0aW9uc0Vjb25vbXkvZGF0YSIsImhpYmJzLmRhdCIpLCBoZWFkZXI9VFJVRSkKaGVhZChoaWJicykKYGBgCgojIyMjIFN0ZXAgMTogQ3JlYXRpbmcgdGhlIHByZXRlbmQgd29ybGQKCmBgYHtyIH0KYSA8LSA0Ni4yCmIgPC0gMy4xCnNpZ21hIDwtIDMuOAp4IDwtIGhpYmJzJGdyb3d0aApuIDwtIGxlbmd0aCh4KQpgYGAKCiMjIyMgU3RlcCAyOiBTaW11bGF0aW5nIGZha2UgZGF0YQoKYGBge3IgfQpzZXQuc2VlZCgxKQp5IDwtIGEgKyBiKnggKyBybm9ybShuLCAwLCBzaWdtYSkKZmFrZSA8LSBkYXRhLmZyYW1lKHgsIHkpCmBgYAoKIyMjIyBTdGVwIDM6IEZpdHRpbmcgdGhlIG1vZGVsIGFuZCBjb21wYXJpbmcgZml0dGVkIHRvIHByZXRlbmQgdmFsdWVzCgpgYGB7ciByZXN1bHRzPSdoaWRlJ30KZml0IDwtIHN0YW5fZ2xtKHkgfiB4LCBkYXRhID0gZmFrZSwgcmVmcmVzaCA9IDApCmBgYApgYGB7ciB9CnByaW50KGZpdCkKCnBpNTAgPC0gcG9zdGVyaW9yX2ludGVydmFsKGZpdCwgcHJvYiA9IDAuNSlbJ3gnLF0KcGk5MCA8LSBwb3N0ZXJpb3JfaW50ZXJ2YWwoZml0LCBwcm9iID0gMC45KVsneCcsXQpjb3Zlcl81MCA8LSBhcy5udW1lcmljKGIgPiBwaTUwWzFdICYgYiA8IHBpNTBbMl0pCmNvdmVyXzkwIDwtIGFzLm51bWVyaWMoYiA+IHBpOTBbMV0gJiBiIDwgcGk5MFsyXSkKY2F0KHBhc3RlKCI1MCUgY292ZXJhZ2U6ICIsIGNvdmVyXzUwLCAiXG4iKSkKY2F0KHBhc3RlKCI5MCUgY292ZXJhZ2U6ICIsIGNvdmVyXzkwLCAiXG4iKSkKYGBgCgojIyMjIFN0ZXAgNDogIEVtYmVkZGluZyB0aGUgc2ltdWxhdGlvbiBpbiBhIGxvb3AKCmBgYHtyIHJlc3VsdHM9J2hpZGUnfQpuX2Zha2UgPC0gMTAwMApjb3Zlcl81MCA8LSByZXAoTkEsIG5fZmFrZSkKY292ZXJfNjggPC0gcmVwKE5BLCBuX2Zha2UpCmNvdmVyXzY4YiA8LSByZXAoTkEsIG5fZmFrZSkKY292ZXJfOTAgPC0gcmVwKE5BLCBuX2Zha2UpCmNvdmVyXzk1IDwtIHJlcChOQSwgbl9mYWtlKQpjb3Zlcl85NWIgPC0gcmVwKE5BLCBuX2Zha2UpCnBiIDwtIHR4dFByb2dyZXNzQmFyKG1pbj0wLCBtYXg9bl9mYWtlLCBpbml0aWFsPTAsIHN0eWxlPTMpCmZvciAocyBpbiAxOm5fZmFrZSl7CiAgc2V0VHh0UHJvZ3Jlc3NCYXIocGIsIHMpCiAgc2V0LnNlZWQocykKICB5IDwtIGEgKyBiKnggKyBybm9ybShuLCAwLCBzaWdtYSkKICBmYWtlIDwtIGRhdGEuZnJhbWUoeCwgeSkKICBmaXQgPC0gc3Rhbl9nbG0oeSB+IHgsIGRhdGEgPSBmYWtlLCB3YXJtdXAgPSA1MDAsIGl0ZXIgPSAxNTAwLCByZWZyZXNoID0gMCwKICAgICAgICAgICAgICAgICAgICAgIHNhdmVfd2FybXVwID0gRkFMU0UsIGNvcmVzID0gMSwgb3Blbl9wcm9ncmVzcyA9IEZBTFNFLAogICAgICAgICAgICAgICAgICAgICAgc2VlZCA9IHMpCiAgcGk1MCA8LSBwb3N0ZXJpb3JfaW50ZXJ2YWwoZml0LCBwcm9iID0gMC41KVsneCcsXQogIHBpNjggPC0gcG9zdGVyaW9yX2ludGVydmFsKGZpdCwgcHJvYiA9IDAuNjgpWyd4JyxdCiAgcGk5MCA8LSBwb3N0ZXJpb3JfaW50ZXJ2YWwoZml0LCBwcm9iID0gMC45KVsneCcsXQogIHBpOTUgPC0gcG9zdGVyaW9yX2ludGVydmFsKGZpdCwgcHJvYiA9IDAuOTUpWyd4JyxdCiAgY292ZXJfNTBbc10gPC0gYXMubnVtZXJpYyhiID4gcGk1MFsxXSAmIGIgPCBwaTUwWzJdKQogIGNvdmVyXzY4W3NdIDwtIGFzLm51bWVyaWMoYiA+IHBpNjhbMV0gJiBiIDwgcGk2OFsyXSkKICBjb3Zlcl85MFtzXSA8LSBhcy5udW1lcmljKGIgPiBwaTkwWzFdICYgYiA8IHBpOTBbMl0pCiAgY292ZXJfOTVbc10gPC0gYXMubnVtZXJpYyhiID4gcGk5NVsxXSAmIGIgPCBwaTk1WzJdKQogIGNvdmVyXzY4YltzXSA8LSBhcy5udW1lcmljKGFicyhiLWNvZWYoZml0KVsneCddKSA8IHNlKGZpdClbJ3gnXSkKICBjb3Zlcl85NWJbc10gPC0gYXMubnVtZXJpYyhhYnMoYi1jb2VmKGZpdClbJ3gnXSkgPCAyKnNlKGZpdClbJ3gnXSkKfQpjbG9zZShwYikKYGBgCmBgYHtyIH0KY2F0KHBhc3RlKCI1MCUgY292ZXJhZ2U6ICIsIG1lYW4oY292ZXJfNTApLCAiXG4iKSkKY2F0KHBhc3RlKCI2OCUgY292ZXJhZ2U6ICIsIG1lYW4oY292ZXJfNjgpLCAiXG4iKSkKY2F0KHBhc3RlKCI5MCUgY292ZXJhZ2U6ICIsIG1lYW4oY292ZXJfOTApLCAiXG4iKSkKY2F0KHBhc3RlKCI5NSUgY292ZXJhZ2U6ICIsIG1lYW4oY292ZXJfOTUpLCAiXG4iKSkKY2F0KHBhc3RlKCI2OCUgY292ZXJhZ2U6ICIsIG1lYW4oY292ZXJfNjhiKSwgIlxuIikpCmNhdChwYXN0ZSgiOTUlIGNvdmVyYWdlOiAiLCBtZWFuKGNvdmVyXzk1YiksICJcbiIpKQpgYGAKCg==