2.1 Notation and Estimands

We consider a finite population of N individuals indexed $i=1,\ldots ,N$ . We observe the outcome $Y_i$ for units that respond to the survey. Define a binary variable $R_i$ that denotes inclusion in the survey, where $R_i = 1$ indicates that unit i responds; $n = \sum _i R_i$ is the total number of respondents. This response variable can include respondents to a probability sample or a convenience sample. These response variables $R_i$ are the only random component of our setup, and all expectations, variances, and probability statements will be with respect to the randomness in the response variables. In addition, each individual is also associated with a set of d categorical covariates $X_{i1},\ldots ,X_{id}$ , where the $\ell $ th covariate is categorical with $J_\ell $ levels, so that the vector of covariates is $\boldsymbol {X}_i \in [J_1] \times \cdots \times [J_d]$ .Footnote ⁴

Rather than consider these variables individually, we rewrite the vector $X_i$ as a single categorical covariate, the cell for unit i, $S_i \in [J]$ , where $J \leq J_1\times \cdots \times J_d$ is the number of unique levels in the target population.Footnote ⁵ While we primarily consider a fixed population size N and total number of distinct cells J, in Section C of the Supplementary Material, we extend this setup to an asymptotic framework where both the population size and the number of cells can grow. With these cells, we can summarize the covariate information. We denote $\boldsymbol {N}^{\mathcal {P}} \in \mathbb {N}^J$ as the population count vector with $N^{\mathcal {P}}_s = \sum _i \mathbb {1}\{S_i = s\}$ , and $\boldsymbol {n}^{\mathcal {R}} \in \mathbb {N}^J$ as the respondent count vector, with $n_s^{\mathcal {R}} = \sum _i R_i \mathbb {1}\{S_i = s\}$ . While the population count $N^{\mathcal {P}}_s> 0$ for every cell s, the respondent count $n_s^{\mathcal {R}}$ may be equal to zero. We assume that we have access to these cell counts for both the respondent sample and the population.

Finally, for each cell s, we consider a set of binary vectors $\boldsymbol {D}_s^{(k)}$ that denote the cell in terms of its kth order interactions, and collect the vectors into matrices $\boldsymbol {D}^{(k)} = [\boldsymbol {D}_1^{(k)} \ldots \boldsymbol {D}_J^{(k)}]'$ , and into one combined $J \times J$ design matrix $\boldsymbol {D} = [\boldsymbol {D}^{(1)},\ldots , \boldsymbol {D}^{(d)}]$ . This is the usual design matrix for interaction terms in linear models, with the rows corresponding to unique cells, rather than units. It can be constructed using the model.matrix command in the R programming language. Figure A.1 in the Supplementary Material shows an example of $\boldsymbol {D}$ with three covariates from our running example: a binary for self-reported female, discretized age, and party identification.

Our goal is to estimate the population average outcome, which we can write as a cell-size weighted average of the within-cell averages, that is,

(1) $$ \begin{align} \mu \equiv \frac{1}{N}\sum_{i=1}^N Y_i = \sum_{s=1}^J \frac{N^{\mathcal{P}}_s}{N} \mu_s, \;\;\; \text{ where } \;\;\; \mu_s \equiv \frac{1}{N_s^{\mathcal{P}}}\sum_{S_{i} = s} Y_i. \end{align} $$

To estimate the population average, we rely on the average outcomes we observe within each cell. For cell s, the responder average is

(2) $$ \begin{align} \bar{Y}_s \equiv \frac{1}{n_s^{\mathcal{R}}}\sum_{S_i = s} R_i Y_i. \end{align} $$

We invoke the assumption of missingness at random within cells, so that the cell responder averages are unbiased for the true cell averages (Rubin Reference Rubin1976).

We denote the propensity score as $P(R_i = 1) \equiv \pi _i$ and the probability of responding conditional on being in cell s as $\pi (s) \equiv \frac {1}{N_s^{\mathcal {P}}}\sum _{S_{i} = s} \pi _i$ . For a probability sample, $\pi _i$ denotes the joint probability of both selection into the survey and responding. The analyst knows and controls the selection probabilities but does not know the probability of response given selection. For a convenience sample, $\pi _i$ is the unknown probability of inclusion in the sample. For both cases, the overall propensity score $\pi _i$ is unknown. We assume that this probability is nonzero for all cells.

These assumptions allow us to identify the overall population average using only the observed data. However, in order for Assumption 1 to be plausible, we will need the cells to be very fine-grained and have covariates that are quite predictive of nonresponse and the outcome. As we will see below, this creates a trade-off between identification and estimation: the former becomes more plausible with more fine-grained information, whereas the latter becomes more difficult (see also D’Amour et al. Reference D’Amour, Ding, Feller, Lei and Sekhon2020).

2.2 Review: Raking and Post-Stratification

We first consider estimating $\mu $ by taking a weighted average of the respondents’ outcomes, with weights $\hat {\gamma }_i$ for unit i. Because the cells are the most fine-grained information we have, we will restrict the weights to be constant within each cell. Researchers typically face a trade-off, discussed below, between feasibility of the cell-level estimator, known as post-stratification, and imposing modeling assumptions to address sparsity concerns, such as through calibration weighting. Our goal is to parsimoniously balance the feasibility advantage of calibration with the nonparametric bias-reduction of post-stratification. In particular, we leverage the fact that first-order coefficients typically explain more variation in outcome and response models than higher-order interactions (Cox Reference Cox1984), and so mitigating bias can be accomplished by enforcing balance on these first-order terms and allowing approximate balance on higher-order interactions.

We start by denoting the estimated weight for cell s as $\hat {\gamma }(s)$ and estimate the population average $\mu $ via

(3) $$ \begin{align} \hat{\mu}\left(\hat{\gamma}\right) \equiv \frac{1}{N}\sum_{i = 1}^N R_i \hat{\gamma}_i Y_i = \frac{1}{N}\sum_s n_s^{\mathcal{R}} \hat{\gamma}(s) \bar{Y}_s. \end{align} $$

If the individual probabilities of responding were known, we could choose to weight cell S by the inverse of the propensity score, $\frac {1}{\pi (s)}$ . Since the propensity score $\pi (s)$ is unknown, researchers can estimate $\pi (s)$ via $\hat {\pi }(s)$ and weight cell s by the inverse estimated propensity score $\hat {\gamma }(s) = \frac {1}{\hat {\pi }(s)}$ .

Post-stratification weights estimate the propensity score as the proportion of the population in cell s that responded, $\frac {n_s^{\mathcal {R}}}{N_s^{\mathcal {P}}}$ , leading to $\hat {\gamma }^{\text {ps}}(s) = \frac {N_s^{\mathcal {P}}}{n_s^{\mathcal {R}}}$ . Under Assumption 1, these post-stratification weights give an unbiased estimator for the population average $\mu $ . However, in practice post-stratification faces feasibility constraints due to sparsity—this estimator is only defined if there is at least one responder within each cell, which is unlikely with even a moderate number of covariates. For example, in Figure 1, cells corresponding to nearly one quarter of the population are empty when post-stratifying on only five of our eight covariates.Footnote ⁶

An alternative, calibration, chooses weights so that the weighted distribution of covariates exactly matches that of the full population, which identifies the population average under a modeling assumption called linear ignorability (Hartman, Hazlett, and Sterbenz Reference Hartman, Hazlett and Sterbenz2021). A common implementation, raking on margins, matches the marginal distribution by solving a convex optimization problem that finds the minimally “disperse” weights that satisfy this balance constraint (Deming and Stephan Reference Deming and Stephan1940; Deville and Särndal Reference Deville and Särndal1992; Deville, Särndal, and Sautory Reference Deville, Särndal and Sautory1993). For example, raking on the margins ensures that the percentage of female, the percentage of Democrats, Independents, and Republicans, and so forth for our eight covariates all match the target population; however, it would not ensure that the pairwise interactions are matched. Bias is mitigated so long as these higher-order interactions are not important to the outcome or response model. Calibration thus addresses feasibility concerns by focusing on balancing lower-order interactions, which are less likely to contain empty cells. For example, in Figure 1, focusing on three-way interactions and below would retain nearly the whole population.

Specifically, starting with baseline weights $q(s)$ for cell s, calibration finds weights that solve

(4) $$ \begin{align} \begin{aligned} \min_{\gamma} \;\; & \sum_{s} n_s^{\mathcal{R}} f(\gamma(s), q(s))\\ \text{subject to } & \sum_s \boldsymbol{D}_s^{(1)} n_s^{\mathcal{R}} \gamma(s) = \sum_s \boldsymbol{D}_s^{(1)} N_s^{\mathcal{P}}\\ & L \leq \gamma(s) \leq U \;\;\; \forall \;\; s = 1,\ldots,J, \end{aligned} \end{align} $$

where the function $f:\mathbb {R} \times \mathbb {R} \to \mathbb {R}^+$ , in the first equation, is a measure of dissimilarity. Some choices include the scaled and squared distance $f(\gamma (s), q(s)) = \frac {1}{2q(s)}(\gamma (s) - q(s))^2$ and the KL divergence $f(\gamma (s), q(s)) = \gamma (s) \log \frac {\gamma (s)}{q(s)}$ . Common choices for baseline weights, $q(s)$ , include probability weights when they are known (e.g., Deville and Särndal Reference Deville and Särndal1992), or estimated inverse propensity weights. Below we will use uniform weights $q(s) = \frac {N}{n}$ , for which the dissimilarity measures reduce to standard notions of the spread of the weights: the scaled and squared distance becomes the variance of the weights and the KL divergence becomes the entropy of the weights. As we will see below, the variance of the weights is directly related to the variance of the weighting estimate, and so is a natural choice of penalty.

The second equation encodes the calibration constraints, or the weighted survey moments that must exactly match the target population. In addition to the dispersion penalty in the objective, in the third equation, we constrain the weights to be between a lower bound $L = 0$ and an upper bound $U = \infty $ , which restricts the normalized cell weights, $\frac {1}{N}\gamma (s) n_s^{\mathcal {R}}$ to be in the $J-1$ simplex. This ensures that the imputed cell averages are in the convex hull of the respondents’ values and so do not extrapolate, and that the resulting estimator $\hat {\mu }(\hat {\gamma })$ is between the minimum and maximum sample outcomes.Footnote ⁷

Note that the general calibration procedure can also include exact constraints on higher-order interaction terms, and including all interactions recovers the post-stratification weights. We use the raking and post-stratification nomenclature to distinguish between calibration weighting with first-order margins and with all interactions. Several papers have proposed “soft” or “penalized” calibration approaches to relax the exact calibration constraint in Equation (4), allowing for approximate balance in some covariates (see, e.g., Rao and Singh Reference Rao and Singh1997; Park and Fuller Reference Park and Fuller2009; Guggemos and Tillé Reference Guggemos and Tillé2010; Gao, Yang, and Kim Reference Gao, Yang and Kim2023). Our multilevel calibration approach below can be seen as adapting the soft calibration approach to full post-stratification.

Before turning to our proposal for balancing higher-order interactions, we briefly describe some additional approaches. Chen, Valliant, and Elliott (Reference Chen, Valliant and Elliott2019) and McConville et al. (Reference McConville, Breidt, Lee and Moisen2017) discuss model-assisted calibration approaches which rely on the LASSO for variable selection. Caughey and Hartman (Reference Caughey and Hartman2017) select higher-order interactions to balance using the LASSO. Hartman et al. (Reference Hartman, Hazlett and Sterbenz2021) provide a kernel balancing method for matching joint covariate distributions between non-probability samples and a target population. Linzer (Reference Linzer2011) provides a latent class model for estimating cell probabilities and marginal effects in highly-stratified data.

Finally, in Section 4, we also discuss outcome modeling strategies as well as approaches that combine calibration weights and outcome modeling. A small number of papers have previously explored this combination for non-probability samples. Closest to our setup is Chen et al. (Reference Chen, Li and Wu2020), who combine inverse propensity score weights with a linear outcome model and show that the resulting estimator is doubly robust. Related examples include Yang, Kim, and Song (Reference Yang, Kim and Song2020), who give high-dimensional results for a related setting; and Si et al. (Reference Si, Trangucci, Gabry and Gelman2020), who combine weighting and outcome modeling in a Bayesian framework.

3.1 Estimation Error

We begin by inspecting the estimation error $\hat {\mu }(\gamma ) - \mu $ for weights $\gamma $ . Define the residual for unit i as $\varepsilon _i \equiv Y_i - \mu _{S_i}$ , and the average respondent residual in cell s as $\bar {\varepsilon }_s = \frac {1}{n_s^{\mathcal {R}}}\sum _{S_i = s}R_i \varepsilon _i$ . The estimation error decomposes into a term due to imbalance in the cell distributions and a term due to idiosyncratic variation within cells:

(5) $$ \begin{align} \hat{\mu}\left(\hat{\gamma}\right) - \mu = \underbrace{\frac{1}{N}\sum_{s}\left(n_s^{\mathcal{R}}\hat{\gamma}(s) - N_s^{\mathcal{P}}\right) \times \mu_s}_{\text{imbalance in cell distribution}} + \underbrace{\frac{1}{N}\sum_{s}n_s^{\mathcal{R}} \hat{\gamma}(s) \bar{\varepsilon}_s}_{\text{idiosyncratic error}}. \end{align} $$

By Assumption 1, which states that outcomes are missing at random within cells, the idiosyncratic error will be zero on average, and so the bias will be due to imbalance in the cell distribution. By Hölder’s inequality, we can see that the MSE, given the number of respondents in each cell, is

(6) $$ \begin{align} \begin{aligned} \mathbb{E}\left[\left(\hat{\mu}\left(\hat{\gamma}\right) - \mu\right) ^ 2 \mid n^{\mathcal{R}} \right] & = \underbrace{\frac{1}{N^2}\left(\sum_s \left(n_s^{\mathcal{R}} \hat{\gamma}(s) - N_s^{\mathcal{P}}\right)\mu_s\right)^2}_{\text{bias}^2} + \underbrace{\sum_{s} \left(\frac{n_s^{\mathcal{R}}}{N}\right)^2\hat{\gamma}(s)^2\sigma_s^2}_{\text{variance}}\\ & \leq \frac{1}{N^2}\sum_s\mu_s^2 \times \underbrace{\sum_s\left(n_s^{\mathcal{R}} \hat{\gamma}(s) - N_s^{\mathcal{P}}\right)^2}_{\text{imbalance in cell distribution}} + \underbrace{\sigma^2 \sum_{s} \left(\frac{n_s^{\mathcal{R}}}{N}\right)^2\hat{\gamma}(s)^2}_{\text{noise}}, \end{aligned} \end{align} $$

where $\sigma ^2_s = \text {Var}(\bar {Y}_s \mid n^{\mathcal {R}})$ and $\sigma ^2 = \max _s \sigma ^2_s$ . We therefore have two competing objectives if we want to control the MSE for any given realization of our survey. To minimize the bias, we want to find weights that control the imbalance between the true and weighted proportions within each cell. To minimize the variance, we want to find “diffuse” weights so that the sum of the squared weights is small.

The decomposition above holds for imbalance measures across all of the strata, without taking into account their multilevel structure. In practice, we expect cells that share features to have similar outcomes on average. We can therefore have finer-grained control by leveraging our representation of the cells into their first-order marginals $\boldsymbol {D}_s^{(1)}$ and interactions of order k, $\boldsymbol {D}_s^{(k)}$ . To do this, consider the infeasible population regression using the $\boldsymbol {D}_s^{(k)}$ representation as regressors,

(7) $$ \begin{align} \min_{\boldsymbol{\eta}} \sum_{i=1}^N \left(Y_i - \sum_{k=1}^d \boldsymbol{D}_{S_i}^{(k)} \cdot \boldsymbol{\eta}_k \right)^2. \end{align} $$

With the solution to this regression, $\boldsymbol {\eta }^\ast = (\boldsymbol {\eta }_1^\ast ,\ldots ,\boldsymbol {\eta }_d^\ast )$ , we can decompose the population average in cell s based on the interactions between the covariates, $\mu _s = \sum _{k=1}^d \boldsymbol {D}_s^{(k)} \cdot \boldsymbol {\eta }_k^\ast $ .Footnote ⁸

This decomposition in terms of the multilevel structure allows us to understand the role of imbalance in lower- and higher-order interactions on the bias. Plugging this decomposition into Equation (6), we see that the bias term in the conditional MSE further decomposes into the level of imbalance for the kth-order interactions weighted by the strength of the interactions:

(8) $$ \begin{align} \begin{aligned} \mathbb{E}\left[\left(\hat{\mu}\left(\hat{\gamma}\right) - \mu\right) ^ 2 \mid n^{\mathcal{R}} \right] & = \frac{1}{N^2}\left(\sum_{k=1}^d \boldsymbol{\eta}^{\ast}_k \cdot \sum_s \left(n_s^{\mathcal{R}} \hat{\gamma}(s) - N_s^{\mathcal{P}} \right)\boldsymbol{D}_s^{(k)}\right)^2 + \sum_{s} \left(\frac{n_s^{\mathcal{R}}}{N}\right)^2\hat{\gamma}(s)^2\sigma_s^2\\ & \leq \frac{1}{N^2}\left(\sum_{k=1}^d \left\|\boldsymbol{\eta}^{\ast}_k\right\|_2 \left\|\sum_s \left(n_s^{\mathcal{R}} \hat{\gamma}(s) - N_s\right)\boldsymbol{D}_s^{(k)}\right\|_2\right)^2 + \sigma^2\sum_{s} \left(\frac{n_s^{\mathcal{R}}}{N}\right)^2 \hat{\gamma}(s)^2. \end{aligned} \end{align} $$

Equation (8) formalizes the benefits of raking on the margins as in Equation (4). If there is an additive functional form with no influence from higher-order interaction terms—so $\boldsymbol {\eta }^\ast _k = 0$ for all $k \geq 2$ —then raking will yield an unbiased estimator. Even if the “main effects” are stronger than any of the interaction terms and so the coefficients on the first-order terms, $\|\boldsymbol {\eta }^\ast _1\|_2$ , are large relative to the coefficients for higher-order terms, raking can remove a large amount of the bias and so it is often seen as a good approximation (Mercer, Lau, and Kennedy Reference Mercer, Lau and Kennedy2018). However, ignoring higher-order interactions entirely can lead to bias. We therefore propose to find weights that prioritize main effects while still minimizing imbalance in interaction terms when feasible.

3.2 Multilevel Calibration

We now design a convex optimization problem that controls the conditional MSE on the right-hand side of Equation (8). To do this, we apply the ideas and approaches developed for approximate balancing weights (e.g., Hirshberg, Maleki, and Zubizarreta Reference Hirshberg, Maleki and Zubizarreta2019; Ning, Sida, and Imai Reference Ning, Sida and Imai2020; Wang and Zubizarreta Reference Wang and Zubizarreta2020; Xu and Yang Reference Xu and Yang2022; Zubizarreta Reference Zubizarreta2015) to the problem of controlling for higher-order interactions, using our MSE decomposition as a guide. We find weights that control the imbalance in all interactions in order to control the bias, while penalizing the sum of the squared weights to control the variance. Specifically, we solve the following optimization problem:

(9) $$ \begin{align} \begin{aligned} \min_{\gamma \in \mathbb{R}^J} \;\; & \underbrace{\sum_{k=2}^d \frac{1}{\lambda_k} \left\|\sum_s \boldsymbol{D}_s^{(k)} n_s^{\mathcal{R}} \gamma(s) - \boldsymbol{D}_s^{(k)} N_s^{\mathcal{P}} \right\|_2^2}_{\text{approximate higher-order balance}} + \underbrace{\sum_{s} n_s^{\mathcal{R}} \left(\gamma(s) - \frac{N}{n}\right)^2}_{\text{variance penalty}}\\[0.5em] \text{subject to} \;\; & \underbrace{\sum_s \boldsymbol{D}_s^{(1)} n_s^{\mathcal{R}} \gamma(s) = \sum_s \boldsymbol{D}_s^{(1)} N_s^{\mathcal{P}}}_{\text{raking constraint}}\\[0.5em] & L \leq \gamma(s) \leq U\;\;\; \forall s=1,\ldots J, \end{aligned} \end{align} $$

where the $\lambda _k$ are hyper-parameters. Note that the optimization problem is on the scale of the population counts rather than the population proportions. This follows the bound on the bias in Equation (8); categories and interactions with larger population counts $\boldsymbol {D}^{(k)'} \boldsymbol {N}^{\mathcal {P}}$ contribute more to the bias. With a common hyper-parameter, this means that higher-order interactions or smaller categories—which have lower population counts—will have less weight by design.

We can view this optimization problem as adding an additional objective optimizing for higher-order balance to the usual raking estimator in Equation (4) with the scaled squared dissimilarity function and uniform baseline weights. As with the raking estimator, the multilevel calibration weights are constrained to exactly balance first-order margins. Subject to this exact marginal constraint, the weights then minimize the imbalance in kth-order interactions for all $k=2,\ldots ,d$ . In this way, the multilevel calibration weights approximately post-stratify by optimizing for balance in higher-order interactions rather than requiring exact balance in all interactions as the post-stratification weights do. Following the bias–variance decomposition in Equations (6) and (8), this objective is penalized by the sum of the squared weights. As with the raking estimator, we can also replace the variance penalty with a different penalty function, including penalizing deviations from known sampling weights following Deville and Särndal (Reference Deville and Särndal1992). A benefit of the variance penalty is that Equation (9) is a quadratic program (QP) that can be solved efficiently via first-order methods (Boyd et al. Reference Boyd, Parikh, Chu, Peleato and Eckstein2010). We use OSQP, an efficient first-order method developed by Stellato et al. (Reference Stellato, Banjac, Goulart, Bemporad and Boyd2020).

Variations on the measure of approximate higher-order balance are also possible. Equation (8) uses the Cauchy–Schwarz inequality to relate the sum of the squared imbalances in kth-order interactions to the bias. We can instead use Hölder’s inequality to relate the maximal imbalance via the $L^\infty $ norm. Using this measure in Equation (9) would find weights that minimize the imbalance in the worst-balanced kth-order interaction, related to the proposal from Wang and Zubizarreta (Reference Wang and Zubizarreta2020). We could also solve a variant of Equation (9) without the multilevel structure encoded by the $\boldsymbol {D}_s^{(k)}$ variables. This would treat cells as entirely distinct and perform no aggregation across cells while approximately post-stratifying. From our discussion in Section 3.1, this would ignore the potential bias gains from directly leveraging the multilevel structure. Finally, in Section B of the Supplementary Material, we inspect the Lagrangian dual of this optimization problem and show that the weights are a form of propensity score weights with a multilevel GLM propensity score model.

An important component of the optimization problem are the hyper-parameters $\lambda _k$ for $k=2,\ldots ,d$ . They control the relative priority that balancing the higher-order interactions receives in the objective in an inverse relationship, and define a bias–variance trade-off. If $\lambda _k$ is large, then the weights will be more regularized and the kth-order interaction terms will be less prioritized. In the limit as all $\lambda _k \to \infty $ , no weight is placed on any interaction terms, and Equation (9) reduces to raking on the margins. Conversely, if $\lambda _k$ is small, more importance will be placed on balancing kth-order interactions. For example, if $\lambda _2 = 0$ , then the optimization problem will rake on margins and second-order interactions. As all $\lambda _k \to 0$ , we recover post-stratification weights, if they exist.

The trade-off between bias and variance can be viewed as one between balance and effective sample size. Improving the balance decreases the bias, but comes at the expense of increasing variance by decreasing the effective sample size. In practice, we suggest explicitly tracing out this trade-off. For a sequence of potential hyper-parameter values $\lambda ^{(1)},\lambda ^{(2)},\ldots $ , set all of the hyper-parameters to be $\lambda _k = \lambda ^{(j)}$ . We can then look at the two components of the objective in Equation (9), plotting the level of imbalance $\sum _{k=2}^d\left \|\sum _s \boldsymbol {D}_s^{(k)} n_s^{\mathcal {R}} \gamma (s) - \boldsymbol {D}_s^{(k)} N_s^{\mathcal {P}} \right \|_2^2$ against the effective sample size $n^{\text {eff}} = \frac {\left (\sum _s n_s^{\mathcal {R}} \hat {\gamma }(s)\right )^2}{\sum _s n_s^{\mathcal {R}} \hat {\gamma }(s)^2}$ . After fully understanding this trade-off, practitioners can choose a common $\lambda $ somewhere along the curve. For example, in our analysis in Section 5, we choose $\lambda $ to achieve 95% of the potential balance improvement in higher-order terms of $\lambda = 0$ relative to raking.

4.1 Using an Outcome Model for Bias Correction

A common alternative to the weighting approach above is to estimate the population average $\mu $ using an outcome regression model $m(x_i)$ to predict the outcome given the covariates, averaging over the predictions $\hat {m}(x_i)$ . When observations are in discrete cells, this is equivalent to taking the modeled regression estimates of the cell averages, $\hat {\mu }_s$ , and post-stratifying them to the population totals as (Gelman and Little Reference Gelman and Little1997):

(10) $$ \begin{align} \hat{\mu}^{\text{omp}} = \frac{1}{N} \sum_i \hat{m}(x_i) = \frac{1}{N}\sum_{s}N_s^{\mathcal{P}} \hat{\mu}_s, \end{align} $$

where $\hat {\mu }_s = 1 / n_s \sum _{i: S_i = s} \hat {m}(x_i)$ is the cell-level average prediction, and where “omp” denotes outcome modeling and post-stratification. For example, we could obtain a prediction from our model of vote choice in each eight-way interacted cell in our running example, and obtain an overall estimate for vote choice by weighting these by the population proportion in each cell. By smoothing estimates across cells, outcome modeling gives estimates of $\hat {\mu }_s$ even for cells with no respondents, thus sidestepping the primary feasibility problem of post-stratification. We discuss particular choices of outcome regression model in Section 4.2.

Heuristically, for weights $\hat {\gamma }$ , we can use an outcome regression model to estimate the bias (conditional on the cell counts) due to imbalance in higher-order interactions by taking the difference between the outcome regression model estimate for the population and a hypothetical estimate with population cell counts $n_s^{\mathcal {R}} \hat {\gamma }(s)$ :

(11) $$ \begin{align} \widehat{\text{bias}} = \hat{\mu}^{\text{omp}} - \frac{1}{N}\sum_sn_s^{\mathcal{R}}\hat{\gamma}(s) \hat{\mu}_s = \frac{1}{N} \sum_s \hat{\mu}_s \times \left(N_s^{\mathcal{P}} - n_s^{\mathcal{R}} \hat{\gamma}(s)\right). \end{align} $$

This uses the outcome model to collapse the imbalance in the J cells into a single diagnostic. Our main proposal is to use this diagnostic to correct for any remaining bias from the multilevel calibration weights. We refer to the estimator as DRP, as it incorporates two forms of “regression”—a regression of the outcome $\hat {\mu }(s)$ and a regression of response $\hat {\gamma }(s)$ through the dual problem, as we discuss in Section B of the Supplementary Material. We construct the estimator using weights $\hat {\gamma }(s)$ and cell estimates $\hat {\mu }_s$ as

(12) $$ \begin{align} \begin{aligned} \hat{\mu}^{\text{drp}}\left(\hat{\gamma}\right) & = \hat{\mu}\left(\hat{\gamma}\right) & + &\;\frac{1}{N}\sum_s \hat{\mu}_s \times \underbrace{\left(N_s^{\mathcal{P}} - n_s^{\mathcal{R}} \hat{\gamma}(s) \right)}_{\text{imbalance in cell } s}\\ & = \hat{\mu}^{\text{omp}} & + & \;\frac{1}{N}\sum_s n_s^{\mathcal{R}} \hat{\gamma}(s) \times \underbrace{(\bar{Y}_s - \hat{\mu}_s)}_{\text{error in cell }s}. \end{aligned} \end{align} $$

The two lines in Equation (12) give two equivalent perspectives on how the DRP estimator adjusts for imbalance. The first line begins with the multilevel calibration estimate $\hat {\mu }(\hat {\gamma })$ and then adjusts for the estimate of the bias using the outcome model $\frac {1}{N}\sum _s \hat {\mu }_s(N_S^{\mathcal {P}} - n_s^{\mathcal {R}}\hat {\gamma }(s))$ . If the population and re-weighted cell counts are substantially different in important cells, the adjustment from the DRP estimator will be large. On the other hand, if the population and re-weighted sample counts are close in all cells then $\hat {\mu }^{\text {drp}}(\hat {\gamma })$ will be close to $\hat {\mu }(\hat {\gamma })$ . In the limiting case of post-stratification where all the counts are equal, the two estimators will be equivalent, $\hat {\mu }^{\text {drp}}(\hat {\gamma }^{\text {ps}}) = \hat {\mu }(\hat {\gamma }^{\text {ps}})$ . The second line instead starts with the outcome regression estimate, $\hat {\mu }^{\text {omp}}$ , and adjusts the estimate based on the error within each cell. If the outcome model has poor fit in cells that have large weight, then the adjustment will be large. This estimator is a special case of augmented approximate balancing weights estimators (e.g., Hirshberg and Wager Reference Hirshberg and Wager2021) and is closely related to generalized regression estimators (Cassel, Sarndal, and Wretman Reference Cassel, Sarndal and Wretman1976), augmented IPW estimators (Chen et al. Reference Chen, Li and Wu2020), and bias-corrected matching estimators (Rubin Reference Rubin1976).

This DRP approach uses outcome information to reduce bias by adjusting for imbalance remaining after weighting. To see this, we can again inspect the estimation error. Analogous to Equation (5), the difference between the DRP estimator and the true population average is

(13) $$ \begin{align} \hat{\mu}^{\text{drp}}\left(\hat{\gamma}\right) - \mu = \frac{1}{N}\sum_{s}\underbrace{\left(n_s^{\mathcal{R}}\hat{\gamma}(s) - N_s^{\mathcal{P}}\right)}_{\text{imbalance in cell }s } \times \underbrace{(\hat{\mu}_s - \mu_s)}_{\text{error in cell } s} + \underbrace{\frac{1}{N}\sum_{s=1}^sn_s^{\mathcal{R}} \hat{\gamma}(s) \bar{\varepsilon}_s}_{\text{noise}}. \end{align} $$

Comparing to Equation (5), where the estimation error depends solely on the imbalance and the true cell averages, we see that the estimation error for DRP depends on the product of the imbalance from the weights and the estimation error from the outcome model. Therefore, if the model is a reasonable predictor for the true cell averages, the estimation error will decrease.

In Section C of the Supplementary Material, we formalize this intuition via finite-population asymptotic theory. We find that as long as the modeled cell averages estimate the true cell averages well enough, the model and the calibration weights combine to ensure that the bias will be small enough to conduct normal-based asymptotic inference. Furthermore, the asymptotic variance will depend on the variance of the residuals $\varepsilon _i$ , which we expect to have much lower variance than the raw outcomes. So asymptotically, the DRP estimator will also have lower variance than the oracle Horvitz–Thompson estimator that uses the true response probabilities, similar to other model-assisted estimators (Breidt and Opsomer Reference Breidt and Opsomer2017). We construct confidence intervals for the population total $\mu $ based on these asymptotic results. First, we start with a plug-in estimate for the variance,

(14) $$ \begin{align} \hat{V} = \frac{1}{N^2} \sum_{i=1}^n R_i \hat{\gamma}(S_i)^2 (Y_i - \hat{\mu}_{S_i})^2. \end{align} $$

We then construct approximate level $\alpha $ confidence intervals via $\hat {\mu }^{\text {drp}}(\hat {\gamma }) \pm z_{1-\alpha /2}\sqrt {\hat {V}}$ , where $z_{1-\alpha /2}$ is the $1-\alpha /2$ quantile of a standard normal distribution.

4.2 Choosing an Outcome Model for Bias Correction

The choice of outcome model is crucial for bias correction. As we discussed in Section 3, we often believe that the strata have important hierarchical structure where main effects and lower-order interactions are more predictive than higher-order interactions. We consider two broad classes of outcome model that accommodate this structure: multilevel regression outcome models, which explicitly regularize higher-order interactions; and tree-based regression models, which implicitly regularize higher-order interactions. Section C of the Supplementary Material gives technical conditions on the outcome regression model under the finite-population asymptotic framework above. We require these regularized models to estimate the true relationship sufficiently well—that is, the regression error of the cell-level estimates $\sum _s(\mu _s - \hat {\mu }_s)^2$ must decrease to zero at a particular rate as the population and sample size increases—depending on how the total number of cells and the population size relate.

4.2.1 Multilevel Outcome Model

We first consider multilevel models, which have a linear form as $\hat {\mu }_s^{\text {mr}} = \hat {\boldsymbol {\eta }}^{\text {mr}} \cdot \boldsymbol {D}_s$ , where $\hat {\boldsymbol {\eta }}^{\text {mr}}$ are the estimated regression coefficients (Gelman and Little Reference Gelman and Little1997; Ghitza and Gelman Reference Ghitza and Gelman2013). MRP directly post-stratifies these model estimates, using the coefficients to predict the value in the population:

$$\begin{align*}\hat{\mu}^{\text{mrp}} = \hat{\boldsymbol{\eta}}^{\text{mr}} \cdot \frac{1}{N}\sum_s N_s^{\mathcal{P}} \boldsymbol{D}_s. \end{align*}$$

This is closely related to the multilevel calibration approach in Section 3.2. If we set $L = -\infty $ and $U = \infty $ in Equation (9)—and so allow for unbounded extrapolation—the resulting estimator will be equivalent to using the maximum a posteriori (MAP) estimate of $\hat {\boldsymbol {\eta }}^{\text {mr}}$ , with regularization hyper-parameters $\lambda ^{(k)}$ (Ben-Michael et al. Reference Ben-Michael, Feller and Rothstein2021). In contrast, the DRP estimator only uses the coefficients to adjust for any remaining imbalance after weighting,

$$\begin{align*}\hat{\mu}^{\text{drp}}(\hat{\gamma}) = \hat{\mu}(\hat{\gamma}) + \hat{\boldsymbol{\eta}}^{\text{mr}} \cdot \left(\frac{1}{N} \sum_s \boldsymbol{D}_s \left(N_s^{\mathcal{P}} - n_s^{\mathcal{R}} \hat{\gamma}(s)\right)\right). \end{align*}$$

This performs bias correction. When we use the MAP estimate of a multilevel outcome model, the corresponding DRP estimator is itself a weighting estimator where the outcome regression model directly adjusts the weights:

$$\begin{align*}\tilde{\gamma}(s) = \hat{\gamma}(s) + \left(\boldsymbol{N}^{\mathcal{P}} - \text{diag}(n^{\mathcal{R}}) \hat{\gamma}\right)'\boldsymbol{D}\left(\boldsymbol{D}'\text{diag}(\boldsymbol{n}^{\mathcal{R}})\boldsymbol{D} + \boldsymbol{Q} \right)^{-1} \boldsymbol{D}_s, \end{align*}$$

where $\boldsymbol {Q}$ is the prior covariance matrix associated with the multilevel model (Breidt and Opsomer Reference Breidt and Opsomer2017). While the multilevel calibration weights $\hat {\gamma }$ are constrained to be nonnegative, DRP weights allow for extrapolation outside the support of the data (Ben-Michael et al. Reference Ben-Michael, Feller and Rothstein2021).

4.2.2 Trees and General Weighting Methods

More generally, we can consider an outcome regression model that smooths out the cell averages, using a weighting function between cells s and $s'$ , $W(s, s')$ , to estimate the population cell average, $\hat {\mu }_s = \sum _{s'}W(s,s')n_{s'}^{\mathcal {R}} \bar {Y}_{s'}$ . A multilevel model is a special case that smooths the cell averages by partially pooling together cells with the same lower-order features. In general, the DRP estimator is again a weighting estimator, with adjusted weights

$$\begin{align*}\tilde{\gamma}(s) = \hat{\gamma}(s) + \sum_{s'}W(s,s')(N_{s'}^{\mathcal{P}} -n_{s'}^{\mathcal{R}}\hat{\gamma}(s')).\\[-24pt] \end{align*}$$

Here the weights are adjusted by a smoothed average of the imbalance in similar cells. In the extreme case, where the weight matrix is diagonal with elements $\frac {1}{n^{\mathcal {R}}}$ , the DRP estimate reduces to the post-stratification estimate, as above.

One important special case is tree-based methods such as those considered by Montgomery and Olivella (Reference Montgomery and Olivella2018) and Bisbee (Reference Bisbee2019). These methods estimate the outcome via bagged regression trees, such as random forests, gradient boosted trees, or Bayesian additive regression trees. These approaches can be viewed as data-adaptive weighting estimators where the weight for cell s and $s'$ , $W(s,s')$ , is proportional to the fraction of trees where cells s and $s'$ share a leaf node (Athey, Tibshirani, and Wager Reference Athey, Tibshirani and Wager2019). Therefore, the DRP estimator will smooth the weights by adjusting them to correct for the imbalance in cells that share many leaf nodes.