  Non-linear hierarchical models are commonly used in many disciplines.
However, inference in the presence of non-nested effects and on large datasets
is challenging and computationally burdensome. This paper provides two
contributions to scalable and accurate inference. First, I derive a new
mean-field variational algorithm for estimating binomial logistic hierarchical
models with an arbitrary number of non-nested random effects. Second, I propose
"marginally augmented variational Bayes" (MAVB) that further improves the
initial approximation through a step of Bayesian post-processing. I prove that
MAVB provides a guaranteed improvement in the approximation quality at low
computational cost and induces dependencies that were assumed away by the
initial factorization assumptions.
  I apply these techniques to a study of voter behavior using a
high-dimensional application of the popular approach of multilevel regression
and post-stratification (MRP). Existing estimation took hours whereas the
algorithms proposed run in minutes. The posterior means are well-recovered even
under strong factorization assumptions. Applying MAVB further improves the
approximation by partially correcting the under-estimated variance. The
proposed methodology is implemented in an open source software package.
