Predictive uncertainty for congressional elections. See Chapters 10 and 15 in Regression and Other Stories.


Load packages

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

Load data

congress <- read.csv(root("Congress/data","congress.csv"))
inconsistent <- (congress$inc88==(-1) & congress$v86 > 0.5) |  (congress$inc88==1 & congress$v86 < 0.5)
head(congress)
  inc86 inc88 inc90       v86       v88       v90   v86_adj   v88_adj   v90_adj
1     1     1     1 0.7450362 0.7724427 0.7140294 0.7450362 0.7724427 0.7140294
2     1     1     1 0.6738455 0.6361816 0.5970501 0.6738455 0.6361816 0.5970501
3     1     1     0 0.6964566 0.6649283 0.5210433 0.6964566 0.6649283 0.5210433
4    -1    -1    -1 0.4645901 0.2738342 0.2343770 0.4645901 0.2738342 0.2343770
5    -1    -1     0 0.3910945 0.2636131 0.4774393 0.3910945 0.2636131 0.4774393
6    -1    -1    -1 0.3582454 0.3341927 0.2562970 0.3582454 0.3341927 0.2562970

Regression predicting 1988 from 1986

data88 <- data.frame(vote=congress$v88_adj, past_vote=congress$v86_adj, inc=congress$inc88)
fit88 <- stan_glm(vote ~ past_vote + inc, data=data88, refresh=0)
print(fit88, digits=2)
stan_glm
 family:       gaussian [identity]
 formula:      vote ~ past_vote + inc
 observations: 435
 predictors:   3
------
            Median MAD_SD
(Intercept) 0.24   0.02  
past_vote   0.52   0.03  
inc         0.10   0.01  

Auxiliary parameter(s):
      Median MAD_SD
sigma 0.07   0.00  

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

Simulation for inferences and predictions of new data points

Predict from 1988 to 1990

data90 <- data.frame(past_vote=congress$v88_adj, inc=congress$inc90)

Simulate predictive simulations of the vector of new outcomes

pred90 <- posterior_predict(fit88, newdata=data90)

Simulate the number of elections predicted to be won by the Democrats in 1990

dems_pred <- rowSums(pred90 > 0.5)

Alternately calculate that sum in a loop

n_sims <- 4000
dems_pred <- rep(NA, n_sims)
for (s in 1:n_sims) {
  dems_pred[s] <- sum(pred90[s,] > 0.5)
}

Our posterior mean and sd of how many districts the Dems will win

print(c(mean(dems_pred), sd(dems_pred)), digits=2)
[1] 260.0   2.5

Histogram of how many districts the Dems will win

hist(dems_pred)

Graphs

par(mar=c(3,0,1,0), mgp=c(1.5, .5, 0), tck=-.01)
v88_hist <- ifelse(congress$v88<.1, .0001, ifelse(congress$v88>.9, .9999, congress$v88))
hist(v88_hist, breaks=seq(0,1,.05),
     xlab="Democratic share of the two-party vote", ylab="", yaxt="n", main="")
mtext("Congressional elections in 1988", 3, 0)

if (savefigs) dev.off()

if (savefigs) pdf("cong.pdf", height=4, width=4, colormodel="gray")
par(pty="s")
par(mar=c(3,0,3,0), mgp=c(1.7, .5, 0), tck=-.01)
plot(0, 0, xlim=c(0,1), ylim=c(0,1), type="n", xaxs="i", yaxs="i",
  xlab="Democratic vote share in 1986", ylab="Democratic vote share in 1988")
abline(0,1, lwd=.5)
jitt <- function(vote){
  n <- length(vote)
  ifelse(vote<0.1, runif(n, 0.01, 0.04), ifelse(vote>0.9, runif(n, 0.96, 0.99), vote))
}
j_v86 <- jitt(congress$v86)
j_v88 <- jitt(congress$v88)
points(j_v86[congress$inc88==0], j_v88[congress$inc88==0], pch=1, cex=1)
points(j_v86[congress$inc88==1], j_v88[congress$inc88==1], pch=16, cex=.8)
points(j_v86[congress$inc88==-1], j_v88[congress$inc88==-1], pch=4, cex=.8)
mtext("Raw data", 3, .7)

par(pty="s")
par(mar=c(3,0,3,0), mgp=c(1.7, .5, 0), tck=-.01)
plot(0, 0, xlim=c(0,1), ylim=c(0,1), type="n", xaxs="i", yaxs="i",
  xlab="Adjusted Dem. vote share in 1986", ylab="Adjusted Dem. vote share in 1988")
abline(0,1, lwd=.5)
points(congress$v86_adj[congress$inc88==0], congress$v88_adj[congress$inc88==0], pch=1, cex=1)
points(congress$v86_adj[congress$inc88==1], congress$v88_adj[congress$inc88==1], pch=16, cex=.8)
points(congress$v86_adj[congress$inc88==-1], congress$v88_adj[congress$inc88==-1], pch=4, cex=.8)
mtext("Adjusted data", 3, .7)

LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogQ29uZ3Jlc3MiCmF1dGhvcjogIkFuZHJldyBHZWxtYW4sIEplbm5pZmVyIEhpbGwsIEFraSBWZWh0YXJpIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiAyCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tClByZWRpY3RpdmUgdW5jZXJ0YWludHkgZm9yIGNvbmdyZXNzaW9uYWwgZWxlY3Rpb25zLiBTZWUgQ2hhcHRlcnMgMTAKYW5kIDE1IGluIFJlZ3Jlc3Npb24gYW5kIE90aGVyIFN0b3JpZXMuCgotLS0tLS0tLS0tLS0tCgoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRSwgY29tbWVudD1OQSkKIyBzd2l0Y2ggdGhpcyB0byBUUlVFIHRvIHNhdmUgZmlndXJlcyBpbiBzZXBhcmF0ZSBmaWxlcwpzYXZlZmlncyA8LSBGQUxTRQpgYGAKCiMjIyMgTG9hZCBwYWNrYWdlcwoKYGBge3IgfQpsaWJyYXJ5KCJycHJvanJvb3QiKQpyb290PC1oYXNfZmlsZSgiLlJPUy1FeGFtcGxlcy1yb290IikkbWFrZV9maXhfZmlsZSgpCmxpYnJhcnkoInJzdGFuYXJtIikKYGBgCgojIyMjIExvYWQgZGF0YQoKYGBge3IgfQpjb25ncmVzcyA8LSByZWFkLmNzdihyb290KCJDb25ncmVzcy9kYXRhIiwiY29uZ3Jlc3MuY3N2IikpCmluY29uc2lzdGVudCA8LSAoY29uZ3Jlc3MkaW5jODg9PSgtMSkgJiBjb25ncmVzcyR2ODYgPiAwLjUpIHwgIChjb25ncmVzcyRpbmM4OD09MSAmIGNvbmdyZXNzJHY4NiA8IDAuNSkKaGVhZChjb25ncmVzcykKYGBgCgojIyBSZWdyZXNzaW9uIHByZWRpY3RpbmcgMTk4OCBmcm9tIDE5ODYKCmBgYHtyIH0KZGF0YTg4IDwtIGRhdGEuZnJhbWUodm90ZT1jb25ncmVzcyR2ODhfYWRqLCBwYXN0X3ZvdGU9Y29uZ3Jlc3Mkdjg2X2FkaiwgaW5jPWNvbmdyZXNzJGluYzg4KQpmaXQ4OCA8LSBzdGFuX2dsbSh2b3RlIH4gcGFzdF92b3RlICsgaW5jLCBkYXRhPWRhdGE4OCwgcmVmcmVzaD0wKQpwcmludChmaXQ4OCwgZGlnaXRzPTIpCmBgYAoKIyMgU2ltdWxhdGlvbiBmb3IgaW5mZXJlbmNlcyBhbmQgcHJlZGljdGlvbnMgb2YgbmV3IGRhdGEgcG9pbnRzCiMjIyMgUHJlZGljdCBmcm9tIDE5ODggdG8gMTk5MAoKYGBge3IgfQpkYXRhOTAgPC0gZGF0YS5mcmFtZShwYXN0X3ZvdGU9Y29uZ3Jlc3Mkdjg4X2FkaiwgaW5jPWNvbmdyZXNzJGluYzkwKQpgYGAKCiMjIyMgU2ltdWxhdGUgcHJlZGljdGl2ZSBzaW11bGF0aW9ucyBvZiB0aGUgdmVjdG9yIG9mIG5ldyBvdXRjb21lcwoKYGBge3IgfQpwcmVkOTAgPC0gcG9zdGVyaW9yX3ByZWRpY3QoZml0ODgsIG5ld2RhdGE9ZGF0YTkwKQpgYGAKCiMjIyMgU2ltdWxhdGUgdGhlIG51bWJlciBvZiBlbGVjdGlvbnMgcHJlZGljdGVkIHRvIGJlIHdvbiBieSB0aGUgRGVtb2NyYXRzIGluIDE5OTAKCmBgYHtyIH0KZGVtc19wcmVkIDwtIHJvd1N1bXMocHJlZDkwID4gMC41KQpgYGAKCkFsdGVybmF0ZWx5IGNhbGN1bGF0ZSB0aGF0IHN1bSBpbiBhIGxvb3AKCmBgYHtyIH0Kbl9zaW1zIDwtIDQwMDAKZGVtc19wcmVkIDwtIHJlcChOQSwgbl9zaW1zKQpmb3IgKHMgaW4gMTpuX3NpbXMpIHsKICBkZW1zX3ByZWRbc10gPC0gc3VtKHByZWQ5MFtzLF0gPiAwLjUpCn0KYGBgCgojIyMjIE91ciBwb3N0ZXJpb3IgbWVhbiBhbmQgc2Qgb2YgaG93IG1hbnkgZGlzdHJpY3RzIHRoZSBEZW1zIHdpbGwgd2luCgpgYGB7ciB9CnByaW50KGMobWVhbihkZW1zX3ByZWQpLCBzZChkZW1zX3ByZWQpKSwgZGlnaXRzPTIpCmBgYAoKIyMjIyBIaXN0b2dyYW0gb2YgaG93IG1hbnkgZGlzdHJpY3RzIHRoZSBEZW1zIHdpbGwgd2luCgpgYGB7ciB9Cmhpc3QoZGVtc19wcmVkKQpgYGAKCiMjIEdyYXBocwoKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBwZGYoImhpc3Q4OC5wZGYiLCBoZWlnaHQ9My4yLCB3aWR0aD00LjEsIGNvbG9ybW9kZWw9ImdyYXkiKQpgYGAKYGBge3IgfQpwYXIobWFyPWMoMywwLDEsMCksIG1ncD1jKDEuNSwgLjUsIDApLCB0Y2s9LS4wMSkKdjg4X2hpc3QgPC0gaWZlbHNlKGNvbmdyZXNzJHY4ODwuMSwgLjAwMDEsIGlmZWxzZShjb25ncmVzcyR2ODg+LjksIC45OTk5LCBjb25ncmVzcyR2ODgpKQpoaXN0KHY4OF9oaXN0LCBicmVha3M9c2VxKDAsMSwuMDUpLAogICAgIHhsYWI9IkRlbW9jcmF0aWMgc2hhcmUgb2YgdGhlIHR3by1wYXJ0eSB2b3RlIiwgeWxhYj0iIiwgeWF4dD0ibiIsIG1haW49IiIpCm10ZXh0KCJDb25ncmVzc2lvbmFsIGVsZWN0aW9ucyBpbiAxOTg4IiwgMywgMCkKaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKCmlmIChzYXZlZmlncykgcGRmKCJjb25nLnBkZiIsIGhlaWdodD00LCB3aWR0aD00LCBjb2xvcm1vZGVsPSJncmF5IikKcGFyKHB0eT0icyIpCnBhcihtYXI9YygzLDAsMywwKSwgbWdwPWMoMS43LCAuNSwgMCksIHRjaz0tLjAxKQpwbG90KDAsIDAsIHhsaW09YygwLDEpLCB5bGltPWMoMCwxKSwgdHlwZT0ibiIsIHhheHM9ImkiLCB5YXhzPSJpIiwKICB4bGFiPSJEZW1vY3JhdGljIHZvdGUgc2hhcmUgaW4gMTk4NiIsIHlsYWI9IkRlbW9jcmF0aWMgdm90ZSBzaGFyZSBpbiAxOTg4IikKYWJsaW5lKDAsMSwgbHdkPS41KQpqaXR0IDwtIGZ1bmN0aW9uKHZvdGUpewogIG4gPC0gbGVuZ3RoKHZvdGUpCiAgaWZlbHNlKHZvdGU8MC4xLCBydW5pZihuLCAwLjAxLCAwLjA0KSwgaWZlbHNlKHZvdGU+MC45LCBydW5pZihuLCAwLjk2LCAwLjk5KSwgdm90ZSkpCn0Kal92ODYgPC0gaml0dChjb25ncmVzcyR2ODYpCmpfdjg4IDwtIGppdHQoY29uZ3Jlc3Mkdjg4KQpwb2ludHMoal92ODZbY29uZ3Jlc3MkaW5jODg9PTBdLCBqX3Y4OFtjb25ncmVzcyRpbmM4OD09MF0sIHBjaD0xLCBjZXg9MSkKcG9pbnRzKGpfdjg2W2NvbmdyZXNzJGluYzg4PT0xXSwgal92ODhbY29uZ3Jlc3MkaW5jODg9PTFdLCBwY2g9MTYsIGNleD0uOCkKcG9pbnRzKGpfdjg2W2NvbmdyZXNzJGluYzg4PT0tMV0sIGpfdjg4W2NvbmdyZXNzJGluYzg4PT0tMV0sIHBjaD00LCBjZXg9LjgpCm10ZXh0KCJSYXcgZGF0YSIsIDMsIC43KQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZigiY29uZ2NsZWFuLnBkZiIsIGhlaWdodD00LCB3aWR0aD00LCBjb2xvcm1vZGVsPSJncmF5IikKYGBgCmBgYHtyIH0KcGFyKHB0eT0icyIpCnBhcihtYXI9YygzLDAsMywwKSwgbWdwPWMoMS43LCAuNSwgMCksIHRjaz0tLjAxKQpwbG90KDAsIDAsIHhsaW09YygwLDEpLCB5bGltPWMoMCwxKSwgdHlwZT0ibiIsIHhheHM9ImkiLCB5YXhzPSJpIiwKICB4bGFiPSJBZGp1c3RlZCBEZW0uIHZvdGUgc2hhcmUgaW4gMTk4NiIsIHlsYWI9IkFkanVzdGVkIERlbS4gdm90ZSBzaGFyZSBpbiAxOTg4IikKYWJsaW5lKDAsMSwgbHdkPS41KQpwb2ludHMoY29uZ3Jlc3Mkdjg2X2Fkaltjb25ncmVzcyRpbmM4OD09MF0sIGNvbmdyZXNzJHY4OF9hZGpbY29uZ3Jlc3MkaW5jODg9PTBdLCBwY2g9MSwgY2V4PTEpCnBvaW50cyhjb25ncmVzcyR2ODZfYWRqW2NvbmdyZXNzJGluYzg4PT0xXSwgY29uZ3Jlc3Mkdjg4X2Fkaltjb25ncmVzcyRpbmM4OD09MV0sIHBjaD0xNiwgY2V4PS44KQpwb2ludHMoY29uZ3Jlc3Mkdjg2X2Fkaltjb25ncmVzcyRpbmM4OD09LTFdLCBjb25ncmVzcyR2ODhfYWRqW2NvbmdyZXNzJGluYzg4PT0tMV0sIHBjaD00LCBjZXg9LjgpCm10ZXh0KCJBZGp1c3RlZCBkYXRhIiwgMywgLjcpCmBgYApgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIGRldi5vZmYoKQpgYGAKCg==