
==== Front
8215016
7188
Stat Med
Stat Med
Statistics in medicine
0277-6715
1097-0258

37965978
10.1002/sim.9956
nihpa2021454
Article
Embedded multilevel regression and poststratification: Model-based inference with incomplete auxiliary information
Li Katherine 1
http://orcid.org/0000-0001-8707-7374
Si Yajuan 2
1 Department of Biostatistics, School of Public Health, University of Michigan, Ann Arbor, Michigan, USA
2 Survey Research Center, Institute for Social Research, University of Michigan, Ann Arbor, Michigan, USA
Correspondence: Yajuan Si, Survey Research Center, Institute for Social Research, University of Michigan, ISR 4014, 426 Thompson St, Ann Arbor, MI 48104, USA. yajuan@umich.edu
16 9 2024
30 1 2024
15 11 2023
23 9 2024
43 2 256278
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
Health disparity research often evaluates health outcomes across demographic subgroups. Multilevel regression and poststratification (MRP) is a popular approach for small subgroup estimation as it can stabilize estimates by fitting multilevel models and adjust for selection bias by poststratifying on auxiliary variables, which are population characteristics predictive of the analytic outcome. However, the granularity and quality of the estimates produced by MRP are limited by the availability of the auxiliary variables’ joint distribution; data analysts often only have access to the marginal distributions. To overcome this limitation, we embed the estimation of population cell counts needed for poststratification into the MRP workflow: embedded MRP (EMRP). Under EMRP, we generate synthetic populations of the auxiliary variables before implementing MRP. All sources of estimation uncertainty are propagated with a fully Bayesian framework. Through simulation studies, we compare different methods of generating the synthetic populations and demonstrate EMRP’s improvements over alternatives on the bias-variance tradeoff to yield valid subpopulation inferences of interest. We apply EMRP to the Longitudinal Survey of Wellbeing and estimate food insecurity prevalence among vulnerable groups in New York City. We find that all EMRP estimators can correct for the bias in classical MRP while maintaining lower standard errors and narrower confidence intervals than directly imputing with the weighted finite population Bayesian bootstrap (WFPBB) and design-based estimates. Performances from the EMRP estimators do not differ substantially from each other, though we would generally recommend using the WFPBB-MRP for its consistently high coverage rates.

Bayesian bootstrap
incomplete poststratifiers
sequential imputation
synthetic population
==== Body
pmc1 | INTRODUCTION

Health disparity research often evaluates health outcomes across demographic subgroups, with a focus on vulnerable minors.1 Identification of such at-risk groups provides essential information for health policy research. Multilevel regression and poststratification (MRP) is a method that has recently become popular for subgroup estimation, as it can extrapolate sample inferences to the target population with either probability or nonprobability samples.2–6 MRP has two key components: (1) small subgroup estimation by fitting a predictive multilevel outcome model with a large number of covariates and regularizing with Bayesian prior specifications; and (2) poststratification to adjust for selection bias. The flexible modeling of analytic outcomes can capture complex data structures conditional on poststratification cells—which are determined by the cross-tabulation of categorical variables that affect the sample inclusion (selection and response)—and the use of population control information in poststratification can balance the sample discrepancy.7,8

The availability of population control information strongly related to the analytic variables affects the inferential validity for both model-based and design-based approaches. Poststratification is crucial for sampling selection and nonresponse bias adjustment and requires the use of predictive auxiliary variables with their joint distribution in the population. Practical applications solicit population information from either census records or large studies with minimal errors. For example, across MRP application studies, Si et al3 obtain the joint population control distribution from the American Community Survey (ACS); Wang et al9 use aggregated exit polls; Zhang et al10 use census records; Yougov11 uses the Current Population Survey; Kuriwaki et al12 use the ACS and election official turnout statistics; and Ghitza and Gelman13 turn to large-scale voter registration databases to directly obtain such information for the poststratification adjustment.

However, the population joint distribution of poststratification variables (ie, poststratifiers) is often unavailable, resulting in unknown cell counts; we may only have the marginal distributions of partial variables. We formalize an extension of MRP by embedding the estimation of population cell counts with incomplete auxiliary information, a framework we refer to as “embedded MRP” (EMRP). These estimated cell frequencies are derived from synthetic populations generated from nonparametric bootstrap and sequential imputation approaches, and all sources of estimation uncertainty are propagated under a Bayesian paradigm. EMRP is a class of model-based estimators that accoun for survey weights in modeling and a data integration framework for combining multiple sources of data.

Our primary contributions are that we consolidate different methods under the EMRP class of estimators, implement them to demonstrate the EMRP integrative workflow, and compare their performances with alternatives. We consider three methods for estimating the population cell counts:(1) the weighted finite population Bayesian bootstrap (WFPBB),14 (2) draws from a multinomial distribution informed by observed sampled cell frequencies,2,15 and (3) predictive values from a logistic or multinomial-logit regression that explicitly models the distribution of the missing poststratifier on fully observed variables.16,17 While multinomial draws and logistic regression predictions have established literature use in cell count estimation, this is the first instance (to our knowledge) that the WFPBB has been used for this purpose.

Our methodological research is motivated by the practical operation of an ongoing survey: the New York City (NYC) Longitudinal Survey of Wellbeing (LSW),18 which aims to provide assessments of poverty, material hardship, and child and family general health and wellbeing of the NYC residents. The survey organizers are particularly interested in the life quality aspects of minority groups. The survey collects data from NYC adult residents by oversampling from low-income neighborhoods and following up every three months. We have calibrated the baseline samples to the 2011 ACS-NYC records, assuming all variables affecting the sample inclusion are available in the ACS data: sex, age, race, education, and income.19

We are interested in applying MRP to estimate food insecurity prevalence for sociodemographic subgroups. Since food insecurity is associated with an increased risk of diabetes and hypertension in adults, it is of interest to healthcare professionals and policymakers to identify subpopulations at risk for food insecurity to improve health outcomes.20 Weighted estimation can inflate estimation variability—especially for small groups—and we would likely obtain more stable estimates with MRP. As residents who often visit food acquisition agencies in NYC tend to suffer from food insecurity and material hardship in addition to having a different sample inclusion rate, MRP inferences with the LSW sample should adjust for visiting frequency. However, we are unable to do so as the population distribution of visiting frequency is unknown.

The problems we face in the LSW baseline survey are reflective of problems in most survey applications. When combining multiple data sources, we often encounter incomplete auxiliary information. We would like to integrate the estimation of poststratifier distributions into the survey outcome modeling under EMRP. We consider nonparametric bootstrap algorithms and parametric models for the population cell count estimation.

The paper structure is organized as follows. Section 2 outlines the methods that fall within the EMRP framework. Section 3 provides a simulation study to evaluate the performances of different EMRP methods. We apply EMRP to estimate food insecurity prevalence using the LSW dataset in Section 4. Finally, Section 5 summarizes the findings and extensions.

2 | METHODS

Suppose the outcome in the population is Yi, and the auxiliary variables of the population are categorical (or discretized continuous variables) denoted by Xi and Zi, for i=1,…,N, where N is the population size and either Xi or Zi can be multivariate. For illustrative purposes, we assume that Xi is univariate and that Zi represents multiple variables. The sample of size n drawn from the population includes zi,xi,yi, for i=1,…,n.

When the population distribution of Z is known, the cross-tabulation of Z results in M cells with known cell sizes Nmz′s with ∑m=1MNmz=N. If our goal is to estimate the overall population mean of the outcome θ=∑iYi/N, then the classical MRP estimator can be described as follows: (1) θ^MRPz=∑m=1MNmzNθ^mz,

where θˆmz is the model-based estimate in cell m for m=1,…,M. The outcome model (Y∣Z) fitted to the sample data yi,zi can be a Bayesian multilevel regression. Examples include a Bayesian hierarchical model with weakly informative or informative prior specifications, or a flexible prediction algorithm.3,9,21,22 With a binary outcome yi∈{0,1}, the model can be a Bayesian logistic regression, Gaussian process regression model,2 or stacked regression23 by assuming that yi~Bernoulliθj[i], where θj[i]=Pryi=1=fαiz is a function of parameters αiz corresponding to the level of z for unit i and cell j[i] that unit i belongs to.

However, the poststratifiers’ information can be incomplete. The cross-tabulation of all the poststratifiers Zi,Xi results in J cells. If the joint population distribution of Zi,Xi is unknown, then the population cell sizes must be estimated by Nˆj ’s, with ∑j=1JNˆj=N. We can embed the estimation of the poststratification weights into the classical MRP workflow and present the conceptual EMRP workflow in Figure 1.

(2) θ^EMRP=∑jJN^jNθ^j.

The differences between MRP (1) and EMRP (2) estimators are as follows.

The poststratification variables of EMRP are (Z,X), while MRP poststratifies only on Z. Poststratifiers should be predictive of either the analytical outcome (primarily) or the response propensity (secondarily).4,24 We assume that the sample inclusion mechanism depends only on Z such that the data are missing at random (MAR). If sample inclusion depends on incomplete auxiliary variables, then the data are missing not at random (MNAR) and both MRP and EMRP estimators are subject to bias. Nevertheless, EMRP is expected to reduce the bias of MRP by leveraging the correlation structure between X and Z to generate the joint distribution, which will be elaborated on below.

EMRP estimates the population cell size Nˆj and propagates its uncertainty, while MRP treats Nmz as fixed. The propagated variance estimate of EMRP is larger than that of the MRP estimator with known Nj, which is expressed in the following decomposition: var(θ^EMRP)=var[E(∑jN^jNθ^j|N^j)]+E[var(∑jN^jNθ^j|N^j)]=var[∑jN^jNE(θ^j∣N^j)]+E[var(∑jN^jNθ^j|N^j)],

where the second term is the variance of the MRP estimator poststratifying on (Z,X):varθˆMRPzx=var∑jNjNθˆj with known Nj.

The cell estimate θˆj under EMRP is based on the model (Y∣Z,X) fitted to the sample data (y,z,x), while classical MRP fits the model (Y∣Z) to the sample data (y,z). Considering the subpopulation mean E(Y∣Z)=∑XE(Y∣Z,X)Pr(X∣Z), the MRP estimator can be similar to the EMRP estimator asymptotically or when the sample resembles the population distribution of X given Z in the subgroup to achieve the equality. However, because of small subgroup sizes in the sample data, it is possible that the observed frequency of X values within a subgroup defined by Z is different from the population distribution and the equality will not hold. We expect that the estimation of E(Y∣Z,X) will be more accurate than E(Y∣Z) when X is a strong predictor of Y.

Hence, we expect that the EMRP estimator can reduce the bias of the MRP estimator, especially for small subgroups. We compare three strategies to estimate Nˆj based on Nmz and the sample data (y,z,x) : (1) drawing synthetic populations of (Z,X) from the weighted finite population Bayesian bootstrap (WFPBB-MRP), (2) drawing Nˆj from a multinomial distribution (Multinomial MRP), and (3) predicting X values in the population by regressing X on Z with logistic or multinomial-logit models (Two-stage MRP). In all of our methods, the estimation uncertainty of Nˆj is propagated under a Bayesian framework.

2.1 | The weighted finite population Bayesian bootstrap (WFPBB-MRP)

Dong et al14 propose the WFPBB as a nonparametric method of generating synthetic populations that can be analyzed as simple random samples by “undoing” the complex sampling design and accounting for the sampling weights. We use the WFPBB to estimate the joint distribution (Z,X) in the population. Let nmz denote the number of sampled units with z=m, and we have ∑m=1Mnmz=n. We construct the sampling weights for unit i,wi=Nm[i]z/nm[i]z, where m[i] represents the value of z that is assigned to unit i, for i=1,…,n. The idea is to draw from the posterior predictive distribution of non-observed (nob) data given the observed (obs) data and weights: zi,xinob∣zi,xiobs,wi. Assume that there is a finite number of unique pairs zi,xi in the sample, the population is also comprised of these unique pairs, and the corresponding counts for each pair in the population follow a multinomial distribution. Given a non-informative Dirichlet prior distribution on the multinomial probabilities, the Pólya distribution can be used in place of the Dirichlet-multinomial distribution to draw predictive samples and reduce the computational burden. We adapt and embed the WFPBB in the MRP implementation as “WFPBB-MRP” shown in Figure 2, following the steps below.

Resample via Bayesian bootstrap (BB):25 To capture the sampling variability of the “parent” (original) sample, we generate L number of BB samples: B1,…,BL, each of size n.

Recalibrate weights: For each BB sample Bl, we recalibrate the bootstrap weights by multiplying the base weights by the number of replicates for unit i in Blril and normalizing the weights to sum to the population size N, so that wil=Nwiril∑jwirjl.

Use the WFPBB to incorporate weights: Construct the initial Pólya urn based on zi,xi with their corresponding replicate weights wil’s and draw N-n units with probability (3) wil-1+li,k-1(N-n)/nN-n+(k-1)(N-n)/n,

for the kth draw, k∈{1,…,(N-n)}, where li,k-1 is the number of bootstrap selections of zi,xi among the elements present in our urn at the k-1 draw. The draws form the WFPBB sample Sl,f of size N. We repeat this step F times, yielding multiple samples Sl,1,…,Sl,F-which we pool to create the synthetic population Sl— of size F*N, where Sl is considered a single draw from the WFPBB.

Make inference: For each of the samples Sl,f∈Sl,1,…,Sl,F we obtain estimates Nˆj(l,f). The estimate of Nˆj associated with synthetic population Sl is the average of the estimates from the F samples. We normalize the sum of Nˆj to be N and obtain L samples of Nˆj. As a parallel step, we obtain posterior samples of θˆj from the multilevel regression. The embedding of the two estimates yields the posterior samples of θˆEMRP; see Figure 1.

In the WFPBB implementation, we have constructed weights based on the population distribution of Z. In practice, if the sampling weights are available from complex survey data, the weighted counts of the joint cells of (Z,X) can be used to estimate their population distributions, but these point estimates have ignored the sampling uncertainty of the survey. Similar to replication approaches that are often recommended to propagate the sampling variance, WFPBB undoes the weights in the generation of multiple synthetic populations and obtains the posterior samples of the population counts.

An alternative method for the population mean estimation is directly imputing the outcome of Y with the WFPBB in a manner similar to that of Dong et al.14 The motivation of our adaptation is to improve estimation accuracy and precision using predictive auxiliary information. Direct WFPBB-imputation inference tends to become unstable with large standard errors when the sampled cell count is small or the replicate weights are highly variable; in such cases, the synthetic populations are mainly informed by a few outcome values observed in the sample and require a large number of replicates to be able to represent the true population. By using predictive auxiliary information in a multilevel model, we borrow inference across cells to improve estimation accuracy and precision. We compare the performance of the WFPBB-MRP to that of the direct WFPBB implementation in the simulation studies.

2.2 | Multinomial draws (multinomial MRP)

Using the WFPBB can be computationally demanding. An alternative is to approximate the weighted Bayesian bootstrap by drawing the cell counts from a probability distribution—for example, a Poisson or multinomial distribution—as in Si et al2 and Makela et al.15 Suppose the possible values of the categorical variable X (which can be multivariate) are (1,…,C), where C is the total number of levels for X. Here, we consider a multinomial distribution for X conditional on the fully observed Z.

(4) X=1,…,C∣Z=m~MultinomialNmz;p1,…,pC,pc=nm,czx/nmz,c=1,…,C.

We use the marginal population counts Nmz and the observed frequency in the sample nm,cz,x/hmz, where nm,cz,x denotes the sample count of units with zi=m,xi=c, for m=1,…,M and c=1,…,C. This distribution accounts for the correlation among variables based on the sample data. We extend MRP with synthetic poststratification (Multinomial MRP) to address the estimation of the joint distribution of (Z,X).

Leemann and Wasserfallen26 assume both population margins of (Z,X) are available (whereas any population information for X is missing in our setting) and fit a multinomial model similar to (4), except that the observed frequencies nm,cz,x/nmz are updated to match the known margins of X. The approach taken by Leemann and Wasserfallen26 is similar to raking when population margins of multiple variables are available,27–29 but draws repeated bootstrap samples to account for uncertainty in the synthetic poststratification. Multinomial MRP is a fully-Bayesian procedure that uses posterior samples to propagate all sources of estimation uncertainty.

2.3 Sequential regression models (two-stage MRP)

Reilly et al16 apply regression models to predict the unknown population poststratifier. Since MRP essentially predicts outcomes in the population, Kastellec et al17 propose a sequential estimation procedure by using two separate MRP procedures, that is, Two-stage MRP, where the first stage estimates the population cell counts Nˆj using a multilevel model for (X∣Z), and the second estimates the cell means θˆj with (Y∣Z,X). The model (X∣Z) assumes the probability pc in Model (4) is a function of pre-specified main effects or high-order interaction terms of Z variables, but not necessarily the cell-wise observed frequencies across the cross-tabulation, resulting in regularized estimates with a strong dependency on model specification. If X is a binary indicator (with values 0/1), the first stage regresses X on Z. (5) logitPrXi=1∣Zi=βiz,

where the coefficient βiz denotes the main or high-order effect corresponding to the level of z for unit i. The synthetic predictions (Xi∣Zi) would yield an estimate of the joint frequency in cell j[i], that is, Nˆj, for j=1,…,J. This method becomes cumbersome when X has more than two levels, demanding a multinomial regression model.

The three EMRP strategies presented rely on the availability of population frequencies of Z variables to generate synthetic populations of (Z,X) and assume that the relationships between X and Z are the same between the sample and the population. The WFPBB-MRP constructs its base weights from the cross-tabulation of Z variables and draws nonparametric Bayesian bootstrap samples. The Multinomial MRP assumes that the conditional distribution (X∣Z) within each cell in the cross-tabulation of Z variables is a multinomial distribution with probabilities set as the observed frequencies. The Two-stage MRP applies a Bayesian multilevel model to predict X given Z, where the model in practice often includes only the main effects of Z.

The Two-stage MRP is subject to model misspecification and has stronger modeling assumptions than the WFPBB-MRP and Multinomial MRP, which automatically consider high-order interaction terms between X and Z. Both the WFPBB-MRP and Multinomial MRP use the observed conditional distributions and may exhibit robust performances with large sample sizes. However, if the sample cell counts in the contingency table of Z are sparse the estimated Nˆj from these two methods may have large variances, though the parametric Multinomial MRP has intrinsic variances that may be smaller than those from the WFPBB-MRP. It is possible that some Z values are not present in the sample, resulting in empty cells. While WFPBB-MRP and Multinomial MRP omit empty cells, Two-stage MRP can generate the predictions of empty cells in the population since it relies on the model structure rather than the observed frequencies. All methods propagate the uncertainty in estimating the unknown population counts Nˆj in the EMRP variance estimator varθˆEMRP, but there is a slight difference: WFPBB-MRP accounts for both Z and the modeling uncertainty of estimating (X∣Z). Multinomial-MRP and Two-stage MRP only account for the uncertainty of estimating (X∣Z) while treating the external distribution of Z as fixed. Previous work16,17,26 also ignores the sampling error of Z and only accounts for the modeling uncertainty of (X∣Z). Since the external data Z are often of large sample size—as is the case in the ACS—the counts are approximated as the population data with negligible sampling error. However, the sampling error would become substantial if the external sample size is small.

3 | SIMULATION STUDIES

We conduct simulation studies to compare embedded MRP methods, direct WFPBB imputation of the outcome, and classical MRP using (1) for the overall population and subdomain inferences, both between the classes of estimators (WFPBB vs. classical MRP vs. EMRP overall) and within the class of EMRP estimators (WFPBB-MRP vs. Multinomial MRP vs. Two-stage MRP). The simulation code is publicly available here.

3.1 | Setup

We simulate a population of size N=10,000 with a binary outcome Y and four categorical variables for poststratification: three Z variables with known population distributions and one binary X variable (0/1) without population information. The first two Z variables (Za and Zb) have five levels and are generated from multinomial distributions, and the third variable Zc is binary and drawn from a binomial distribution, where probabilities are normalized random numbers drawn from the uniform distribution, Uniform (0.25, 1). The cross-tabulation of Z variables results in M=50(=5×5×2) cells. The variable X is generated conditional on Z with probability (6) PrXi=1∣Zia,Zib,Zic=expitβ0+βiZa+βiZb+βiZc+βiZa,Zc+βiZb,Zc,

where the function expit(μ)=exp(μ)/(1+exp(μ)),Zia,Zib,Zic denote the values of Za,Zb,Zc for unit i∈{1,…,N}, respectively, and βivar represents the value of βvar corresponding to the level of var for unit i, for var∈Za,Zb,Zc. We assume that Model (6) has both main effects and interaction terms (INT) with β0=-0.5,βZa=(1.7,0.25,0.2,-0.75,-1.7)⊤,βZb=(2.3,1.5,0.15,0.2,0.9)⊤,βZc=(0,-1)⊤,βZaZc=(0,-0.6,0.5,0.35,-0.4)⊤, and βZb,Zc=(0,1.7,0.1,2,-0.75)⊤. We also consider the case when Model (6) has only main effects of Z and include the output in Appendix B.

We specify the data generating process (DGP) for Y as: (7) Yi∣Zia,Zib,Zic,Xi~Bernoulliexpitα0+αiZa+αiZb+αiZc+αiX.

The values are assigned as α0=0,αZa=(1.37,-0.56,0.36,0.63,0.40)⊤,αZb=(-0.11,1.51,-0.09,2.02,-0.06)⊤,αZc=(0,0.24)⊤, and αx=(0,-1.3)⊤, where αZa,αZb are drawn from a standard normal distribution.

We assume that the inclusion mechanism depends only on the fully observed Z. The inclusion probabilities Pr(I=1∣Z) are based on the cross-tabulation of three Z variables, Zcat∈{1,…,M}, with values drawn from different ranges. Table 1 gives the cell indices and ranges of inclusion probabilities. We randomly draw values with replacement from each interval of equally spaced probabilities and assign them to the corresponding cells. We draw 200 repeated samples from the population with the pre-specified inclusion mechanism Pr(I=1∣Z). The population cell frequencies range from 3 to 470, with an average of 100. In the repeated studies, the resulting average overall sample size is 4104; we obtained similar results with a smaller average sample size of 2052.

We are interested in the overall population and subgroup mean estimates. We assume that the MAR sample inclusion mechanism depends only on Z. Within each poststratification cell, MRP assumes that the inclusion probability is the same, and the individuals are independently and identically distributed. EMRP uses the correlation between X and Z to infer the cell counts for poststratification with (Z,X). The inclusion probabilities can therefore be treated as a summary statistic of the distribution of Z and their correlation with the inclusion mechanism. We use the intervals of inclusion probabilities to create subgroups of interest.

We consider two methods to create subgroups. First, we create four subgroups based on the percentiles of the inclusion probabilities of the J cells and the distribution of (X∣Z). Each subgroup contains 20 cells: the first subgroup includes cells with inclusion probabilities lower than or equal to the 40th percentile; the second contains those with moderate inclusion probabilities between the 20th and 60th percentiles; the third group includes those with medium inclusion probabilities between the 40th and 80th percentiles, and the fourth group covers high inclusion probabilities above or equal to the 60th percentile. The X categories have different frequencies across subgroups, where for the first, second, and fourth subgroup, 15 of the 20 cells have X=0 and 5 have X=1; for the third subgroup, 5 cells have X=0, and 15 have X=1. The average sampled cell sizes in the (Z,X) cross-tabulation table for subgroups with low pI:0.03-0.27, moderate pI:0.21-0.38, medium pI:0.30-0.54, and high pI:0.40-0.95 are (9, 24, 53, 87) with average subgroup sample sizes (198, 491, 1064, 1744). The simulation scenarios cover cases with sparse cells and small groups.

Second, the four subgroups are based on cells in the cross-tabulation of only Z variables that are fully observed poststratifiers. All parameters are identical to those in the first subgrouping definition except for the subsampling procedure in each inclusion probability bracket: subgroups are defined by a random sample of 10 cells based on the cross-tabulation of Z variables in each bracket instead of creating an imbalance of X categories. The average sampled cell sizes in the Z cross-tabulation table for subgroups with low pI:0.03-0.25, moderate pI:0.25-0.38, medium pI:0.34-0.51, and high pI:0.40-0.93 are (13, 28, 40, 62) with average subgroup sample sizes (262, 567, 815, 1251).

Table 2 presents the differences between the observed values in one random sample and the population values of the missing poststratifying variable’s frequency distributions Pr(X=1) within the subgroups of two cases: (1) subgroup membership is defined based on the joint (Z,X) distribution and (2) membership is defined based on categories of Z only, which shows that the first scenario generally has larger differences than the second scenario. We expect that larger differences in terms of Pr(X=1) lead to more different EMRP and MRP estimates.

For all three EMRP methods, the outcome model fitted to the sample data is identical to the DGP of Y (see (7)).

The outcome model for classical MRP omits the main effect for X. The estimation model for (X∣Z) in the Two-stage MRP only accounts for the main effects, as misspecified: logitPrXi=1∣Zi=β0+βiZa+βiZb+βiZc.

To implement the WFPBB, we use the polyapost package version 1.630 to draw from the weighted Pólya posterior distributions and generate L=1000 synthetic populations of size F*N=20*10,000 per sample for inferences. We use Stan for the fully Bayesian posterior computation with MRP and perform convergence diagnostics.31 For each sample, we fit the outcome model using two Markov chain Monte Carlo (MCMC) chains with 2000 iterations and keep the last 500 iterations from each for a total of 1000 draws (permuted and merged across chains) for estimation. Regression coefficients of multiple categories are assigned weakly informative priors:32 normal distributions with mean 0 and unknown standard deviation parameters that are assigned hyperpriors of Cauchy+(0,1), where Cauchy+(0,1) is the half-Cauchy distribution restricted to positive values with a standard deviation 1. For the intercept and coefficients of binary predictors, we use noninformative priors.

We assess bias, root mean squared error (rMSE), the average length of 95% confidence intervals (CI length), and the nominal coverage rate of 95% CIs (coverage rates) for the finite population quantities of interest. To quantify the uncertainty surrounding the Nˆj estimation in EMRP, we pair each of the 1000 sets of cell mean estimates θˆj with one of 1000 sets of Nˆj. This means we construct the sets of Nˆj from 1000 synthetic populations of size F*N for WFPBB-MRP, 1000 draws of Nˆj from the multinomial distribution for Multinomial MRP, and 1000 posterior draws from the sampler fitting the logistic regression model for (X∣Z) for Two-stage MRP. We compute 95% CIs of the estimates of all five methods for a given sample by taking the 2.5th and 97.5th percentiles of their respective 1000 posterior estimates.

We also present the survey-weighted and unweighted mean estimators. We use the same base weights as those in WFPBB (wi=Nm[i]z/nm[i]z) and the survey package version 4.1–133 to obtain the design-based estimates, standard errors with the finite population correction, and 95% confidence intervals.

3.2 | Results

We present the simulation results by plotting heatmaps in Figure 3 and reporting detailed values in Tables 3 and 4, for the two subgrouping methods, respectively.

When the subgroup is defined by the joint distribution of (Z,X), as shown in the left column of Figure 3, the classical MRP yields high bias values and near-zero coverage rates due to the differences between the observed distribution of (X∣Z) and the population distribution in the subgroups. The EMRP estimators correct for these deficiencies: for subgroup estimates, classical MRP has absolute bias values of at least 0.033 and coverage rates of at most 0.005 while the EMRP methods produce absolute bias values of at most 0.011 and coverage rates of 0.870 or above. The EMRP methods do not differ substantially; the rMSE values differ by at most 0.003 and the bias values by 0.006. While WFPBB-MRP has a coverage rate of at least 0.99 in all subdomains, Multinomial MRP and Two-stage MRP have slight undercoverage in subgroups with moderate and medium inclusion probabilities pI:0.21-0.38 and pI:0.30-0.54:(0.885,0.895) and (0.910, 0.870), respectively.

Table 3 shows that the survey-weighted estimator has competitive bias and coverage rates near or above 95% except in the subgroup with low inclusion probabilities. Its rMSE and bias values are comparable to those from imputing the outcome with WFPBB, but its CIs are much narrower. The stabilizing benefit of Bayesian multilevel modeling is apparent when we compare direct imputation with WFPBB with EMRP methods. While the former yields comparable bias, the large variation between synthetic population inflates its rMSE and CI width compared to EMRP methods. This results in conservative coverage rates in most subgroups, though the bias in the low inclusion probability subgroup pI:0.03-0.27 incurred by the variable weights and sparse sampling results in slight undercoverage (0.940) despite wide intervals: the CI of the WFPBB estimate for the low group spans 0.202 compared to 0.105 from the WFPBB-MRP and 0.073 from the Multinomial-MRP and Two-stage MRP.

The right column of Figure 3 shows that the MRP and EMRP estimates are similar for subgroups defined by the categories of Z only, where the adjustment of incomplete X in EMRP does not provide additional gains. Consistent with Table 2, the largest difference in the Pr(X=1) is in the group with low inclusion probabilities (pI:0.03-0.25), for which the EMRP and MRP estimates have the most prominent dissimilarity among the five estimates. Table 4 shows that WFPBB and survey-weighted estimators tend to have larger variances than the EMRP and MRP estimates, especially for the subgroup with low inclusion probabilities.

Figure 4 presents the performance metrics for the EMRP population cell frequency Nj estimation in comparison of different methods. Estimating Nj with the WFPBB results in the widest CIs among the EMRP methods for Nˆj and, subsequently, the subgroup estimates. For Two-stage MRP, misspecification of the (X∣Z) logistic model results in Nj estimates that are substantially more biased than the non-regression methods and contaminates the performance of Two-stage MRP, as shown in Figure 3. Multinomial MRP has lower bias values and interval lengths comparable to the Two-stage MRP. However, both methods underestimate the uncertainty in estimating Nˆj and thus lead to low coverage rates of the corresponding EMRP estimators.

Overall, EMRP estimators have higher precision than the design-based methods and smaller bias values than the classical MRP estimator when the distribution (X∣Z) in the observed sample is different from the population. WFPBB-MRP has larger variances and conservative CI coverage compared to Multinomial MRP and Two-stage MRP.

4 | APPLICATION TO THE LONGITUDINAL SURVEY OF WELLBEING

The LSW dataset is comprised of two different samples: a phone sample of 2002 residents contacted by random digit dialing, and a face-to-face sample of 226 residents visiting food acquisition agencies for a total of 2228 respondents. We analyze the publicly released data, which do not distinguish the two samples even though they have different sample inclusion mechanisms. Integrating the phone and face-to-face samples will be discussed as a future extension in Section 5. The study oversamples residents from low-income neighborhoods. The publicly released survey weights have been calibrated to the ACS-NYC 2011 weighted totals and account for unequal probabilities of selection, undercoverage, and nonresponse. While these weights are available in this particular composite sample, the combination is ad hoc and requires improvement with rigorous data integration methods. The frequency of residents visiting a food acquisition agency is related to their food insecurity status; however, we do not have access to the population distribution of the agency visit frequency and will need to estimate it for our analysis.

We classify a respondent to be food insecure if it is “often” the case that they “worried whether [their] food would run out before [they] got money to buy more.” We use the binary proxy of agency visit frequency to indicate whether they have visited a food acquisition agency within the last 12 months. The outcome model in the EMRP methods accounts for age in years (18–35, 36–50, 50+), sex (male, female), race (White, Black, other), the highest level of education achieved (less than high school, high school or equivalent, some college or associate’s degree, bachelor’s degree or higher), annual pre-tax cash income for the household (<$35k, $35–55k, $55–100k, >$100k), and agency visitor status (“visitor” if the respondent has visited a food acquisition agency within the last 12 months, “non-visitor” otherwise). There are 11 participants with a missing response variable and 20 missing their education values which we impute by randomly sampling from the corresponding observed values, and the effect of the small amount of item nonresponse is negligible.

Table 5 gives the distributions of sociodemographics and the food insecurity prevalence of the LSW sample and two groups stratified by the indicator of agency visits. Agency visitors tend to be younger, non-White, less educated, and more food insecure compared to non-visitors. Interestingly, the frequency of agency visitors who are low-income residents (<$35k) is lower than that of nonvisitors (32.6% vs. 57.5%). The population distribution of sociodemographics (age, sex, race, education, and income) is available in the ACS-NYC 2011, but that of agency visitors is not, reflecting the EMRP setting in Figure 1.

We apply the EMRP methods to the LSW study to estimate the prevalence of food insecurity among NYC adult residents with different income levels. Two sets of analyses are conducted: one set focuses on the income groups (<$35k,n=884;$35-55k,n=273;$55-100k,n=456; and >$100k,n=615) to compare EMRP and classical MRP, and the other set is defined by cross-tabulations of income levels and agency visitor status. In both sets of subgroups, we also include results from directly imputing the outcome with the WFPBB and from the survey-weighted estimator (using publicly released weights).

For the WFPBB-MRP, we use observed sociodemographic frequencies from the sample nmz and their weighted totals from the ACS Nmz to obtain the initial weights in the Pólya urn wm=Nmz/nmz and apply WFPBB to estimate Nˆj (where cell j is from the contingency table based on (Z,X)) with L=5000 synthetic populations and F=20 draws from the weighted Pólya posterior. To reduce the computational burden, we modify the size of each draw f(=1,…,F) as T*n=30*2228-which is large enough to overwhelm the sample size-instead of synthesizing the entire population. We set L as a large value of 5000 to match the number of 5000 posterior samples from EMRP. The same parameter settings and base weights are used when directly imputing with the WFPBB.

The outcome model in EMRP includes sociodemographic variables and visit status: logitPryi=1∣xi,zi=α0+αiage+αisex+αirace+αieduc+αiincome+αivisit+αivisit:income,

where yi=1 indicates food insecurity and yi=0 otherwise, for i=1,…,n. The model includes all main effects and the two-way interaction between visit frequency and income. Classical MRP omits the covariate of visit frequency: logitPryi=1∣zi=α0+αiage+αisex+αirace+αieduc+αiincome, and the same set of covariates are used in the estimation model for visit frequency (X∣Z) in Two-stage MRP.

We assign the same weakly informative prior distributions to the coefficients as those in Section 3.1. The model runs four MCMC chains with 10,000 iterations each, with 8750 for burn-in and the last 1250 iterations permuted and merged across chains for a total of 5000 iterations for estimation. The procedures for conducting subgroup inferences from the EMRP and WFPBB estimators are the same as those described in the simulations. The diagnostics indicate model convergence.

We use Bayesian leave-one-out cross-validation and posterior predictive check to evaluate the goodness of fit,34,35 neither of which raises concerns about model performances. We compare the pointwise out-of-sample prediction accuracy between the EMRP and MRP models. Pareto k estimates are less than 0.7 for both models, implying that all leave-one-out posteriors are similar to the full posterior. The expected log predictive density for the EMRP model is greater than that of the MRP model (18.4 difference); the EMRP model is preferred for prediction. Details are given in Appendix C.

The estimated coefficient for the agency visitor indicator is αˆvisit=0.83 (95% CI: 0.48, 1.18), and estimates for its interaction terms with income at the <$35k, $35–55k, $55–100k, >$100k levels are −0.24 (95% CI: −1.34, 0.88), 0.31(95% CI: −0.78, 1.67), 0.54 (95% CI: −0.44, 1.95), and −0.60 (95% CI: −2.29, 0.45), respectively. The main effect of X is larger than 0, but the interaction effects have large variability.

Figures 5 and 6 compare the food insecurity prevalence estimates from the unweighted and survey-weighted analysis, classical MRP, direct imputation of the outcome using the WFPBB, and EMRP for the overall NYC adult population and the four subdomains defined by: (1) annual income intervals and (2) the cross-tabulation between income and agency visit status, respectively.

In Figure 5, the survey-weighted, classical MRP, and EMRP overall prevalence estimates are around 9%, which is slightly lower than the unweighted estimate of 10.1% and the WFPBB direct imputation estimate of 11.4% (see Table A1 for detailed values). The unweighted estimates are generally higher than the weighted estimates with the exception of the $55–100k group. The model-based estimates show an inverse relationship between food insecurity and income: around 16% of those with <$35k annual income are food insecure compared to 2% of those in the >$100k group. Incorporating weights in the analysis is crucial for bias adjustment. The variances of EMRP estimators are smaller than those of the weighted and WFPBB estimators, where predictive models improve precision. Different EMRP estimators yield similar estimates and overlapping intervals. Classical MRP point estimates are similar to the EMRP point estimates for the income categories. This is possibly due to small differences in the agent visit frequency between the sample and the population within each income group, the sample sizes of which are large enough.

Figure 6 presents the prevalence estimates for the subgroups defined by interactions of annual income and agency visitor status. Agency visitors have higher food insecurity compared to those who don’t visit, and this gap lessens as annual income increases. The largest gap occurs for individuals with annual income lower than $35k: approximately 11% versus 25%. Here, the conditional estimates given (Z,X) are substantially different from those only given Z. EMRP is beneficial if we are interested in conducting inference on subgroups defined by the missing poststratifier X. Agency visitors tend to be more food insecure than the overall study population; the weighted estimate is 20.1%, more than double the 9.3% estimated for the overall population (Table A2).

The application study shows that accounting for design features is important for correcting the bias in the food insecurity prevalence estimation. The agency visit status is an important predictor of the food insecurity outcome, the distributional imbalance of which between the sample and population within subgroups will affect the mean estimates.

5 | DISCUSSION

Motivated by health disparity research, we focus on estimates for minority groups. MRP has become a popular subgroup estimation method due to its ability to stabilize estimates and adjust for selection bias, but these properties are restricted by whether the population joint distribution of poststratifying auxiliary variables is known. Analysts rarely have access to the population joint distribution of the comprehensive set of predictive auxiliary variables. In such situations, classical use of MRP may require the omission of predictive variables without complete information and produce inaccurate estimates as a result. We have developed the EMRP framework to incorporate variables with incomplete information into MRP by generating synthetic populations to estimate their joint distribution before proceeding with MRP estimation.

Through simulation studies, we compared design-based estimators and direct imputation with the WFPBB with EMRP methods We found that all EMRP estimators can correct for the bias in classical MRP while maintaining lower standard errors and narrower confidence intervals than directly imputing with the WFPBB or using design-based estimators. Performances from the EMRP estimators do not differ substantially from each other, though we would generally recommend the WFPBB-MRP for its consistently high coverage rates. As a benefit of fitting multilevel models and stabilizing small group estimates, the WFPBB-MRP yields bias and rMSE values that are comparable to the Multinomial MRP and Two-stage MRP while producing narrower confidence intervals than directly imputing with the WFPBB. Estimating population cell frequencies using the multinomial distribution leads to small bias values and rMSE for the Nˆj estimates, but low coverage rates of the frequencies in sparsely sampled cells lead to undercoverage for the Multinomial MRP in a few subgroup inferences. Conversely, using MRP to recover Nj as in the Two-stage MRP gives reasonable coverage in most cells, but generates the largest bias and rMSE values among all methods. Model misspecification for (X∣Z) can introduce bias in domain inferences, especially for domains with few observations. The WFPBB-MRP avoids this issue by weighting observed cases, accounts for sampling uncertainty when estimating the joint (Z,X) distribution, and combines with MRP to improve the inferences for (Y∣Z,X) with a predictive model.

Under settings with incomplete poststratifier X information, differences in bias between classical MRP and EMRP methods are contingent on an imbalance in the missing poststratifying variable’s frequency distributions for the inferential subgroup between the sample and population (eg, a special case is that not all legitimate X values for a given Z in the population are observed in the sample). When the inclusion mechanism is MAR given Z, we would expect that the overall mean estimates of MRP and EMRP are similar, but the subgroup estimates could be substantially different. EMRP improves the estimates for subgroups defined by (Z,X). We have compared the differences of the missing post-stratifying variable’s frequency distributions Pr(X=1) between the population and one randomly drawn sample across the four subgroups constructed under two scenarios in Section 3.1: (1) based on the joint distribution (Z,X), where four subgroups are based on the percentiles of the inclusion probabilities of the J cells and the distribution of (X∣Z); and (2) based on only Z, where four subgroups by a random sample of 10 cells based on the cross-tabulation of Z variables in each inclusion probability bracket. Table 2 shows that the first scenario generally has larger differences than the second scenario. The largest difference in the Pr(X=1) is in the group with the lowest inclusion probabilities, where the EMRP and MRP estimates have the most prominent dissimilarity. If the analyst anticipates that their inferential subgroup will have a balanced (X∣Z) distribution between the sample and population, then the MRP estimator will have similar bias values to those of the EMRP estimators. The utility of EMRP is best showcased when there is an imbalance of (X∣Z) in the subgroup between the sample and the population.

The EMRP framework has a few interesting directions for future extensions. First, the EMRP framework can handle general problems of data integration. Datasets from different sources may have incongruous study measures, resulting in incomplete auxiliary information. When combining two datasets obtained through different sampling mechanisms, we modify the Nˆj estimation procedure such that we use only one of the datasets—the one with a selection mechanism that is independent of the incomplete auxiliary variables—to estimate the (X∣Z) distribution. To illustrate, the LSW dataset in our application is a composite sample of phone and face-to-face surveys. Given the indicator of the two survey components in the restricted dataset, we would regress X on Z using only the phone sample for the Two-stage MRP, use sociodemographic cell frequencies from the phone sample as the initial sampling weights for the WFPBB-MRP, and use the phone sample to estimate the probabilities for the Multinomial MRP. Integrating data from multiple studies can also present multivariate incomplete auxiliary variables. Such a task would require an iterative or sequential estimation process to estimate the joint distribution, where the choice of WFPBB, multinomial, or MRP models can be tailored to each auxiliary variable and combined under a framework that is similar to multiple imputation.

Second, in situations where marginal distributions of the multivariate incomplete auxiliary variables are available but their joint distribution is unknown—similar to the raking setting—the WFPBB needs to account for the known constraints when constructing the base weights. Model-based estimation approaches under known margins can be applied.29

Third, EMRP can be extended to settings where the data are MNAR because the inclusion mechanism depends on the incomplete auxiliary variables. For example, the inclusion of face-to-face samples in the LSW study depends on the agency visit frequency. When evaluating the performance of EMRP methods under MNAR in our simulation studies, we found that all methods yield bias, but EMRP reduces the bias of classical MRP. Enhancing EMRP methods to handle informative inclusion would further broaden the circumstances under which MRP can be applied successfully.

Finally, the wide use of the EMRP approaches calls for scalable and efficient software development. User-friendly implementations will facilitate broad applications.

Supplementary Material

supplement

ACKNOWLEDGEMENTS

This work is supported by grants from the National Science Foundation (SES1760133) and the National Institutes of Health (U01MD017867). The authors are grateful to Dr Michael R. Elliott for his assistance in the theory and implementation of the weighted finite population Bayesian bootstrap.

DATA AVAILABILITY STATEMENT

The simulation code is publicly available on GitHub: https://github.com/likat/EMRP. The data used in the application study are openly available from the New York City Longitudinal Survey of Wellbeing: https://cprc.columbia.edu/content/new-york-city-longitudinal-survey-wellbeing.

Abbreviations:

EMRP embedded multilevel regression and poststratification

LSW longitudinal survey of wellbeing

MRP multilevel regression and poststratification

WFPBB weighted finite population Bayesian bootstrap

FIGURE 1 Conceptual illustration of the embedded multilevel regression and poststratification (EMRP) workflow.

FIGURE 2 The illustration of the weighted finite population Bayesian bootstrap (WFPBB). From the original, or “parent sample” (PS), we generate Bayesian bootstrap samples B1,…,BL, and for each bootstrapped sample we pool F populations drawn from the weighted Pólya to produce a single synthetic population Sl of size F*N.

FIGURE 3 Comparing the simulation cases for the overall and subdomain mean estimates, where subdomains are defined by inclusion probability ranges and either the joint distribution of (Z,X) (left) or the levels of (Z) (right). We compare the root mean squared error (rMSE), absolute bias, average 95% confidence interval (CI) length, and 95% CI coverage rates between the direct imputation of the outcome using the weighted finite population Bayesian bootstrap (WFPBB), classical multilevel regression and poststratification (MRP), and embedded multilevel regression and poststratification (EMRP) methods (WFPBB-MRP, Multinomial MRP, Two-stage MRP). Darker colors correspond to higher values.

FIGURE 4 Simulation results for the population cell frequency Nˆj estimates from the EMRP methods (WFPBB-MRP, Multinomial MRP, Two-stage MRP). We report the root mean squared error (rMSE), absolute bias, average 95% confidence interval (CI) length, and 95% CI coverage rate. Tiles correspond to cells created by the (Z,X) cross-tabulation and are ordered by inclusion probabilities pI: cells with smaller inclusion probabilities are at the bottom of the y-axis, with the lowest in the bottom left; cells with greater inclusion probabilities are at the top, with the highest on the upper right. Darker colors correspond to higher values.

FIGURE 5 Estimates of food insecurity prevalence for the overall NYC adult population and subgroups defined by annual income: <$35k(n=884),$35-55k(n=273),$55-100k(n=456), and >$100k(n=615). We report results from the unweighted and survey-weighted estimates, direct imputation with the weighted finite population Bayesian bootstrap (WFPBB), classical multilevel regression and poststratification (MRP), and embedded multilevel regression and poststratification (EMRP) techniques (WFPBB-MRP, Multinomial MRP, Two-stage MRP). Error bars refer to the 95% confidence intervals.

FIGURE 6 Estimates of food insecurity prevalence for the subgroups defined by interactions of annual income and agency visitor status. We report results from the unweighted and survey-weighted estimates, direct imputation with the weighted finite population Bayesian bootstrap (WFPBB), classical multilevel regression and poststratification (MRP), and embedded multilevel regression and poststratification (EMRP) techniques (WFPBB-MRP, Multinomial MRP, Two-stage MRP). Error bars refer to the 95% confidence intervals.

TABLE 1 Value ranges of the simulated inclusion probabilities pI across cells based on the cross-tabulation of fully observed auxiliary variables.

Cell	1–5	6–20	21–40	41–45	46–50	
Range of pI	(0.01, 0.10)	(0.11, 0.40)	(0.21, 0.60)	(0.51, 0.80)	(0.80, 0.99)	

TABLE 2 Differences of the observed values in one random sample and the population values of Pr(X=1) within the subgroups of two simulation cases: (1) Subgroup membership is defined based on the joint (Z,X) distribution and (2) membership is defined based on categories of (Z) only.

(Z,X)	pI:0.03-0.27	pI:0.21-0.38	pI:0.30-0.54	pI:0.40-0.95	
	0.205	−0.046	−0.013	0.048	
(Z)	pI:0.03-0.25	pI:0.25-0.38	pI:0.34-0.51	pI:0.40-0.93	
	0.122	−0.031	0.003	0.054	

TABLE 3 Simulation results for the overall and subdomain mean estimates, where subdomains are defined by inclusion probability ranges and the joint distribution of (Z,X).

		Overall	pI:0.40-0.95	pI:0.30-0.54	pI:0.21-0.38	pI:0.03-0.27	
Population Y¯		0.577	0.581	0.427	0.584	0.618	
rMSE	Unweighted Est.	0.015	0.036	0.022	0.026	0.035	
	Weighted Est.	0.010	0.007	0.010	0.017	0.053	
	WFPBB	0.010	0.007	0.011	0.017	0.054	
	Classical MRP	0.007	0.033	0.059	0.053	0.068	
	WFPBB-MRP	0.007	0.007	0.011	0.011	0.018	
	Multinomial MRP	0.007	0.007	0.012	0.011	0.018	
	Two-stage MRP	0.007	0.008	0.013	0.014	0.017	
Bias	Unweighted Est.	0.014	0.035	−0.019	−0.019	0.016	
	Weighted Est.	−0.001	<0.001	<0.001	<0.001	−0.007	
	WFPBB	<0.001	−0.001	0.001	−0.001	0.003	
	Classical MRP	−0.002	−0.033	0.058	−0.052	−0.067	
	WFPBB-MRP	0.001	0.004	−0.008	−0.005	0.007	
	Multinomial MRP	0.000	0.004	−0.008	−0.006	0.006	
	Two-stage MRP	−0.002	0.006	−0.010	−0.011	0.005	
95% CI length	Unweighted Est.	0.023	0.035	0.046	0.068	0.104	
	Weighted Est.	0.035	0.038	0.047	0.069	0.164	
	WFPBB	0.045	0.050	0.061	0.088	0.202	
	Classical MRP	0.033	0.033	0.037	0.040	0.062	
	WFPBB-MRP	0.038	0.041	0.049	0.059	0.105	
	Multinomial MRP	0.033	0.035	0.040	0.041	0.073	
	Two-stage MRP	0.033	0.035	0.041	0.042	0.073	
Coverage rate	Unweighted Est.	0.305	0.005	0.640	0.805	0.835	
	Weighted Est.	0.945	0.990	0.975	0.955	0.880	
	WFPBB	0.950	0.995	1.000	0.990	0.940	
	Classical MRP	0.980	0.005	0.000	0.000	0.005	
	WFPBB-MRP	1.000	1.000	0.990	0.995	1.000	
	Multinomial MRP	0.970	0.985	0.885	0.910	0.955	
	Two-stage MRP	0.975	0.975	0.895	0.870	0.985	
Note: We report root mean squared error (rMSE), absolute bias, average 95% confidence interval (CI) length, and 95% CI coverage rate from the direct imputation of the outcome using the weighted finite population Bayesian bootstrap (WFPBB), classical multilevel regression and poststratification (MRP), embedded multilevel regression and poststratification (EMRP) methods (WFPBB-MRP, Multinomial MRP, Two-stage MRP), and the unweighted and survey-weighted estimators.

TABLE 4 Simulation results for the overall and subdomain mean estimates, where subdomains are defined by inclusion probability ranges and the levels of (Z).

		Overall	pI:0.40-0.93	pI:0.34-0.51	pI:0.25-0.38	pI:0.03-0.25	
Population Y¯		0.572	0.507	0.535	0.564	0.660	
rMSE	Unweighted Est.	0.015	0.034	0.014	0.025	0.067	
	Weighted Est.	0.010	0.008	0.014	0.017	0.043	
	WFPBB	0.010	0.008	0.012	0.016	0.046	
	Classical MRP	0.007	0.012	0.009	0.009	0.017	
	WFPBB-MRP	0.007	0.007	0.010	0.011	0.019	
	Multinomial MRP	0.007	0.007	0.010	0.011	0.019	
	Two-stage MRP	0.007	0.011	0.009	0.011	0.017	
Bias	Unweighted Est.	0.014	0.033	<0.001	−0.018	0.063	
	Weighted Est.	<0.001	<0.001	0.002	0.001	−0.001	
	WFPBB	<0.001	<0.001	−0.001	<0.001	0.004	
	Classical MRP	−0.002	0.010	<0.001	0.001	−0.006	
	WFPBB-MRP	0.001	0.005	0.001	−0.004	0.009	
	Multinomial MRP	<0.001	0.004	<0.001	−0.005	0.008	
	Two-stage MRP	−0.002	0.009	−0.001	−0.006	<0.001	
95% CI length	Unweighted Est.	0.023	0.042	0.053	0.063	0.083	
	Weighted Est.	0.035	0.045	0.053	0.064	0.140	
	WFPBB	0.045	0.058	0.069	0.082	0.174	
	Classical MRP	0.033	0.036	0.046	0.045	0.076	
	WFPBB-MRP	0.038	0.046	0.051	0.057	0.104	
	Multinomial MRP	0.033	0.036	0.046	0.045	0.077	
	Two-stage MRP	0.033	0.036	0.046	0.046	0.078	
Coverage rate	Unweighted Est.	0.340	0.055	0.935	0.775	0.170	
	Weighted Est.	0.895	1.000	0.950	0.950	0.890	
	WFPBB	0.950	0.995	0.995	1.000	0.945	
	Classical MRP	0.980	0.935	0.990	0.980	0.975	
	WFPBB-MRP	1.000	1.000	0.995	0.995	0.990	
	Multinomial MRP	0.975	0.995	0.980	0.970	0.955	
	Two-stage MRP	0.975	0.960	0.990	0.975	0.980	
Note: We report root mean squared error (rMSE), absolute bias, average 95% confidence interval (CI) length, and 95% CI coverage rate from the direct imputation of the outcome using the weighted finite population Bayesian bootstrap (WFPBB), classical multilevel regression and poststratification (MRP), embedded multilevel regression and poststratification (EMRP) methods (WFPBB-MRP, Multinomial MRP, Two-stage MRP), and the unweighted and survey-weighted estimators.

TABLE 5 Descriptive summary of food insecurity and sociodemographics for respondents from the longitudinal survey of wellbeing, stratified by whether the respondent has visited a food acquisition agency in the last 12 months (visitor) or not (nonvisitor).

	Visitor (%)	Nonvisitor (%)	Overall (%)	
Sample size	626	1602	2228	
Age				
 18–35	33.9	26.3	28.4	
 36–50	28.0	24.7	25.6	
 50+	38.2	49.1	46.0	
Sex				
 Male	34.8	37.1	36.4	
 Female	65.2	62.9	63.6	
Race				
 White	13.7	37.8	31.1	
 Black	34.8	28.3	30.1	
 Other	51.4	33.9	38.8	
Education				
 Less than high school	22.5	11.0	14.2	
 High school	30.4	20.1	23.0	
 Some college	27.5	22.2	23.7	
 Bachelors or higher	19.6	46.7	39.1	
Income				
 <$35k	32.6	57.7	39.7	
 $35–55k	11.5	14.1	12.3	
 $55–100k	22.3	15.7	20.5	
 >$100k	33.5	12.6	27.6	
Food security status				
 Food secure	79.4	94.0	89.9	
 Food insecure	20.6	6.0	10.1	
Note: Values are reported as percentages.

CONFLICT OF INTEREST STATEMENT

The authors declare no potential conflict of interest.

SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.
==== Refs
REFERENCES

1. Jacobs JA , Jones E , Gabella BA , Spring B , Brownson RC . Tools for implementing an evidence-based approach in public health practice. Prev Chronic Dis. 2012;9 :110324. doi:10.5888/pcd9.110324
2. Si Y , Pillai NS , Gelman A . Bayesian nonparametric weighted sampling inference. Bayesian Anal. 2015;10 (3 ):605–625.
3. Si Y , Trangucci R , Gabry JS , Gelman A . Bayesian hierarchical weighting adjustment and survey inference. Surv Methodol. 2020;46 (2 ):181–214.
4. Si Y On the use of auxiliary variables in multilevel regression and poststratification. Under Review. https://arxiv.org/abs/2011.00360 2022.
5. Covello L , Gelman A , Si Y , Wang S . Routine hospital-based SARS-CoV-2 testing outperforms state-based data in predicting clinical burden. Epidemiology. 2021;32 (6 ):792–799.34432721
6. Si Y , Covello L , Wang S , Covello T , Gelman A . Beyond vaccination rates: a synthetic random proxy metric of Total SARS-CoV-2 immunity seroprevalence in the community. Epidemiology. 2022;33 (4 ):457–464.35394966
7. Holt D , Smith TMF . Post stratification. J R Stat Soc Ser A (General). 1979;142 (1 ):33–46.
8. Gelman A , Carlin JB . Poststratification and weighting adjustments. In: Groves RM , Dillman DA , Eltinge JL , Little RJA , eds. Survey Nonresponse. NY: Wiley; 2002:289–302.
9. Wang W , Rothschild D , Goel S , Gelman A . Forecasting elections with non-representative polls. Int J Forecast. 2015;31 (3 ):980–991.
10. Zhang X , Holt JB , Yun S , Lu H , Greenlund KJ , Croft JB . Validation of multilevel regression and poststratification methodology for small area estimation of health indicators from the behavioral risk factor surveillance system. Am J Epidemiol. 2015;182 (2 ):127–137.25957312
11. Yougov Inc. Introducing the YouGov referendum model. https://yougov.co.uk 2017.
12. Kuriwaki S , Ansolabehere S , Dagonel A , Yamauchi S . The geography of racially polarized voting: calibrating surveys at the district level. Am Polit Sci Rev. 2023;1–18. doi:10.1017/s0003055423000436
13. Ghitza Y , Gelman A . Voter registration databases and MRP: toward the use of large-scale databases in public opinion research. Politic Anal. 2020;28 :507–531.
14. Dong Q , Elliott MR , Raghunathan TE . A nonparametric method to generate synthetic populations to adjust for complex sampling design features. Surv Methodol. 2014;40 (1 ):29–46.29200608
15. Makela S , Si Y , Gelman A . Bayesian inference under cluster sampling with probability proportional to size. Stat Med. 2018;37 (26 ):3849–3868.29974495
16. Reilly C , Gelman A , Katz J . Poststratication without population level information on the poststratifying variable, with application to political polling. J Am Stat Assoc. 2001;96 :1–11.
17. Kastellec JP , Lax JR , Malecki M , Phillips JH . Polarizing the electoral connection: partisan representation in supreme court confirmation politics. J Polit. 2015;77 (3 ):787–804.
18. Wimer C , Garfinkel I , Gelblum M , Poverty Tracker—Monitoring Poverty and Well-Being in NYC. NY: Columbia Population Research Center and Robin Hood Foundation; 2014.
19. Si Y , Gelman A . Survey Weighting for New York Longitudinal Survey on Poverty Measure. Tech. rep NY: Columbia University; 2014.
20. Gundersen C , Ziliak JP . Food insecurity and health outcomes. Health Aff. 2015;34 (11 ):1830–1839.
21. Ghitza Y , Gelman A . Deep interactions with MRP: election turnout and voting patterns among small electoral subgroups. Am J Polit Sci. 2013;57 (3 ):762–776.
22. Downes M , Gurrin LC , English DR , Multilevel regression and poststratification: a modeling approach to estimating population quantities from highly selected survey samples. Am J Epidemiol. 2018;187 (8 ):1780–1790.29635276
23. Ornstein JT . Stacked regression and poststratification. Polit Anal. 2020;28 (2 ):293–301.
24. Little R , Vartivarian S . Does weighting for nonresponse increase the variance of survey means? Surv Methodol. 2005;31 (2 ):161–168.
25. Rubin DB . The Bayesian bootstrap. Ann Stat. 1981;9 (1 ):130–134.
26. Leemann L , Wasserfallen F . Extending the use and prediction precision of subnational public opinion estimation. Am J Polit Sci. 2017;61 (4 ):1003–1022.
27. Deming WE , Stephan FF . On a least squares adjustment of a sampled frequency table when the expected marginal totals are known. Ann Math Stat. 1940;11 (4 ):427–444.
28. Little R , Wu MM . Models for contingency tables with known margins when target and sampled populations differ. J Am Stat Assoc. 1991;86 :87–95.
29. Si Y , Zhou P . Bayes-raking: Bayesian finite population inference with known margins. J Surv Stat Methodol. 2021;9 (4 ):833–855.
30. Meeden G , Lazar R , Geyer CJ . R package polyapost: simulating from the Polya posterior. https://cran.r-project.org/web/packages/polyapost/index.html 2020.
31. Stan Development Team. Stan: a C++ library for probability and sampling. http://mc-stan.org 2021.
32. Gelman A , Jakulin A , Pittau MG , Su YS . A weakly informative default prior distribution for logistic and other regression models. Ann Appl Stat. 2008;2 (4 ):1360–1383.
33. Lumley T Survey: analysis of complex survey samples. https://cran.r-project.org/web/packages/survey/index.html 2021.
34. Vehtari A , Gelman A , Gabry JS . Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput 2017;27 :1413–1432.
35. Gelman A , Meng XL , Stern HS . Posterior predictive assessment of model fitness via realized discrepancies (with discussion). Stat Sin. 1996;6 :733–807.
