Using data on voting behavior in European elections from http://www.jkarreth.net/files/full_cses3.dta, estimate a Bayesian logistic regression that predicts whether someone voted or not in 2007 as a function of their age, gender, education, and income. Present point estimates and measures of uncertainty for coefficient posterior distributions in a table. Calculate the predicted probability of voting as a function of age, and present in graphical form.
If you donâ€™t want to write your own Stan model, thereâ€™s one available on my GitHub. Johannes Karreth has helpfully written a function that automatically performs predicted probability calcuations on the output of a Bayesian model. Combining the summary()
function with kable()
can produce a serviceable regression table, but Iâ€™ve written a function which makes much nicer looking tables with minimal effort.
library(rstan) # Bayesian inference
rstan_options(auto_write = TRUE) # avoid recompiling model
options(mc.cores = parallel::detectCores()) # parallelizing multiple chains
library(xtable) # table output
# read in data from Johannes Karreth's website and subset to 2007
elections <- rio::import('http://www.jkarreth.net/files/full_cses3.dta')
elections <- elections[grepl('2007', elections$election), ]
# drop NAs in right and left side variables
elections <- na.omit(elections[, c("voted", "age", "female", "education",
"income")])
# response variable
Y <- elections$voted
# predictors
X <- model.matrix(voted ~ scale(age) + female + as.factor(education)
+ as.factor(income), data = elections)
# create named list of data to pass to Stan
stan_data <- list(Y = Y, X = X, N = nrow(X), K = ncol(X))
# download Stan model
download.file('https://raw.githubusercontent.com/jayrobwilliams/JRWmisc/master/logit.stan',
destfile = 'logit.stan')
# fit model
vote_fit <- stan(file = 'logit.stan', data = stan_data, pars = 'beta',
iter = 1000, chains = 4)
# rename betas
names(vote_fit)[1:14] <- gsub('as.factor', '', colnames(X))
# source MCMC regression table function
source('https://raw.githubusercontent.com/jayrobwilliams/RWmisc/master/mcmcreg.R')
# create regression table
mcmcreg(vote_fit, pars = 'beta', format = 'html', caption = 'Bayesian logistic regression.
Point estimates are posterior means.')
Model 1 | ||
---|---|---|
(Intercept) | 0.47 | |
[-0.02; 0.98] | ||
scale(age) | 0.40^{*} | |
[0.35; 0.45] | ||
female | -0.28^{*} | |
[-0.37; -0.19] | ||
(education)2 | 0.57 | |
[-0.03; 1.16] | ||
(education)3 | 0.97^{*} | |
[0.45; 1.47] | ||
(education)4 | 0.52 | |
[-0.01; 1.02] | ||
(education)5 | 0.70^{*} | |
[0.16; 1.19] | ||
(education)6 | 0.53^{*} | |
[0.01; 1.02] | ||
(education)7 | 1.27^{*} | |
[0.71; 1.81] | ||
(education)8 | 1.33^{*} | |
[0.79; 1.86] | ||
(income)2 | 0.40^{*} | |
[0.28; 0.52] | ||
(income)3 | 0.66^{*} | |
[0.52; 0.79] | ||
(income)4 | 0.91^{*} | |
[0.77; 1.06] | ||
(income)5 | 1.10^{*} | |
[0.94; 1.25] | ||
^{*} 0 outside 95% credible interval |
# source predicted probabilities function
source('https://raw.githubusercontent.com/jkarreth/JKmisc/master/MCMC_simcase_probs.R')
# extract post-warmup beta draws from all four chains
vote_betas <- extract(vote_fit, pars = 'beta')$beta
# calculate predicted probabilities based on median of observed variables
vote_pp <- MCMC_simcase_probs(model_matrix = X, mcmc_out = vote_betas, x_col = 2,
x_range_vec = seq(min(X[, 2]), max(X[, 2]), length.out = 100),
link = 'logit',
lower = .025, upper = .975)
# source clean ggplot theme
source('https://raw.githubusercontent.com/jayrobwilliams/RWmisc/master/theme_rw.R')
# plot predicted probabilities, with density plot of observed values
ggplot(data = vote_pp, aes(x = predictor, y = median_pp, ymin = lower_pp, ymax = upper_pp)) +
geom_ribbon(alpha = .5) +
geom_line() +
geom_density(data = elections, aes(x = as.numeric(scale(age))), inherit.aes = F, color = NA,
fill = 'darkblue', alpha = .25) +
coord_cartesian(ylim = c(0,1)) +
labs(x = '(scaled) age', y = 'Pr(voting)') +
theme_rw()