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==