02-mrp_noncensus.Rmd

# MRP with Noncensus Variables

<style type="text/css">
h1.title {
  font-size: 32px;
  text-align: center;
}

h2 {
  padding-bottom: 4px;
}

h4.author {
  padding-top: 22px;
  text-align: center;
  font-style: italic;
}
h4.date {
  padding-top: 14px;
  font-size: 14px;
  text-align: center;
  font-style: italic;
  padding-bottom: 20px;
}
</style>

<script>
function myFunction() {
    var x = document.getElementById("myDIV");
    if (x.style.display === "none") {
        x.style.display = "block";
    } else {
        x.style.display = "none";
    }
}
function myFunction2() {
    var x = document.getElementById("myDIV2");
    if (x.style.display === "none") {
        x.style.display = "block";
    } else {
        x.style.display = "none";
    }
}
function myFunction3() {
    var x = document.getElementById("myDIV3");
    if (x.style.display === "none") {
        x.style.display = "block";
    } else {
        x.style.display = "none";
    }
}
</script>
<style>
#myDIV {
  width: 100%;
  padding: 20px 30px;
  background-color: rgba(192,192,192,0.15);
  margin-top: 10px;
  border-radius: 4px;
}

#myButton{
border-color: #008CBA;
background-color: rgba(192,192,192,0.05);
color: #008CBA;
border-radius: 4px;
}

#myDIV2 {
  width: 100%;
  padding: 20px 30px;
  background-color: rgba(192,192,192,0.15);
  margin-top: 10px;
  border-radius: 4px;
}

#myButton2{
border-color: #008CBA;
background-color: rgba(192,192,192,0.05);
color: #008CBA;
border-radius: 4px;
}

#myDIV3 {
  width: 100%;
  padding: 20px 30px;
  background-color: rgba(192,192,192,0.15);
  margin-top: 10px;
  border-radius: 4px;
}

#myButton3{
border-color: #008CBA;
background-color: rgba(192,192,192,0.05);
color: #008CBA;
border-radius: 4px;
}
</style>

```{r packages-2, message=FALSE, echo = FALSE}
library(brms)
library(rstanarm)
library(data.table)
library(dplyr)
library(forcats)
library(tidyr)
library(reshape2)
library(kableExtra)

library(ggplot2)
library(bayesplot)
library(gridExtra)
library(ggalt)
library(scales)
library(usmap)

theme_set(bayesplot::theme_default())

Sys.setenv(LOCAL_CPPFLAGS = '-march=corei7 -mtune=corei7')
options(mc.cores = parallel::detectCores(logical = FALSE))
```

When our sample population is different than our target population, MRP can only adjust for the predictors included in the model. As these are restricted by the variables in the poststratification table, which in turn are limited by the questions asked in the census, the characteristics that we can use for poststratification are quite reduced. This is the reason researchers tend to use simple demographic and geographic variables, which unfortunately do not provide much help if the bias in the survey originates from non-response in voters of a certain party, for instance. As a potential solution, @kastellec2015polarizing propose extending the postratification table using a survey that contains one or multiple non-census variables that could help adjusting for the differences between the sample and the target population. For instance, if our survey asked for partisanship, we could use the CCES to extend the poststratification table such as that it also contains this variable. The extension is done in two steps. First, we fit a multilevel model in which we try to predict partisanship in the CCES based on the variables available in the census. Second, we use this model to predict, for each cell in the original poststratification table, what proportion of subjects are Democrats, Republicans, or Independents. This extended poststratification table that contains partisanship will allow us to (a) generate MRP estimates that adjust for differential party nonresponse in the original survey; and/or (b) obtain estimates outcome of interest by party.

For this case study we will continue using the previous example of studying support for the right of employers to exclude abortion coverage.

```{r, echo = FALSE}
state_abb <- datasets::state.abb
state_fips <- c(1,2,4,5,6,8,9,10,12,13,15,16,17,18,19,20,21,22,23,24,
                25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,
                44,45,46,47,48,49,50,51,53,54,55,56)
recode_fips <- function(column) {
    factor(column, levels = state_fips, labels = state_abb)
}

get_se_bernoulli <- function(p, n){
  return(sqrt(p*(1-p)/n))
}

# Reed CCES data again, but this time we also read the partisanship variable
clean_cces2 <- function(df, remove_nas = TRUE){
  
  ## Abortion -- dichotomous (0 - Oppose / 1 - Support)
  df$abortion <- abs(df$CC18_321d-2)
  
  ## State -- factor
  df$state <- recode_fips(df$inputstate)
  
  ## Gender -- dichotomous (-0.5 Female, +0.5 Male)
  df$male <- abs(df$gender-2)-0.5
  
  ## ethnicity -- factor
  df$eth <- factor(df$race,
                   levels = 1:8,
                   labels = c("White", "Black", "Hispanic", "Asian", "Native American", "Mixed", "Other", "Middle Eastern"))
  df$eth <- fct_collapse(df$eth, "Other" = c("Asian", "Other", "Middle Eastern", "Mixed", "Native American"))
  
  ## Age -- cut into factor
  df$age <- 2018 - df$birthyr
  df$age <- cut(as.integer(df$age), breaks = c(0, 29, 39, 49, 59, 69, 120), 
                labels = c("18-29","30-39","40-49","50-59","60-69","70+"),
                ordered_result = TRUE)
  
  ## Education -- factor
  df$educ <- factor(as.integer(df$educ), 
                    levels = 1:6, 
                    labels = c("No HS", "HS", "Some college", "Associates", "4-Year College", "Post-grad"), ordered = TRUE)
  df$educ <- fct_collapse(df$educ, "Some college" = c("Some college", "Associates"))  
  
  # Party
  df$party <- factor(df$pid3,
                    levels = 1:5,
                    labels = c("Democrat", "Republican", "Independent", "Other", "Not sure"), 
                    ordered = TRUE)
  df$party <- fct_collapse(df$party, "Independent" = c("Independent", "Other", "Not sure"))
  
  # Clean and remove NAs
  df <- df %>% select(abortion, state, eth, male, age, educ, party)
    if (remove_nas){
    df <- df %>% drop_na()
  }
  
  return(df)
  
}
```

```{r, results = 'asis', cache=FALSE, warning=FALSE, message=FALSE}
# Read CCES data with the same outcome variable and predictors, but also including 
# party. This is done by the clean_cces2 function (not shown, see Github code)
cces_all_df <- read_csv("data_public/chapter1/data/cces18_common_vv.csv.gz")
cces_all_df <- clean_cces2(cces_all_df, remove_nas = TRUE)

# Read poststratification table
poststrat_df <- read_csv("data_public/chapter1/data/poststrat_df.csv")

# Read state-level predictors and add them to the CCES and poststratification table
statelevel_predictors_df <- read_csv('data_public/chapter1/data/statelevel_predictors.csv')
cces_all_df <- left_join(cces_all_df, statelevel_predictors_df, 
                    by = "state", keep = FALSE)
poststrat_df <- left_join(poststrat_df, statelevel_predictors_df, 
                          by = "state", keep = FALSE)
```

## Model-based Extension of the Poststratification Table 

As we have described, we start fitting a multilevel model to predict partisanship as a function of the same demographic and geographic variables used in the standard MRP model, which will allow us to predict the proportion of Republicans, Democrats, and Independents in each row of the poststratification table. As there are three levels for partisanship, we use a Bayesian multinomial (i.e. unordered) logistic regression which can be fitted in `brms` (currently, `rstanarm` does not support multinomial logistic regression).

For this extension step we should use a survey that we think is to some degree representative with respect to the variable that we are trying to include in the poststratification table. In our example, if we extended our census-based poststratification table using a highly non-representative survey with respect to party, we would indeed generated a biased poststratification table and ultimately obtain compromised MRP estimates. In other words, this is our opportunity to bring outside information in order to generate a richer poststratification table that can adjust for potential biases in the main survey, so we need to make sure that the survey we use to extend the poststratification table is trustworthy with respect to the non-census variable. 

In this example, we will use a 5,000-person sample of the CCES to extend the poststratification table to include partisanship, which is addressed in the CCES:

> Generally speaking, do you think of yourself as a ...? (Democrat, Republican, Independent, Other, Not Sure)

For simplicity, we included the few respondents that indicated "Other" or "Not Sure" as Independents.

```{r, cache=FALSE, warning=FALSE, message=FALSE, eval=FALSE}
# Setting seed to arbitrary number for reproducibility
set.seed(1010)

# Taking random sample from the CCES survey
cces_unbiased_df <- cces_all_df %>% sample_n(5000)

fit_party <- brm(party ~ (1 | state) + (1 | eth) + (1 | age) + (1 | educ) + male +
                   (1 | male:eth) + (1 | educ:age) + (1 | educ:eth) +
                   repvote + factor(region),
   family = "categorical",
   data = cces_unbiased_df,
   prior = c(prior(normal(0, 5), class = Intercept),
             prior(normal(0, 1), class = b),
             prior(exponential(0.5), class = sd, dpar = muIndependent),
             prior(exponential(0.5), class = sd, dpar = muRepublican)),
   control = list(adapt_delta = 0.9, max_treedepth  = 10),
   seed = 1010)
```

```{r, echo = FALSE, eval=TRUE}
# we save the model for future use. By default we do not retrain the model and
# save it, only retrieving the previously file version. To train the model again,
# simply change eval=TRUE in the previous cell and eval=FALSE in this one.
#saveRDS(fit_party, file = "data_public/chapter2/models/fit_party_example.rds")
fit_party <- readRDS("data_public/chapter2/models/fit_party_example.rds")
```

This model gives us, for each poststratification cell $j$, an estimate for the proportion of Democrats ($\hat{\theta}_{{\rm Democrat}, j}$), Republicans ($\hat{\theta}_{{\rm Republican}, j}$), and Independents ($\hat{\theta}_{{\rm Independent}, j}$). We can multiply these quantities by the number of people in cell $j$ to estimate the number of Democrats ($N_j \: \hat{\theta}_{{\rm Democrat}, j}$), Republicans ($N_j \: \hat{\theta}_{{\rm Republican}, j}$), and Independents ($N_j \: \hat{\theta}_{{\rm Independent}, j}$), obtaining an extended poststratification table in which each cell has been expanded into three. That is, if the original poststratification table had $J$ rows (e.g. 12,000 in our case), the new one will have $3 J$ (e.g. 36,000). There is, however, a certain complication that must be taken into account. The model-based estimates for the proportion of Democrats, Republicans, and Independents are not single numbers, but several draws from the posterior distribution that capture the uncertainty about these estimates. For instance, if we have 500 draws for $\hat{\theta}_{{\rm Democrat}, j}$, $\hat{\theta}_{{\rm Republican}, j}$, and $\hat{\theta}_{{\rm Independent}, j}$, we can imagine 500 poststratification tables with different numbers for each cell.

```{r}
# Use posterior_epred to predict partisanship for original poststratification table
pred_mat <- brms::posterior_epred(fit_party, newdata = poststrat_df, ndraws = 500, transform = TRUE)

# Extend poststratification table
poststrat_df_threefold <- poststrat_df[rep(seq_len(nrow(poststrat_df)), each = 3), ]
poststrat_df_threefold$party <- rep(c("Democrat", "Republican", "Independent"), nrow(poststrat_df))

# Calculate new numbers for the cells of the new poststratification table. K 
# is a matrix containing 36000 rows (one for each cell of the poststratification table) 
# and 500 columns (corresponding to the 500 draws).
K_theta <- apply(pred_mat, 1, function(x){as.vector(t(x))})
K <- K_theta * rep(poststrat_df$n, each = 3)
```

<button id="myButton" onclick="myFunction()" >Show code explanation </button>
<div id="myDIV" style="display:none">

More generally, let's say that for each $\hat{\theta}_{{\rm Democrat}, j}$, $\hat{\theta}_{{\rm Republican}, j}$, and $\hat{\theta}_{{\rm Independent}, j}$ we have a vector that contains $D_1$ draws from the posterior distribution. In the code $D_1 = 500$, as we have randomly selected 500 draws from the `pred_mat` matrix. Therefore, `pp_expect` (the `brms` equivalent to `posterior_epred`) will return an array with $J$ (number of rows of the original poststratification table) $\times D_1$ (number of draws from posterior distribution) $\times 3$ (estimates for Democrats, Republicans, and Independents) dimensions.

The first step is to take this $J \times D_1 \times 3$ array and transform it to a $3J \times D_1$ matrix. Then, we multiply each row of the matrix by the corresponding $N_j$ of the original poststratification table. This gives us a new $3J \times D_1$ matrix that we will call $K$, and which contains the $D_1$ draws for each of the $3J$ cells of the new poststratification table.

</div>

<br>

In sum, we started with a poststratification table with 12,000 rows. Here we can see the first three rows:

```{r, echo=FALSE}
poststrat_df[1:3,1:6] %>% kable() %>% kable_styling(full_width = TRUE)
```

We have used a model-based approach to include partisanship in this poststratification table, that now has 36,000 rows (again, each row in the original table has been split into three). However, in order to consider the uncertainty in these model-based estimates we have actually built 500 different poststratification tables. Here we show the first 9 rows of one of these 500 poststratification tables:

```{r, echo=FALSE}
poststrat_df_threefold[1:9,c(1, 2, 3, 4, 5, 9, 6)] %>% mutate(n = round(K[1:9, 1], 1)) %>% kable() %>% kable_styling(full_width = TRUE)
```

## Adjusting for Nonresponse Bias

We have described how to extend the poststratification table by including partisanship. Now, we will use this poststratification table to adjust for differential party nonresponse.

### Setting up example with an artificially nonrepresentative sample

To demostraty how non-census MRP can adjust for party, we will use a survey that is biased with respect to party. As we are already familiar with the CCES dataset, what we are going to do is to take a different sample of 5,000 respondents that simulates a high nonresponse rate among Republicans and, to a lesser degree, Independents.

```{r}
# Random sample of 5,000 that weights by party
cces_biased_df <- cces_all_df %>% sample_n(5000, weight = I((cces_all_df$party=="Democrat")*1 + 
                                                        (cces_all_df$party=="Independent")*0.75 + 
                                                        (cces_all_df$party=="Republican")*0.5))
```

Previously, we saw that the national average support for requiring companies to cover abortion in their insurance plans was around `r round(mean(cces_all_df$abortion, na.rm = TRUE), 3)*100`% according to the CCES. Comparatively, this biased sample of the CCES gives an estimate of `r round(mean(cces_biased_df$abortion, na.rm = TRUE), 3)*100`%. This is not surprising, as missing Republicans and Independents in the survey should reduce support for the employers' right to decline abortion coverage.

### Standard MRP

We fit a standard MRP (i.e. without including party) on the nonrepresentative sample, using the same model as in the MRP introduction and the non-extended poststratification table.

```{r, cache=FALSE, warning=FALSE, message=FALSE, eval=FALSE}
fit_abortion_standard <- stan_glmer(abortion ~ (1 | state) + (1 | eth) + (1 | age) + (1 | educ) + male +
                      (1 | male:eth) + (1 | educ:age) + (1 | educ:eth) +
                      repvote + factor(region),
    family = binomial(link = "logit"),
    data = cces_biased_df,
    prior = normal(0, 1, autoscale = TRUE),
    prior_covariance = decov(scale = 0.50),
    adapt_delta = 0.99,
    seed = 1010)
```

```{r, echo = FALSE, eval=TRUE}
# we save the model for future use. By default we do not retrain the model and
# save it, only retrieving the previously file version. To train the model again,
# simply change eval=TRUE in the previous cell and eval=FALSE in this one.
#saveRDS(fit_abortion_standard, file = "data_public/chapter2/models/fit_abortion_standard.rds")
fit_abortion_standard <- readRDS("data_public/chapter2/models/fit_abortion_standard.rds")
```

```{r}
standard_epred_mat <- rstanarm::posterior_epred(fit_abortion_standard, newdata = poststrat_df, draws = 500)
standard_mrp_estimates_vector <- (standard_epred_mat %*% poststrat_df$n)/sum(poststrat_df$n)
```

```{r, echo = FALSE}
cat("Standard MRP estimate mean, sd: ", round(mean(standard_mrp_estimates_vector), 3), round(sd(standard_mrp_estimates_vector), 3))
```

The standard MRP with the nonrepresentative sample gives a national-level estimate of `r 100*round(mean(standard_mrp_estimates_vector), 3)`% ($\pm$ `r 100*round(sd(standard_mrp_estimates_vector), 3)`%). As this estimate does not consider partisanship, standard MRP is not being able to adjust for the smaller statement support that results from oversampling Democrats.

### Non-census MRP with partisanship as a predictor

In the first section we have created a poststratification table that contains partisanship. After doing this, the next step of the non-census MRP approach is to fit the same model as we did in the standard MRP, but also including party as a predictor:

$$
Pr(y_i = 1) = logit^{-1}(
\alpha_{\rm s[i]}^{\rm state}
+ \alpha_{\rm a[i]}^{\rm age}
+ \alpha_{\rm r[i]}^{\rm eth}
+ \alpha_{\rm e[i]}^{\rm educ}
+ \beta^{\rm male} \cdot {\rm Male}_{\rm i} 
+ \alpha_{\rm g[i], r[i]}^{\rm male.eth}
+ \alpha_{\rm e[i], a[i]}^{\rm educ.age}
+ \alpha_{\rm e[i], r[i]}^{\rm educ.eth}
+ \alpha_{\rm p[i]}^{\rm party}
)
$$

$$
\begin{align*}
\alpha_{\rm s}^{\rm state} &\sim {\rm Normal}(\gamma^0 + \gamma^{\rm south} \cdot {\rm South}_{\rm s} + \gamma^{\rm midwest} \cdot {\rm Midwest}_{\rm s} + \gamma^{\rm west} \cdot {\rm West}_{\rm s} + \gamma^{\rm repvote} \cdot {\rm RepVote}_{\rm s}, \sigma_{\rm state}) \textrm{ for s = 1,...,50}\\
\alpha_{\rm a}^{\rm age} & \sim {\rm Normal}(0,\sigma_{\rm age}) \textrm{ for a = 1,...,6}\\
\alpha_{\rm r}^{\rm eth} & \sim {\rm Normal}(0,\sigma_{\rm eth}) \textrm{ for r = 1,...,4}\\
\alpha_{\rm e}^{\rm educ} & \sim {\rm Normal}(0,\sigma_{\rm educ}) \textrm{ for e = 1,...,5}\\
\alpha_{\rm g,r}^{\rm male.eth} & \sim {\rm Normal}(0,\sigma_{\rm male.eth}) \textrm{ for g = 1,2 and r = 1,...,4}\\
\alpha_{\rm e,a}^{\rm educ.age} & \sim {\rm Normal}(0,\sigma_{\rm educ.age}) \textrm{ for e = 1,...,5 and a = 1,...,6}\\
\alpha_{\rm e,r}^{\rm educ.eth} & \sim {\rm Normal}(0,\sigma_{\rm educ.eth}) \textrm{ for e = 1,...,5 and r = 1,...,4}\\
\alpha_{\rm p}^{\rm party} & \sim {\rm Normal}(0,\sigma_{\rm party}) \textrm{ for p = 1,2,3}\\
\end{align*}
$$

```{r, cache=FALSE, warning=FALSE, message=FALSE, eval=FALSE}
fit_abortion_noncensus <- stan_glmer(abortion ~ (1 | state) + (1 | eth) + (1 | age) + (1 | educ) + male +
                     (1 | male:eth) + (1 | educ:age) + (1 | educ:eth) +
                     repvote + factor(region) + (1 | party),
   family = binomial(link = "logit"),
   data = cces_biased_df,
   prior = normal(0, 1, autoscale = TRUE),
   prior_covariance = decov(scale = 0.50),
   adapt_delta = 0.99,
   seed = 1010)
```

```{r, echo = FALSE, eval=TRUE}
# we save the model for future use. By default we do not retrain the model and
# save it, only retrieving the previously file version. To train the model again,
# simply change eval=TRUE in the previous cell and eval=FALSE in this one.
#saveRDS(fit_abortion_noncensus, file = "data_public/chapter2/models/fit_abortion_noncensus.rds")
fit_abortion_noncensus <- readRDS("data_public/chapter2/models/fit_abortion_noncensus.rds")
```

Using `posterior_epred` allows us to estimate abortion coverage support for each of the cells in the extended poststratification table. As we set `draws = 500`, we obtain 500 estimates for each cell. In standard MRP, we will weight each the statement support estimates for each poststratification cell by the number of people in that cell according to the model-based estimates obtained in the previous section. However, as in this case the number of people in each cell was estimated with uncertainty, we need to propagate the uncertainty in the first (party prediction) model to the final MRP estimates. Essentially, what we can do is randomly pick one of the 500 statement support estimates for each poststratification cell (i.e. a 36,000 vector) we have just obtained and weight it by one of the 500 poststratification tables that resulted from the first model. Repeating the process for the remaining draws gives us a distribution of 500 MRP estimates for national support that correctly captures the uncertainty in the two models. 

```{r}
# Use posterior_epred to predict stance on abortion insurance coverage for extended poststratification table 
noncensus_epred_mat <- rstanarm::posterior_epred(fit_abortion_noncensus, newdata = poststrat_df_threefold, draws = 500)

# Calculate national MRP estimates propagating uncertainty from the two models
noncensus_mrp_estimates_vector <- colSums(t(noncensus_epred_mat)*K) / sum(K[,1])
```

```{r, echo = FALSE}
cat("Noncensus MRP estimate mean, sd: ", round(mean(noncensus_mrp_estimates_vector), 3), round(sd(noncensus_mrp_estimates_vector), 3))
```

<button id="myButton2" onclick="myFunction2()" >Show code explanation </button>
<div id="myDIV2" style="display:none">

More generally, `posterior_epred` returns a $D_2 \times 3J$ matrix we name $P$, where $D_2$ is the number of draws and $3J$ is, again, the number of cells in the extended poststratification table that includes partisanship.

With the two matrices $P$ a $K$, the poststratification step becomes:

$$\hat{\theta}_{MRP} = \frac{P \times K}{\sum_j^{3J} k_{j d}} {\rm \: where \:} K = (k_{jd})$$

Were $\hat{\theta}_{MRP}$ is a matrix containing $D_2 \times D_1$ estimates for the proportion of statement support.

In case of $D_1 = D_2$, we can do an elementwise operation such that $\hat{\theta}_{MRP} = \frac{\sum_j^{3J} (P^T \odot K)_{jd}}{\sum_j^{3J} k_{jd}}$ in order for the operation to be less computationally expensive, as we only combine each draw (i.e. column) in matrix $P$ with one draw (i.e. column) in matrix $K$ (instead with all the draws in matrix $K$). This more efficient method, which is the one we implement, returns a vector of length $D_1 = D_2$. The only consideration is that we must make sure that $D_1 = D_2$, which in this case required randomly picking 500 draws from $K$ and only obtaining another 500 draws for $P$ (by specifying `draws = 500`).

</div>

<br>

Our national-level estimate for the right to exclude abortion coverage from employer-sponsored insurance resulting from this non-census variable MRP is `r 100*round(mean(noncensus_mrp_estimates_vector), 3)`% (`r 100*round(sd(noncensus_mrp_estimates_vector), 3)`%). Unsurprisingly, this is much closer to the full (unbiased) 60,000 participant survey (`r 100*round(mean(cces_all_df$abortion), 3)` $\pm$ `r 100*round(sqrt(mean(cces_all_df$abortion)*(1-mean(cces_all_df$abortion))/nrow(cces_all_df)), 3)`%) than the standard MRP estimate seen above (`r 100*round(mean(standard_mrp_estimates_vector), 3)` $\pm$ `r 100*round(sd(standard_mrp_estimates_vector), 3)`%). Using an extended poststratification table that contained partisanship allowed us to adjust for differential partisan nonresponse.

Of course, we can also obtain state-level estimates and compare standard MRP with non-census MRP.

```{r}
# Create empty dataframe
state_estimates_df <- data.frame(
  state = state_abb,
  standard_mrp_state_estimate = NA,
  standard_mrp_state_estimate_se = NA,
  noncensus_mrp_state_estimate = NA,
  noncensus_mrp_state_estimate_se = NA,
  full_cces_state_estimate = NA,
  full_cces_state_estimate_se = NA,
  n_full = NA
)

# Loop to populate the dataframe
for(i in 1:nrow(state_estimates_df)) {
  # Filtering condition for standard_epred_mat (12,000 rows)
  filtering_condition <- which(poststrat_df$state == state_estimates_df$state[i])
  # Filtering condition for noncensus_epred_mat (36,000 rows)
  filtering_condition_threefold <- which(poststrat_df_threefold$state == state_estimates_df$state[i])
  
  # Standard MRP estimate
  state_standard_epred_mat <- standard_epred_mat[ ,filtering_condition]
  k_filtered <- poststrat_df[filtering_condition, ]$n
  standard_mrp_state_estimates_vector <- state_standard_epred_mat %*% k_filtered / sum(k_filtered)
  state_estimates_df$standard_mrp_state_estimate[i] <- mean(standard_mrp_state_estimates_vector)
  state_estimates_df$standard_mrp_state_estimate_se[i] <- sd(standard_mrp_state_estimates_vector)
  
  # Noncensus MRP estimate 
  state_noncensus_epred_mat <- noncensus_epred_mat[ ,filtering_condition_threefold]
  K_filtered <- K[filtering_condition_threefold, ]
  noncensus_mrp_state_estimates_vector <- colSums(t(state_noncensus_epred_mat)*K_filtered) / colSums(K_filtered)
  state_estimates_df$noncensus_mrp_state_estimate[i] <- mean(noncensus_mrp_state_estimates_vector)
  state_estimates_df$noncensus_mrp_state_estimate_se[i] <- sd(noncensus_mrp_state_estimates_vector)
  
  # Full survey estimate
  state_estimates_df$full_cces_state_estimate[i] <- mean(filter(cces_all_df, state==state_estimates_df$state[i])$abortion)
  state_estimates_df$n_full[i] <- nrow(filter(cces_all_df, state==state_estimates_df$state[i]))
  state_estimates_df$full_cces_state_estimate_se[i] <- get_se_bernoulli(state_estimates_df$full_cces_state_estimate[i], state_estimates_df$n_full[i])
}
```

```{r, fig.width=10, fig.height=3.5, warning=FALSE, message=FALSE, results = 'hide', echo=FALSE}
states_order <- poststrat_df %>% group_by(state) %>% summarise(repvote = first(repvote)) %>% arrange(repvote) %>% .$state
state_estimates_df$state <- factor(state_estimates_df$state, levels = states_order, ordered = TRUE)

compare1 <- ggplot(data=state_estimates_df) +
  geom_point(aes(x=state, y=standard_mrp_state_estimate), color = "#E37B1C") +
  geom_errorbar(aes(ymin=standard_mrp_state_estimate - 2*standard_mrp_state_estimate_se,
                    ymax=standard_mrp_state_estimate + 2*standard_mrp_state_estimate_se,
                    x=state), alpha=.5, width = 0, color = "#E37B1C") +
  geom_point(data=state_estimates_df, aes(x=state, y=noncensus_mrp_state_estimate), color = "#7B1CE3") +
  geom_errorbar(data=state_estimates_df, aes(ymin=noncensus_mrp_state_estimate - 2*noncensus_mrp_state_estimate_se, 
                                             ymax=noncensus_mrp_state_estimate + 2*noncensus_mrp_state_estimate_se, 
                                             x=state), alpha=.5, width = 0, color = "#7B1CE3") +
  geom_point(aes(x=state, y=full_cces_state_estimate), color = "#1CE37B") +
  geom_errorbar(data=state_estimates_df, aes(ymin=full_cces_state_estimate - 2*full_cces_state_estimate_se, 
                                             ymax=full_cces_state_estimate + 2*full_cces_state_estimate_se, 
                                             x=state), alpha=.5, width = 0, color = "#1CE37B") +
  scale_y_continuous(breaks=c(0,.25,.5,.75,1), 
                     labels=c("0%","25%","50%","75%","100%"), 
                     expand=c(0,0))+
  coord_cartesian(ylim=c(0, 1)) +
  theme_bw()+
  labs(x="States",y="Support")+
  theme(legend.position="none",
        axis.title=element_text(size=10),
        axis.text.y=element_text(size=10),
        axis.text.x=element_text(angle=90,size=8, vjust=0.3),
        legend.title=element_text(size=10),
        legend.text=element_text(size=10))

compare2 <- ggplot(data = state_estimates_df) +
  geom_point(aes(y=mean(standard_mrp_estimates_vector), x = .25), color = "#E37B1C") +
  geom_errorbar(data=state_estimates_df, aes(y = mean(standard_mrp_estimates_vector), 
                x = .25,
                ymin = mean(standard_mrp_estimates_vector) - 2*sd(standard_mrp_estimates_vector),
                ymax = mean(standard_mrp_estimates_vector) + 2*sd(standard_mrp_estimates_vector)),
                width = 0, color = "#E37B1C") +
  geom_text(data = data.frame(), aes(x = Inf, y = mean(standard_mrp_estimates_vector) - 0.03, label = "Standard MRP"), 
            hjust = -.05, size = 4, color = "#E37B1C") +
  geom_point(aes(y = mean(noncensus_mrp_estimates_vector), x = .75), color = "#7B1CE3") +
  geom_errorbar(aes(y = mean(noncensus_mrp_estimates_vector), 
                x = .75, 
                ymin = mean(noncensus_mrp_estimates_vector) - 2*sd(noncensus_mrp_estimates_vector),
                ymax = mean(noncensus_mrp_estimates_vector) + 2*sd(noncensus_mrp_estimates_vector)),
                width = 0, color = "#7B1CE3") +
  geom_text(data = data.frame(), aes(x = Inf, y = mean(noncensus_mrp_estimates_vector) - 0.01, label = "Non-census MRP"), 
            hjust = -.05, size = 4, color = "#7B1CE3") +
  scale_y_continuous(breaks=c(0,.25,.5,.75,1),
                     labels=c("0%","25%","50%","75%","100%"),
                     limits=c(0,1),expand=c(0,0)) +
  geom_point(data = data.frame(), aes(y=mean(cces_all_df$abortion), x = .5), color = "#1CE37B") +
  geom_errorbar(data = data.frame(), aes(y = mean(cces_all_df$abortion), 
                x = .5,
                ymin = mean(cces_all_df$abortion) - 2*sqrt(mean(cces_all_df$abortion)*(1-mean(cces_all_df$abortion))/nrow(cces_all_df)),
                ymax = mean(cces_all_df$abortion) + 2*sqrt(mean(cces_all_df$abortion)*(1-mean(cces_all_df$abortion))/nrow(cces_all_df))),
                width = 0, color = "#1CE37B") +
  geom_text(data = data.frame(), aes(x = Inf, y = mean(cces_all_df$abortion)+0.05, label = "Complete Survey"), 
            hjust = -.06, size = 4, color = "#1CE37B") +
  scale_y_continuous(breaks=c(0,.25,.5,.75,1),
                     labels=c("0%","25%","50%","75%","100%"),
                     limits=c(0,1),expand=c(0,0))+
  scale_x_continuous(limits=c(0,1),expand=c(0,0), breaks=c(.25, .75)) +
  coord_cartesian(clip = 'off') +
  theme_bw() +
  labs(x="Population",y="")+
   theme(legend.position="none",
        axis.title.y=element_blank(),
        axis.title.x=element_text(size=10, margin = margin(t = 19, r = 0, b = , l = 0)),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        legend.title=element_text(size=10),
        legend.text=element_text(size=10),
        plot.margin = margin(5.5, 105, 5.5, 5.5, "pt")
        )

bayesplot_grid(compare1,compare2, 
               grid_args = list(nrow=1, widths = c(5,1.4)))
```

In general, we see that the estimates from the standard MRP are upwardly biased with respect to the 60,000 survey estimates. Conversely, the MRP with non-census variables is able to adjust for the differential partisan nonresponse.

## Obtaining Estimates for Non-census Variable Subgroups

Even if we do not suspect that our survey population is different from our target population with respect to a non-census variable, using non-census MRP can allow us to obtain different estimates for the levels of the non-census variable. Here, we obtain and plot support for declining coverage of abortions by state and party within state.

```{r, warning=FALSE}
subgroup_estimates_df <- cces_biased_df %>% expand(state, party) %>%
  mutate(noncensus_mrp_subgroup_estimate = NA,
         noncensus_mrp_subgroup_estimate_se = NA)

for(i in 1:nrow(subgroup_estimates_df)) {
  filtering_condition_threefold <- which(poststrat_df_threefold$state == subgroup_estimates_df$state[i] &
                                           poststrat_df_threefold$party == subgroup_estimates_df$party[i])
  subgroup_noncensus_epred_mat <- noncensus_epred_mat[ ,filtering_condition_threefold]
  K_filtered <- K[filtering_condition_threefold, ]
  noncensus_mrp_subgroup_estimates_vector <- colSums(t(subgroup_noncensus_epred_mat)*K_filtered) / colSums(K_filtered)
  
  subgroup_estimates_df$noncensus_mrp_subgroup_estimate[i] <- mean(noncensus_mrp_subgroup_estimates_vector)
  subgroup_estimates_df$noncensus_mrp_subgroup_estimate_se[i] <- sd(noncensus_mrp_subgroup_estimates_vector)
}
```

```{r, message=FALSE, warning=FALSE, echo=FALSE, cache=FALSE, fig.height=3.25, fig.width=12, fig.align = "center"}
states_map <- us_map(regions = "states")
subgroup_estimates_df_melted <- subgroup_estimates_df %>% select(state, noncensus_mrp_subgroup_estimate, party)
states_map  <- left_join(states_map, subgroup_estimates_df_melted, by = c("abbr" = "state")) %>% drop_na()

ggplot(states_map, aes(x = x, y = y, group = group)) +
  geom_polygon(colour = "lightgray") +
  geom_polygon(aes(fill = noncensus_mrp_subgroup_estimate)) + theme_void() + facet_grid(cols = vars(party)) + 
  scale_fill_gradient2(midpoint = 0.5, limits = c(0, 1), breaks = c(0, .5, 1),
                       name = "Support", low = muted("blue"), high = muted("red")) + 
  theme(legend.margin=margin(l = 0.5, unit='cm'))

```