FITTING DEEP MRP FAST

The key limitation in fitting MRP with interactions is the speed of estimation. Earlier research has shown that fitting a single deep MRP model can take multiple hours (e.g., Goplerud Reference Goplerud2022). This is because of the presence of high-dimensional integrals that traditional methods either numerically approximate or address using Bayesian methods.

Variational inference provides a different approach for fast estimation; the goal is to find the best approximating distribution to the posterior given some simplifying assumption—usually that blocks of parameters are independent (Grimmer Reference Grimmer2011). However, the accuracy of this approximation can depend heavily on the specific problem, and thus needs extensive testing to ensure its reliability. Goplerud (Reference Goplerud2022) derived a new general algorithm for binomial hierarchical models and conducted extensive explorations of its performance on the single dataset considered in Ghitza and Gelman (Reference Ghitza and Gelman2013). Those algorithms fit an extremely complex hierarchical model in around 1 minute—versus hour(s) for existing approaches. It demonstrated excellent performance by recovering posterior means on coefficients and predictions that closely aligned with the gold standard approach of Bayesian estimation.Footnote ² Appendix A of the Supplementary Material provides a full exposition of the variational algorithm.

To illustrate the extensions in this paper, I focus on a simplified MRP model: Equation 1 shows a hierarchical model without fixed effects and with a random intercept for state and a random intercept for race, where $ {y}_i $ is the (binary) response for observation i. The notation follows Gelman and Hill (Reference Gelman and Hill2006) where $ {\alpha}_{g[i]}^{\mathrm{state}} $ selects the random effect for the state g of which observation i is a member. Appendix A of the Supplementary Material shows the generalization to random slopes, arbitrary numbers of random effects, and fixed effects.

(1) $$ {\displaystyle \begin{array}{l}{\displaystyle \begin{array}{l}{y}_i\sim \mathrm{B}\mathrm{e}\mathrm{r}\mathrm{n}({p}_i);\hskip0.7em {p}_i=\frac{\mathrm{exp}({\psi}_i)}{1+\mathrm{exp}({\psi}_i)};\hskip0.7em {\psi}_i={\beta}_0+{\alpha}_{g[i]}^{\mathrm{st}\mathrm{a}\mathrm{t}\mathrm{e}}+{\alpha}_{g^{\prime }[i]}^{\mathrm{race}};\\ {}\hskip0em {\alpha}_g^{\mathrm{st}\mathrm{a}\mathrm{t}\mathrm{e}}\sim N(0,{\sigma}_{\mathrm{st}\mathrm{a}\mathrm{t}\mathrm{e}}^2);\hskip0.7em {\alpha}_{g^{\prime}}^{\mathrm{race}}\sim N(0,{\sigma}_{\mathrm{race}}^2);\hskip1em p({\beta}_0)\propto 1;\\ {}{\sigma}_j^2\sim {p}_0({\sigma}_j^2)\hskip00.7em \mathrm{f}\mathrm{o}\mathrm{r}\hskip1em j\in \{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e},\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\}.\end{array}}\end{array}} $$

The choice of prior on the variance of the random effect $ {p}_0({\sigma}_j^2) $ is a difficult task. Some inferential techniques assume a flat prior. A risk of this strategy is that point estimates of $ {\sigma}_j^2 $ could be degenerate and equal zero; this sets all random effects estimates equal to zero. This problem is rather common for the Laplace approximation in glmer (Chung et al. Reference Chung, Gelman, Rabe-Hesketh, Liu and Dorie2015). A proper prior prevents this problem and thus is preferable. An Inverse-Gamma prior is a popular choice, but it is difficult to calibrate the strength correctly (Gelman Reference Gelman2006). Figure 1 illustrates this by showing the prior used in Goplerud (Reference Goplerud2022) (“Inverse-Gamma”) against the Huang and Wand (Reference Huang and Wand2013) prior employed in this paper. The latter prior implies the popular half-t prior on the standard deviation $ {\sigma}_j $ (Gelman Reference Gelman2006), and it can be generalized to multidimensional random effects. In that case, it imposes half-t priors on each marginal standard deviation while maintaining (if desired) a uniform prior on the correlations. Appendix A.1 of the Supplementary Material provides more information on this prior such as its density.

Figure 1. Comparing Prior Density for Random Effect Standard Deviation $ {\sigma}_j $

Note: The dashed line shows the prior density on $ {\sigma}_j $ given an Inverse-Gamma prior on $ {\sigma}_j^2 $ with $ {\alpha}_0=1 $ and $ {\beta}_0=1/2 $ . The solid line shows a Huang–Wand prior on $ {\sigma}_j^2 $ with $ \nu =2 $ and $ A=5 $ (i.e., half-t on $ {\sigma}_j $ ).

Figure 1 shows that Goplerud’s (Reference Goplerud2022) prior puts effectively no mass on small values of $ {\sigma}_j $ (e.g., $ P({\sigma}_j\le 0.25)\approx 0.0003 $ ). Thus, in the event where the true value is small (i.e., the random effect is mostly irrelevant), the prior results in too large estimates of $ {\sigma}_j $ , thereby underregularizing the coefficients, which likely results in poorer performance. By contrast, the Huang–Wand prior puts nontrivial weight on very small $ {\sigma}_j $ and thus allows for strong regularization when appropriate. Appendix A.1 of the Supplementary Material provides a stylized example of this phenomenon.

Unfortunately, Appendix A.5 of the Supplementary Material illustrates that naively incorporating the Huang–Wand prior dramatically increases estimation time. While it does increase the time per iteration, the major problem is that estimation requires 5–10 times more iterations to converge. Thus, a key contribution of this paper is to accelerate variational algorithms when this more appropriate prior is employed. Appendices A.2 and A.3 of the Supplementary Material provide, respectively, full explanations of the techniques employed: (i) a squared iterative method and (ii) a novel application of parameter expansion.

The model in Equation 1 can be extended by adding many interactions between geographic and demographic factors (“deep MRP”; Ghitza and Gelman Reference Ghitza and Gelman2013). However, Broniecki, Leemann, and Wüest (Reference Broniecki, Leemann and Wüest2022) note that additional state-level predictors (e.g., unemployment rate) may also provide considerable benefits. Unlike hierarchical models, many machine learning methods can automatically estimate nonlinear effects or interactions between these continuous predictors, whereas they must be specified explicitly for MRP.

I address that scenario by allowing estimation of nonlinear effects using splines as in a generalized additive model. Appendix A.4 of the Supplementary Material demonstrates how splines can be represented as additional hierarchical terms and thus estimated using the same variational algorithms.

Appendix B of the Supplementary Material provides simulations to illustrate the importance of using hierarchical models that include interactions or nonlinear effects. It shows that ignoring important interactions or nonlinearities hurts the performance of hierarchical models vis-à-vis alternative models, especially as the sample size increases. However, after those terms are included, hierarchical models perform well even against machine learning methods.

COMPARING METHODS FOR FITTING MRP

To compare deep hierarchical models against machine learning systematically, I use Buttice and Highton’s (Reference Buttice and Highton2013) popular dataset for validating new methods for MRP (e.g., Bisbee Reference Bisbee2019; Broniecki, Leemann, and Wüest Reference Broniecki, Leemann and Wüest2022; Ornstein Reference Ornstein2020). It consists of 89 policy questions that are collected from multiple years of the National Annenberg Election Studies (2000, 2004, and 2008) and the Cooperative Congressional Election Studies (2006 and 2008). The benefit of these large samples is that it is possible to use the entire dataset to get a “ground truth” by taking the observed average in each state while drawing a smaller subsample (e.g., 1,500 respondents) to mimic the conditions under which a researcher would need to apply MRP to obtain reliable state-level estimates.

Existing comparisons, however, only rely on a simple hierarchical model outlined below (Equation 2), following the original specification in Buttice and Highton (Reference Buttice and Highton2013). The model includes random effects for age, education (educ), gender–race combination (gXr), state, and region. The state-level continuous predictors pvote (state-level Republican presidential two-party vote share) and relig (share of population identifying as Evangelical Protestant or Mormon) are indexed with $ g[i] $ as they are constant within a state.

(2) $$ {\displaystyle \begin{array}{l}\mathrm{Pr}\left({y}_i=1\right)={\mathrm{logit}}^{-1}\left(\begin{array}{l}{\beta}_0+{\beta}_{\mathrm{pvote}}\cdot {\mathrm{pvote}}_{g\left[i\right]}+{\beta}_{\mathrm{relig}}\cdot {\mathrm{relig}}_{g\left[i\right]}+\\ {}{\alpha}_{g\left[i\right]}^{\mathrm{age}}+{\alpha}_{g\left[i\right]}^{\mathrm{educ}}+{\alpha}_{g\left[i\right]}^{\mathrm{gXr}}+{\alpha}_{g\left[i\right]}^{\mathrm{state}}+{\alpha}_{g\left[i\right]}^{\mathrm{region}}\end{array}\right),\\ {}{\alpha}_g^j\sim N\left(0,{\sigma}_j^2\right)\hskip0.3em \mathrm{for}\hskip0.3em \mathrm{all}\hskip0.3em j\hskip0.3em \mathrm{and}\hskip0.3em g.\end{array}} $$

This model includes no interactions between variables or nonlinear effects on continuous predictors, and thus is likely insufficiently rich to capture the true underlying relationship. It is reasonable to suspect that a “properly specified” MRP model should include at least some interactions to be competitive with methods that can automatically learn interactions or nonlinearities.

I consider three expansions of this model’s hierarchical component. First, I consider a deep MRP where all two-way interactions between demographics and geography are included (e.g., age–education and age–state), as well as a triple interaction between the three demographic variables. Second, I add splines to capture possible nonlinear effects in the state-level continuous variables. A third model includes both extensions. Table 1 summarizes the specifications; Appendix F of the Supplementary Material provides a demonstration of how to fit these models in the accompanying software (Goplerud Reference Ornstein2023).Footnote ³

Table 1. Deep MRP Specifications

Note: The second column indicates how these two variables are included. It is either “Linear” (Equation 2) or “Splines” where a spline is used to allow for nonlinear effects for each variable. The third column indicates how the random effects on age, education, gender × race, state, and region are included. “Additive” refers to five random effects added together (Equation 2). “Interacted” refers to including the interactions noted in the main text alongside the additive terms.

It is important to stress that this paper tracks the existing analyses comparing machine learning and MRP as closely as possible. There are thus other specifications that likely improve upon Table 1, although I show that adding this set of interactions enables MRP to perform competitively against state-of-the-art machine learning techniques.

DEEP MRP IN AN ENSEMBLE

The first comparison explores whether deep MRP adds much benefit when used alongside a suite of machine learning methods. I begin by using a technique known as “stacking” that takes the predictions of many different methods and combines them into a single prediction known as an ensemble. Ornstein (Reference Ornstein2020) applied this method to MRP and found considerable gains over traditional methods. The method performs K-fold cross-validation to get an out-of-sample prediction for each observation in the survey using each constituent model. The out-of-sample predictions are combined to see which weighted average (convex combination) best predicts the outcome, and then these weights are used to combine predictions before post-stratification. It is often the case that the ensemble outperforms any single method (Broniecki, Leemann, and Wüest Reference Broniecki, Leemann and Wüest2022; Ornstein Reference Ornstein2020), although this is not guaranteed and can be empirically assessed by, for example, using a held-out dataset.

A useful property of ensembles is the ability to compare the weights given to the constituent models. The weights reflect both the performance of the method and its “distinctiveness” from the other methods in the ensemble. Using each survey in Buttice and Highton (Reference Buttice and Highton2013), I drew 10 different samples of varying sizes and estimated an ensemble using fivefold cross-validation with the models in Ornstein (Reference Ornstein2020) where I swapped the traditional (“Simple”) MRP model with the deep MRP model from Table 1. Figure 2 summarizes the weights given to each model, averaging across the surveys and simulations.

Figure 2. Weights Given to Models in Ensemble

Note: This figure shows the ensemble weights averaged across all surveys and 10 simulations per survey. Each panel reports the sample size of the survey. Ninety-five percent confidence intervals are shown. The first four methods are from Ornstein (Reference Ornstein2020): LASSO, k-Nearest Neighbors (KNN), Random Forest (Forest), Gradient Boosting Machine (GBM). The final method is “Deep” MRP from Table 1.

The results provide clear support for the importance of deep MRP in an ensemble: it is the highest weighted method when the sample size is over 1,500 and is given over 40% of the total weight when the sample size is 5,000 or higher. The performance of deep MRP is corroborated by the fact that, of the methods in the ensemble, it has the lowest cross-validated error on the survey data when $ N>1,500 $ . In terms of computational time, fitting this deep model on the full survey takes around 30 seconds for the largest sample size of 10,000 observations. Thus, deep MRP can be added to an ensemble with limited cost.

Appendix D of the Supplementary Material explores the trade-off between model complexity and sample size. It examines a larger ensemble that includes all four models from Table 1.Footnote ⁴ While corroborating Figure 2—MRP models collectively receive around 40%–50% of the weight—it shows an expected trade-off between traditional and deep MRP where traditional (noninteractive) methods are given decreasing weight as the sample size increases. This suggests that the ensemble upweights more complex methods as the amount of data increases. The spline-based methods receive relatively low weight, but this may be due to the limited variation in the continuous variables that are measured at the state level.Footnote ⁵ Appendix D of the Supplementary Material also shows that, in terms of raw performance, a well-designed ensemble usually beats any single constituent method. The five-model ensemble beats all of its constituent methods by more than 5% on at least one sample size considered.

MRP AND BART

One limitation of ensembles is appropriately quantifying uncertainty of the post-stratified estimates. This is challenging because it may be difficult to quantify the uncertainty of the estimates from the individual machine learning methods used in the ensemble.Footnote ⁶ It also requires careful work to interpret the effects of the included variables. Thus, researchers often seek to rely on a single model that can incorporate uncertainty and remains highly flexible. To that end, BART (Chipman, George, and McCulloch Reference Chipman, George and McCulloch2010) is an attractive choice. The method is related to popular tree-based methods such as “random forests,” but it is implemented in a Bayesian framework that allows for quantification of uncertainty. Bisbee (Reference Bisbee2019) applies BART to MRP and reports that it substantially outperforms (traditional) MRP. The magnitude of the improvement is large (e.g., around 20%–30% decrease in mean absolute error [MAE]). This motivates an initial question: does BART improve upon deep MRP?

After some preliminary exploration, I discovered an error in Bisbee’s (Reference Bisbee2019) replication archive. Appendix E of the Supplementary Material describes it in detail; see also a corrigendum to Bisbee (Reference Bisbee2019) (Goplerud and Bisbee Reference Goplerud and Bisbee2022). In brief, the error arbitrarily injected random noise into the MRP estimates at the prediction stage. When this is corrected, traditional MRP’s performance increases markedly and is only slightly beaten by BART.

Following the main analysis in Bisbee (Reference Bisbee2019), Figure 3 shows the predictive accuracy on the surveys in Buttice and Highton (Reference Buttice and Highton2013) for a sample with 1,500 observations.Footnote ⁷ For simplicity, I show only three methods: the traditional (“Simple”) MRP following Bisbee’s (Reference Bisbee2019) provided code, a corrected traditional MRP, and deep MRP estimated using variational inference (see Table 1).

Figure 3. Visualizing Performance: MRP versus BART

Note: Each plot compares the mean absolute error (MAE) of BART and MRP, averaged across two-hundred simulations per survey. The points below the 45-degree line indicate that MRP performs worse. The three methods shown are the traditional MRP and prediction method in Bisbee (Reference Bisbee2019) (“MRP (Bisbee Reference Bisbee2019)”), an MRP with the same specifications but a corrected prediction (“MRP”), and a deep MRP fit using variational inference (“Deep”).

Fixing the error shows a noticeably different story; rather than being clearly beaten by BART, traditional MRP looks visually similar to BART in terms of its error across surveys. Table 2 provides a more concise quantitative summary. It shows the percentage gap in MAE versus BART averaged across the 89 surveys: $ \left({\mathrm{MAE}}_k-{\mathrm{MAE}}_{\mathrm{BART}}\right)/{\mathrm{MAE}}_{\mathrm{BART}}\times 100 $ , where $ {\mathrm{MAE}}_{\mathrm{k}} $ indicates the MAE of the model k averaged across two-hundred simulations. A positive number indicates that BART outperforms the other method. This measure is relative as BART and MRP both decrease the observed MAE as sample size increases.

Table 2. Relative Mean Absolute Error versus BART

Note: This table reports percentage gap in mean absolute error between BART and the alternative methods; positive numbers indicate that BART outperforms its competitor. Figure 3 defines the abbreviations.

Table 2 shows that BART outperforms traditional MRP, but it does so by quite small margins (1%–4%) and its relative advantage declines as sample size increases. Deep MRP performs slightly better versus BART; for modest sample sizes (3,000–6,000), it actually slightly outperforms BART, although it does slightly worse at small and very large sample sizes. In terms of performance, across all sample sizes, deep MRP outperforms BART between 45% and 55% of the time and, thus, they can be considered to reach an effective “draw” in terms of performance. The table also suggests a small-but-systematic improvement of deep MRP over traditional MRP. Comparing the traditional (“Simple”) MRP against deep MRP shows that, except for the largest sample sizes, traditional MRP performs around 1%–3% worse than deep MRP in terms of MAE and is beaten around 60%–70% of the time.