Regression and Other Stories: Poststratification

Poststratification after estimation. See Chapter 17 in Regression and Other Stories.

The CBS News poll conducted from 12--16 October 2016 reported that, among likely voters who preferred one of the two major-party candidates, 45% intended to vote for Donald Trump and 55% for Hillary Clinton. Of these respondents, Party ID 33% Republican, 40% Republican, 27% independent.

source: http://www.cbsnews.com/news/cbs-poll-clintons-lead-over-trump-widens-with-three-weeks-to-go/ and https://www.scribd.com/document/327938789/CBS-News-Poll-10-17-toplines

Effective sample size of likely voters

254 Republican, 282 Democrat, 242 Independent

Compare to:

exit polls 2012  32 38 29
exit polls 2016  33 36 31

Republicans:  77% Trump,  8% Clinton (must normalize to 100%)
Democrats:     5% Trump, 89% Clinton (must normalize to 100%)
Independents: 36% Trump, 38% Clinton (must normalize to 100%)

---

Load packages

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

Simulate fake data

n_pid <- c(254, 282, 242)
n <- sum(n_pid)
pid_names <- c("Republican", "Democrat", "Independent")
pid <- rep(pid_names, n_pid)
n_vote <- as.list(rep(NA, 3))
n_vote[[1]] <- round(c(0.77, 0.08)*n_pid[1])
n_vote[[2]] <- round(c(0.05, 0.89)*n_pid[2])
n_vote[[3]] <- round(c(0.36, 0.38)*n_pid[3])
vote <- NULL
y_bar_cells <- rep(NA, 3)
for (j in 1:3){
  n_vote[[j]]<- c(n_vote[[j]], n_pid[j] - sum(n_vote[[j]]))
  vote <- c(vote, rep(c(1, 0, NA), n_vote[[j]]))
  y_bar_cells[j] <- mean(vote[pid==pid_names[j]], na.rm=TRUE)
  round(y_bar_cells[j], 3)
}
poll <- data.frame(vote, pid)
# write.csv(poll, root("Poststrat/data","poll.csv"), row.names=FALSE)
# poll <- read.csv(root("Poststrat/data","poll.csv"))
head(poll)

  vote        pid
1    1 Republican
2    1 Republican
3    1 Republican
4    1 Republican
5    1 Republican
6    1 Republican

summary(poll)

      vote               pid     
 Min.   :0.00   Democrat   :282  
 1st Qu.:0.00   Independent:242  
 Median :0.00   Republican :254  
 Mean   :0.45                    
 3rd Qu.:1.00                    
 Max.   :1.00                    
 NA's   :118                     

Simple poststrat

poststrat_data <- data.frame(pid=c("Republican", "Democrat", "Independent"),
                             N=c(0.33, 0.36, 0.31))
round(sum(poststrat_data$N * y_bar_cells), 3)

[1] 0.469

Linear model

Raw estimate

round(mean(poll$vote, na.rm=TRUE), 3)

[1] 0.45

stan_glm

fit_1 <- stan_glm(vote ~ factor(pid), data = poll, refresh = 0)
print(fit_1, digits=2)

stan_glm
 family:       gaussian [identity]
 formula:      vote ~ factor(pid)
 observations: 660
 predictors:   3
------
                       Median MAD_SD
(Intercept)            0.05   0.02  
factor(pid)Independent 0.43   0.03  
factor(pid)Republican  0.85   0.03  

Auxiliary parameter(s):
      Median MAD_SD
sigma 0.34   0.01  

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

Poststrat using posterior_linpred()

epred_1 <- posterior_epred(fit_1, newdata=poststrat_data)
poststrat_est_1 <- epred_1 %*% poststrat_data$N/sum(poststrat_data$N)
print(c(mean(poststrat_est_1), mad(poststrat_est_1)), digits=2)

[1] 0.469 0.013

Add extra uncertainty

n_sim <- nrow(epred_1)
poststrat_est_2 <- poststrat_est_1 + rnorm(n_sim, 0, 0.02)
print(c(mean(poststrat_est_2), mad(poststrat_est_2)), digits=2)

[1] 0.469 0.025

Logistic model

Fit the regression

fit <- stan_glm(vote ~ factor(pid), family=binomial(link="logit"), data = poll, refresh = 0)
print(fit, digits=2)

stan_glm
 family:       binomial [logit]
 formula:      vote ~ factor(pid)
 observations: 660
 predictors:   3
------
                       Median MAD_SD
(Intercept)            -2.88   0.28 
factor(pid)Independent  2.83   0.31 
factor(pid)Republican   5.17   0.37 

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

Raw estimate

round(mean(poll$vote, na.rm=TRUE), 3)

[1] 0.45

Poststrat using the predict function

X_population <- data.frame(pid=c("Republican", "Democrat", "Independent"))
N_population <- c(0.33, 0.36, 0.31)
predict_poststrat <- colMeans(posterior_epred(fit, newdata=X_population))
poststrat_est_2 <- sum(N_population*predict_poststrat)/sum(N_population)
round(poststrat_est_2, digits=3)

[1] 0.469

Just to compare, poststrat using the original data

predict_a <- predict(fit, type="response")
round(mean(predict_a), 3)

[1] 0.45

This doesn't work--it just spits back the raw estimate--because it's not using the external population info which is what makes poststrat work.