Hierarchical Multivariate Logistic Regression

Problem Set-up

The data set for this problem comes from the chapter on Modeling for accounting in data collection process, where Gelman et al talk about accounting the data collection process into the model set-up. A stratified random sampling opinion poll was conducted and data collected from 16 strata across the country.

There are two parameters of interest for this analysis -

1. Proportion of people voting bush (of the people who had an opinion) - \[ \alpha \ _1j \]
   
2. Proportion of people that gave an opinion - \[ \alpha \ _2j \]

Both of these can be derived from the above data. In addition, we know the total number of people samples (n = 1447), so we use the sample_proportion to calculate how many people were sampled in each strata. (we’ll later use this to convert the data into a long form format)

Sample Survey

library(tidyverse)
library(brms)
library(tidybayes)

Prior Predictive Simulation

Generalized Linear Models are interaction models even if you haven’t specified any interactions. Priors that may make sense on the logit scale don’t always make sense on the probability scale. brms by default uses student_t(3, 0, 10) priors for the population level effects in hierarchical models. These are called “weakly informative priors” and more often than not help with efficient sampling of the posterior space.
Even with small amount data, the likelihood would take over and these “weakly informative priors” won’t affect our inferences, but as you start fitting more complex models, putting sensible informative priors on your parameters is more likely to lead to better sampling.

Here we’ll look at what happens if we simulate data from the default priors and have a look at the plausible outcome space before observing the data

m1.prior <- fit1 <- brm(mvbind(a_1_binary, a_2_binary) ~ (1|p|strata), data = df_long, 
                  family = bernoulli, sample_prior = "only", 
                  chains = 4, cores = 2 ) 
saveRDS(m1.prior, here::here("static", "files", "m1_prior.rds")) 

Before observing the data, what is our prior belief about people voting for bush on average across the whole country?

WOWZA! The default student_t(3, 0, 10) priors puts almost all mass in either extremes, which says, we either expect a large majority to vote for bush or a large majority not to vote for bush. That seems like a pretty unreasonable prior belief to keep.

One of the criticisms that people have of bayesian statistics is that it is “subjective” and an often touted counter point is, to just use non-informative priors. Below is a classic example of why using non-informative priors does more harm than good.

posterior_samples(m1.prior) %>%
  select(1) %>%
  mutate(prior_bush = inv_logit_scaled(b_a1binary_Intercept)) %>%
  ggplot(.,aes(prior_bush)) + geom_density() + 
  labs(title = "Prior Probability of Voting for Bush", 
       x = "Probability",y = "Density")

We could do a similar prior simulation for each region, and unsurprisingly, they have a similar distribution -

tidybayes::add_linpred_draws(newdata = df_long[1, 'strata'], m1.prior) %>%
  filter(.category == "a1binary") %>%
  ggplot(.,aes(.value)) + geom_density()  

Correlation - Before observing the data, what is our prior belief about the correlation between voting bush and having an opinion?

It’s quite uniformly distributed.

posterior_samples(m1.prior) %>%
  select(cor_strata__a1binary_Intercept__a2binary_Intercept) %>%
  ggplot(.,aes(cor_strata__a1binary_Intercept__a2binary_Intercept)) + 
  geom_density() 

m2.prior <- brm(mvbind(a_1_binary, a_2_binary) ~ (1|p|strata), data = df_long, 
                  family = bernoulli, sample_prior = "only", 
                  chains = 4, cores = 2, 
                  prior = c(prior(student_t(3, 0, 1), class = sd, resp = a1binary), 
                            prior(student_t(3, 0, 1), class = sd, resp = a2binary),
                  prior(student_t(3, 0, 1), class = Intercept),
                  prior(lkj(2), class = cor))) 
saveRDS(m2.prior, here::here("static", "files", "m2_prior.rds")) 

posterior_samples(m2.prior) %>%
  select(1) %>%
  mutate(prior_bush = inv_logit_scaled(b_a1binary_Intercept)) %>%
  ggplot(.,aes(prior_bush)) + geom_density() + labs(title = "Prior Probability of Voting for Bush", 
                                                    x = "Probability",
                                                    y = "Density")

tidybayes::add_linpred_draws(newdata = df_long[1, 'strata'], m2.prior) %>%
  filter(.category == "a1binary") %>%
  ggplot(.,aes(.value)) + geom_density() + labs(title = "Probability of Voting Bush in Northeast, I")  

posterior_samples(m2.prior) %>%
  select(cor_strata__a1binary_Intercept__a2binary_Intercept) %>%
  ggplot(.,aes(cor_strata__a1binary_Intercept__a2binary_Intercept)) + 
  geom_density() + labs(title = "Prior Distribution of Correlation") 

m3.prior <- fit1 <- brm(mvbind(a_1_binary, a_2_binary) ~ (1|p|strata), data = df_long, 
                  family = bernoulli, sample_prior = "only", 
                  chains = 4, cores = 2, 
                  prior = c(prior(student_t(3, 0, 0.5), class = sd, resp = a1binary), 
                            prior(student_t(3, 0, 0.5), class = sd, resp = a2binary),
                  prior(student_t(3, 0, 1), class = Intercept),
                  prior(lkj(2), class = cor))) 
saveRDS(m3.prior, here::here("static", "files", "m3_prior.rds")) 

posterior_samples(m3.prior) %>%
  select(1) %>%
  mutate(prior_bush = inv_logit_scaled(b_a1binary_Intercept)) %>%
  ggplot(.,aes(prior_bush)) + geom_density() + 
  labs(title = "Prior Probability of Voting for Bush", 
       x = "Probability", y = "Density")

tidybayes::add_linpred_draws(newdata = df_long[1, 'strata'], m3.prior) %>%
  filter(.category == "a1binary") %>%
  ggplot(.,aes(.value)) + geom_density() + 
  labs(title = "Probability of Voting Bush in Northeast, I")

m4.prior <- brm(mvbind(a_1_binary, a_2_binary) ~ (1|p|strata), data = df_long, 
                  family = bernoulli, sample_prior = "only", 
                  chains = 4, cores = 2, 
                  prior = c(prior(student_t(3, 0, 0.5), class = sd, resp = a1binary), 
                            prior(student_t(3, 0, 0.5), class = sd, resp = a2binary),
                  prior(student_t(3, 0, 1), class = Intercept, resp = a1binary),
                  prior(student_t(3, 2, 1), class = Intercept, resp = a2binary),
                  prior(lkj(2), class = cor))) 
saveRDS(m4.prior, here::here("static", "files", "m4_prior.rds")) 

posterior_samples(m4.prior) %>%
  select(b_a2binary_Intercept) %>%
  mutate(prior_bush = inv_logit_scaled(b_a2binary_Intercept)) %>%
  ggplot(.,aes(prior_bush)) + geom_density() + 
  labs(title = "Prior Probability of Giving an Opinion on Average", 
       x = "Probability",y = "Density")

tidybayes::add_linpred_draws(newdata = df_long[6, 'strata'], m4.prior) %>%
  filter(.category == "a1binary") %>%
  ggplot(.,aes(.value)) + geom_density() + 
  labs(title = "Probability of Giving an Opinion in Northeast, I")  

Posterior Distribution

Let’s go ahead with the priors set in model m4.prior - one of the key elements of bayesian analysis is sensitivity analysis. We can change our priors and see how inferences in the posterior distribution change.

## posterior 
m1.posterior <- brm(mvbind(a_1_binary, a_2_binary) ~ (1|p|strata), data = df_long, 
                  family = bernoulli, 
                  chains = 4, cores = 4, 
                  prior = c(prior(student_t(3, 0, 0.5), class = sd, resp = a1binary), 
                            prior(student_t(3, 0, 0.5), class = sd, resp = a2binary),
                  prior(student_t(3, 0, 1), class = Intercept, resp = a1binary),
                  prior(student_t(3, 2, 1), class = Intercept, resp = a2binary),
                  prior(lkj(2), class = cor))) 
saveRDS(m1.posterior, here::here("static", "files", "m1_posterior.rds")) 

Posterior Checks

pp_check(m1.posterior, nsamples = 100, resp = "a1binary")

pp_check(m1.posterior, nsamples = 100, resp = "a2binary")  

new_d <- df_long %>%
  select(strata) %>%
  distinct()

tidybayes::add_linpred_draws(newdata = new_d, m1.posterior, seed = 123) %>%
  median_qi() %>%
  ggplot(., aes(.value, strata)) + geom_pointintervalh() + 
  facet_wrap(vars(group = .category))  

Shrinkage

Model checking doesn’t stop at the posterior predictive checks, one should ALWAYS be skeptical of the model results even if you’re satisfied at that stage. Another valuable tool is to simply compare the raw data estimates with the posterior predictions, and note where the results are different and where they are the same.

Hierarchical Models “shrink” estimates towards the grand mean, unless there’s strong evidence not to. That in fact is one of the major advantages, as it protects us from concluding groups are different, when in fact they’re not (multiple comparisons discussed in earlier post).

We see that quite clearly in the plot below - for both estimates, extreme values are shrunk towards the center of the population distribution.

An apple doesn’t fall too far away from the tree, so unless we have strong evidence that it is a peach, we treat it as being from the apple tree closest to it.

tidybayes::add_linpred_draws(newdata = new_d, m1.posterior, seed = 123) %>%
  mean_qi() %>%
  filter(.category == "a1binary") %>%
  inner_join(df_raw, by = c("strata" = "stratum")) %>%
  ggplot(., aes(.value, strata)) + geom_pointintervalh() + 
  geom_point(aes(proportion_bush), color = "red") + labs(title = "Estimated Proportion Voting for Bush in Each Region")

tidybayes::add_linpred_draws(newdata = new_d, m1.posterior, seed = 123) %>%
  mean_qi() %>%
  filter(.category == "a2binary") %>%
  inner_join(df_raw, by = c("strata" = "stratum")) %>%
  ggplot(., aes(.value, strata)) + geom_pointintervalh() + 
  geom_point(aes(1 -  proportion_no_opinion), color = "red") + labs(title = "Estimated Proportion having an opinion")

The northeast, I region has a lot more shrinkage than the others, which is a cause for concern. At this point, it is imperative to run some checks to ensure the validity of the model. And this will be our topic for the next post! How to check your model fit by looking at the implied posterior predictions.