Use the simulation approach to generating predicted probabilities to obtain 1000 simulated predicted probabilities for all demographic-state combinations of supporting marriage equality. Use these simulations to generate point estimates and 95% confidence intervals for the proportion of people within each state supporting marriage equality. Use at least 10 simulations.
Hint: the simulate()
function is your friend because it has a method for objects of class merMod
, which allows it to simulate new responses from the fitted model object – but be sure to account for the random effects in the model. You want to carry out the MRP process multiple times, substituting in a different simulated outcome vector each time. You can accomplish this either with a loop, or through a series of calls to apply()
with anonymous functions.
library(foreign) # load .dta files
library(lme4) # multilevel models
# read in megapoll
marriage_data <- read.dta('http://www.princeton.edu/~jkastell/MRP_primer/gay_marriage_megapoll.dta', convert.underscore = TRUE)
# read in state-level dataset
Statelevel <- read.dta('http://www.princeton.edu/~jkastell/MRP_primer/state_level_update.dta', convert.underscore = TRUE)
Statelevel <- Statelevel[order(Statelevel$sstate.initnum), ]
# read in sensus data
Census <- read.dta('http://www.princeton.edu/~jkastell/MRP_primer/poststratification%202000.dta', convert.underscore = TRUE)
Census <- Census[order(Census$cstate), ]
Census$cstate.initnum <- match(Census$cstate, Statelevel$sstate)
# from 1 for white males to 6 for hispanic females
marriage_data$race.female <- (marriage_data$female *3) + marriage_data$race.wbh
# from 1 for 18-29 with low edu to 16 for 65+ with high edu
marriage_data$age.edu.cat <- 4 * (marriage_data$age.cat -1) + marriage_data$edu.cat
# proportion of evangelicals in respondent's state
marriage_data$p.evang.full <- Statelevel$p.evang[marriage_data$state.initnum]
# proportion of mormon's in respondent's state
marriage_data$p.mormon.full <-Statelevel$p.mormon[marriage_data$state.initnum]
# combined evangelical + mormom proportions
marriage_data$p.relig.full <- marriage_data$p.evang.full + marriage_data$p.mormon.full
# kerry's % of 2-party vote in respondent's state in 2004
marriage_data$p.kerry.full <- Statelevel$kerry.04[marriage_data$state.initnum]
# same coding as above
Census$crace.female <- (Census$cfemale *3) + Census$crace.WBH
Census$cage.edu.cat <- 4 * (Census$cage.cat -1) + Census$cedu.cat
Census$cp.evang.full<- Statelevel$p.evang[Census$cstate.initnum]
Census$cp.mormon.full <- Statelevel$p.mormon[Census$cstate.initnum]
Census$cp.relig.full <- Census$cp.evang.full + Census$cp.mormon.full
Census$cp.kerry.full <- Statelevel$kerry.04[Census$cstate.initnum]
individual_model <- glmer(formula = yes.of.all ~ (1 | race.female) + (1 | age.cat)
+ ( 1 |edu.cat) + (1 | age.edu.cat) + (1 | state) + (1 | region)
+ (1 | poll) + p.relig.full + p.kerry.full, data = marriage_data,
family = binomial(link = 'logit'))
library(arm) # inverse logit function
library(doParallel) # parallelize simulation
library(ggplot2) # plots
library(plotly) # interactive plots
# parallel computing setup
registerDoParallel(makeCluster(parallel::detectCores()))
# simulate responses from model
sim <- simulate(individual_model, nsim = 50, re.form = NULL)
# create dataframe to merge with simulated outcome vectors; account for listwise deletion
sim_dat <- na.omit(marriage_data[, c('race.female', 'age.cat', 'edu.cat', 'age.edu.cat',
'state', 'region', 'poll', 'p.relig.full', 'p.kerry.full')])
# check that simulated outcome vectors have same length as dataframe
stopifnot(nrow(sim_dat) == nrow(sim))
# parallel loop for simulations
sim_outcome <- foreach(i = 1:ncol(sim), .packages = c('lme4', 'arm'), .combine = cbind) %dopar% {
# get ith simulated outcome vector
temp_dat <- data.frame(sim_dat, yes_sim = sim[, i])
# run model on simulated outcome vector
temp_mod <- glmer(yes_sim ~ (1 | race.female) + (1 | age.cat) + ( 1 |edu.cat)
+ (1 | age.edu.cat) + (1 | state) + (1 | region) + (1 | poll)
+ p.relig.full + p.kerry.full, data = temp_dat,
family = binomial(link = 'logit'))
# empty vector to hold state random effects
state_ranefs <- array(NA, c(51, 1))
# set state names as row names
dimnames(state_ranefs) <- list(c(Statelevel$sstate), 'effect')
# assign state random effects to array while preserving NAs
for (i in Statelevel$sstate) {
state_ranefs[i, ] <- ranef(temp_mod)$state[i, 1]
}
# set states with missing REs (b/c not in data) to zero
state_ranefs[, 1][is.na(state_ranefs[, 1])] <- 0
# prediction for each cell in census data
temp_pred <- invlogit(fixef(temp_mod)["(Intercept)"]
+ ranef(temp_mod)$race.female[Census$crace.female, 1]
+ ranef(temp_mod)$age.cat[Census$cage.cat, 1]
+ ranef(temp_mod)$edu.cat[Census$cedu.cat, 1]
+ ranef(temp_mod)$age.edu.cat[Census$cage.edu.cat, 1]
+ state_ranefs[Census$cstate, 1]
+ ranef(temp_mod)$region[Census$cregion, 1]
+ (fixef(individual_model)["p.relig.full"] * Census$cp.relig.full)
+ (fixef(temp_mod)["p.kerry.full"] * Census$cp.kerry.full))
# weight prediction by cell frequency
temp_wgt <- temp_pred * Census$cpercent.state
# sum proportion of cells supporting in each state and convert to percent
100 * as.vector(tapply(temp_wgt, Census$cstate, sum))
}
# collect mean, .025 and .975 quantiles, and proportion religious
sim_gg <- data.frame(State = Statelevel$sstate,
Lower = apply(sim_outcome, 1, quantile, probs = .025),
Mean = apply(sim_outcome, 1, mean),
Upper = apply(sim_outcome, 1, quantile, probs = .975),
Religion = Statelevel$p.evang + Statelevel$p.mormon)
# present estimates
sim_gg[, -ncol(sim_gg)]
# plot estimates with uncertainty, colored by religion, label to avoid reorder() in plotly
sim_plot <- ggplot(sim_gg, aes(x = reorder(State, Mean), y = Mean, ymin = Lower,
ymax = Upper, color = Religion, label = State)) +
geom_linerange(alpha = .75, size = 1) +
geom_point() +
labs(x = 'State', y = 'Estimated Support') +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 90),
legend.position = 'right',
plot.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank())
# interactive plot; ggplotly needs colour, not color
ggplotly(sim_plot, tooltip = c('label', 'ymax', 'y', 'ymin', 'colour'))