Nonprobability samples are quickly emerging owing to the rapidly declining response rates and increasing costs of probability samples, and offer detailed outcomes of interest with large sample sizes that are not available in probability surveys. The selection mechanism may be voluntary or deterministic, and the inclusion probabilities are unknown. Our motivating application, the Adolescent Brain Cognitive Development (ABCD) study aims for national representation but is a nonprobability sample (ABCD, 2018) . The 21 research sites across the U.S. geography are selected for convenience with operational constraints, and the sample enrollment is conditional on the school and parental consents. The ABCD study design can result in selection bias. The lack of randomization and sampling frames demolishes the inferential framework, leading the validity to rely on the quality of auxiliary information and underlying population model specification (Smith, 1983) .
Current approaches for inferences with nonprobability surveys rely on calibration with a reference probability sample or population control information (Elliott and Valliant, 2017) . The design-based approach assumes quasi-randomization and constructs pseudo-weights in the combination of the probability and nonprobability samples. Model-based approaches fit a model of the survey outcome based on the sample and predict the outcome for the nonsampled population units (Ghosh and Meeden, 1997) . Doubly robust (DR) methods combine the two approaches and propose a weighted estimator that is also a function of predicted outcomes (Chen et al., 2019; . DR estimators offer protection against the misspecification of the sample inclusion model or the outcome model. Kang and Schafer (2007) find that many DR methods perform better than simple inverse-probability weighting, however, none of the studied has improved upon the performance of simple regression-based prediction of the nonsampled values.
As a prediction approach to generating synthetic populations, multilevel regression and poststratification (MRP, Gelman and Little (1997) ), originally applied to estimate state-level public opinion from sociodemographic subgroups using sample surveys, has become increasingly popular. MRP is now a standard approach for adjusting selection bias and facilitating small area estimation (Pfeffermann, 2013; Rao and Molina, 2015; Ghosh, 2020) to extrapolate sample inferences to the target population. MRP has two key components: 1) multilevel regression for small area estimation by setting up a predictive model with a large number of covariates and regularizing with Bayesian prior specifications; and 2) poststratification to adjust for selection bias and correct for imbalances in the sample composition. The flexible modeling of survey outcomes can capture complex data structures conditional on poststratification cells, which are determined by the crosstabulation of categorical variables that affect the sample inclusion (selection and response) and use population control information to balance the sample discrepancy (Holt and Smith, 1979; Gelman and Carlin, 2000) .
Besides its widespread and successful applications in social sciences, especially in the U.S. and U.K. election forecasting (e.g., Lax and Phillips (2009a,b) ; Wang et al. (2015) ; Lauderdale et al. (2020) ; Zahorski (2020) ), MRP also demonstrates its stability for small area estimation of disease prevalence and adjustment for nonresponse bias in public health (e.g., Zhang et al. (2015) ; Downes et al. (2018) ; Downes and Carlin (2020) ). Notably, MRP has been actively used for the 2020 U.S. presidential election forecasting (The Economist, 2020; Yougov, 2020) and COVID-19 prevalence estimates (Gelman and Carpenter, 2020; Covello et al., 2020) .
We aim to apply MRP to correct for the selection bias and generate valid and stabilized estimates of child groups with diverse socio-demographic characteristics in the ABCD study. The large-scale survey is designed to assess the target population of the U.S. 9-and 10-year old children with diverse biological, familial, social and environmental factors (Garavan et al., 2018) . Population representativeness is crucial for external validity to generalize the results. Calibrating the ABCD study to eliminate the sample discrepancy, MRP has the potential to yield valid estimates, especially for minority groups with small sample sizes.
However, criticisms on treating MRP "as a widely used gold standard" emerge and argue with substantial variation in the performances (Buttice and Highton, 2013; Valliant, 2019) . MRP is a statistical method of adjusting for selection/nonresponse bias and data sparsity to improve survey estimates for overall population and domain inferences, and could be subject to poorly predictive auxiliary information, model misspecification and invalid assumptions. However, there is no guideline but much confusion on the selection and modeling of auxiliary information in practical applications of MRP. This paper fills in the knowledge gap. We develop the statistical foundation on the use of auxiliary information in MRP.
MRP is deeply rooted in survey methodology. Survey inference is essentially a data integration process that calibrates the sample data with auxiliary information to achieve the population representativeness and valid inferences. Accounting for the inclusion mechanism in prediction modeling, MRP combines design-based and model-based approaches for survey inference, similar to the design-adjusted, model-assisted estimation (Sarndal et al., 1992; Breidt and Opsomer, 2017) . As post-adjustments after data collection, MRP can unify inferences for probability and nonprobability samples under a systematic data integration framework. MRP integrates the sample data with the auxiliary or population control information, from census records or large-scale survey data that have low variability, for example, the American Community Survey (ACS) or the Current Population Survey (CPS). The availability and quality of external information and the integrating method affect the inferential validity. Propagating all sources of uncertainty in the synthetic population generation, we generalize MRP as a unified framework for data integration and robust survey inferences with challenging data settings in practice.
The paper structure is organized as below. In Section 2 we develop the theoretical foundation of MRP in the aspects of poststratification and robust survey inferences and propose a systematic framework for data integration. We use simulation studies to illustrate the improvements of the MRP framework in Section 3 and demonstrate the application of inferences with the ABCD nonprobability survey in Section 4. We summarize existing challenges and potential extensions in Section 5. Rubin (1983) and Little (1983) point out that any model for survey outcomes should condition on all information that predicts inclusion probabilities. Suppose the outcome in the population is Y i , the inclusion indicator is I i , and the auxiliary variables of the population are denoted by X i , for i = 1, . . . , N , where N is the population size. We consider the inference framework (Smith, 1983) 
the validity of which relies on the inclusion mechanism f (I i |X i ) and the outcome model f (Y i |I i , X i ). When the outcome Y i is correlated to the inclusion indicator I i , the inclusion mechanism is informative and has to be accounted for in the analysis of the sampled data f (Y i |I i = 1, X i ). Survey practice often assumes that the inclusion mechanism is ignorable, i.e., missing at random (MAR), given the auxiliary information (Rubin, 1976) . However, this demands correct model specification with rich, highly predictive information X i from integrated data sources. By constructing poststratification cells with discretized auxiliary variables X i from the target population, MRP is a post-collection adjustment method that can make inferences of both probability and nonprobability samples and account for the design and response mechanisms. In the paper, we also refer to the included samples as respondents.
We develop statistical guidelines on the use of such categorical auxiliary information in MRP to achieve inferential validity and balance estimation bias and variance. We first examine the properties of MRP in the aspects of poststratification and robust survey inferences, and then present the proposed systematic MRP framework for data integration and inferences. As popular quantities of interest in survey inference, we focus on the descriptive summaries of the population: the overall mean and subdomain inferences.
We start by estimating the population mean of a single survey response:
Assume the population means in the respondents and nonrespondents areȲ R andȲ M , respectively. The population proportion of respondents is ψ. The overall populationȲ is given byȲ = ψȲ R + (1 − ψ)Ȳ M . Suppose units in the population and the sample can be divided into J poststratification cells with population cell size N j and sample cell size n j for each cell j = 1, . . . , J, with the population size N = J j=1 N j and the sample size n = J j=1 n j . Let Y j be the population mean andȳ j be the sample mean within cell j. The overall mean in the population is Y = J j=1 N j N Y j . For subdomain estimation, since the poststratification cells are constructed at the finest level, we will need to group the set of cell-wise estimates that belong to the subdomain. Let the poststratification cell means for respondents and nonrespondents beȲ jR andȲ jM , respectively, and the population cell proportions of respondents beψ j . The population mean can be expressed as
To make inference aboutȲ with the sample data, the unweighted (UnW) estimator is the average of the sample cell meansȳ
The poststratification (PS) estimator accounts for the population cell sizes as a weighted average of the sample cell means,ȳ
Motivated by the model-based perspective of poststratification (Fay, 1979; Little, 1993) , Gelman and Little (1997) propose MRP by modeling cell estimates and predicting nonresponded cases. Assuming the survey outcome follows a normal distribution with cell-specific mean and variance values,
the proposed MRP estimator for the population mean can be expressed in weighted form,
whereθ j is the model-based estimate ofȲ j in cell j. MRP fits multilevel models to borrow information and smooth cell-wise estimatesθ j , i.e., partial pooling (Gelman and Hill, 2007) . Under a Bayesian paradigm, given an exchangeable prior distribution, θ j ∼ N (µ, σ 2 θ ), where the hyperparameters (µ, σ 2 θ ) are assigned with noninformative prior distributions, the posterior mean estimate for MRP isθ
Given the fitted model, the ratio of sums
is a constant that does not depend on j, and can be approximated by ȳ j * n j /n =ȳ s . This is appropriate if the cell-wise inclusion mechanisms are independent of the group sizes N j and the variance parameters (σ 2 j , σ 2 θ ). Hence, the approximated MRP estimator is given bỹ
as a combined estimator ofȳ ps andȳ s with the shrinkage factor δ j . The shrinkage factor depends on the sample cell size n j , within-cell variance σ 2 j and between-cell variance σ 2 θ , the values of which affect the tendency ofθ mrp toward the unweighted estimateȳ s , as illustrated by the empirical work in Buttice and Highton (2013) . When n j → 0, δ j → ∞, i.e., with small cell sizes, the estimate will be pooled toward the overall unweighted meanθ mrp →ȳ s ; when n j → ∞, δ j → 0, i.e., with large cell sizes, the estimate will be lead toward the PS estimator with cell-wise direct estimates θ mrp →ȳ ps .
Next, we examine the bias and variance in the comparison of the UnW, PS and MRP estimators. Assume that the sample respondents are a random sample of the population respondents: E(ȳ j ) = Y jR , similar to Kalton and Kasprzyk (1988) , and we can calculate the bias
N ψ j denotes the overall proportion of respondents in the population. The term A captures the variation between cell-wise response propensities, and the term B is a weighted average of cell-wise covariances between the response propensities and outcomes. The general conditions that the absolute values |bias(ȳ ps )| < |bias(ȳ s )| are if: 1) A and B have the same sign or 2) |A| > 2|B|. The biases can be approximated in a stochastic model (Bethlehem, 2002) bias
where Cov(ψ, y) is the covariance between the response probabilities and the outcome values, Cov j (ψ, y) is the covariance between the response probabilities and the outcome values within cell j, andψ j is the mean of response probabilities in cell j. The expression bias(ȳ s ) is similar to the data defect index discussed in Meng (2018) . The bias of the MRP estimator is given by
If the dataset is a simple random sample, that is, missing completely at random, ψ j ≡ψ, then all three estimators are unbiased. When the cells are homogeneous with respect to either the response probability or the outcome variable, equivalently, based on the model (3),Ȳ jR =Ȳ jM , Cov j (ψ, y) = 0, and thus, B = 0, the PS estimatorȳ ps is unbiased. However, with the exchangeable prior distribution, the MRP estimatorθ mrp will be biased, the second term of which is nonzero, though the bias will be smaller than bias(ȳ s ).
Conditional on the sample cell sizes n = (n 1 , . . . , n J ), the variance estimates are
where s 2 j is the element variance inside cell j. The conditional variance estimate forθ mrp is approximated by
].
With small n j 's, the variance var(ȳ ps | n) could be large, and var(θ mrp | n, δ) reduces the variance with the shrinkage effect. Under a Bayesian paradigm, the variance estimate ofθ mrp will propagate all sources of uncertainty in the posterior computation. Holt and Smith (1979) show that poststratification based on cells that are homogeneous with respect to the target variable both reduces variance and bias. Poststratification based on cells that are homogeneous with respect to the response probabilities reduce the bias but not necessarily the variance (Little, 1986) . Consistent with Little and Vartivarian (2005) , we recommend incorporating auxiliary variables that are predictive of either the survey outcomes (primarily) or response propensities (secondly) in the construction of poststratification cells, and induce prior distributions to shrink the non-predictive terms to be zero. With data-driven shrinkage modeling to find a tradeoff between bias and variance in the comparison withȳ s andȳ ps , we expect that θ mrp under a flexible model structure and prior distributions yields the smallest root mean squared error (RMSE) and robust inferences.
This highlights the incorporation of cell-wise covariates to introduce a conditionally exchangeable model for θ j given the auxiliary variables X j , for example
Here we use the cell indices j in the notation of auxiliary variables X j , rather than the unit indices i, because the poststratification cells are the resulting cross-tabulation of X i . It is appropriate to assume a conditional exchangeable model after incorporating sufficient relevant information in X j that the cells can be thought of as randomly assigned. A nonexchangeable model would be appropriate if information relevant to the outcome were conveyed in the unit indexes rather than by explanatory variables. With cell-wise covariates X j , the MRP estimator, which is commonly used in practice, becomesθ
as a combined estimator ofȳ ps and J j=1 N j N X j β. The MRP estimator is a population average of the composite estimators that are known in the small area estimation (Ghosh, 2020) . Lahiri and Mukherjee (2007) have examined the design-consistency property of a hierarchical Bayes estimator of a finite population stratum mean when the sample size is large. Specifically, they prove that with an exponential family distribution for the outcome and a normal prior distribution for θ j , the posterior mean estimatorθ j →ȳ j as n j → ∞ and n j /N j → f j for some 0 < f j < 1, and the corrected estimatorθ j − (ȳ j −ȳ j,ds ) is design-consistent, whereȳ j,ds is any design-consistent estimator ofȲ j (Theorem 3.1, Lahiri and Mukherjee (2007) ). Extending to MRP, we have the corollary below.
Corollary 1. Assume the following regularity conditions. (R.1) The survey outcome y ij follows an exponential family distribution with cell-specific mean θ j and variance parameters σ 2 j .
The poststratification cell structure fully accounts for the design information. Then the MRP estimatorθ
Here one of the key regularity conditions is the availability of auxiliary information that constructs the poststratification cells and fully accounts for the design information and nonresponse mechanism. Design-consistency is a desirable property in the randomization approach to finite population sampling. Moreover, MRP improves estimation efficiency and balances bias and variance for robust inferences. Baker et al. (2013) emphasize that inference for any probability or nonprobability survey requires reliance on modeling assumptions, a coherent framework and an accompanying set of measures for evaluating the quality. MRP improves overall population inferences and small area estimation by adjusting for selection/nonresponse bias into hierarchical modeling that is robust against misspecification. The crucial guidelines for generalizable and valid inferences of MRP on the use of auxiliary information include 1) poststratification with auxiliary information, and 2) flexible modeling strategies.
Poststratification with auxiliary information: With a rich set of highly predictive covariates, Wang et al. (2015) have applied MRP to obtain estimates of voting behavior in the 2012 U.S. Presidential election based on a non-representative poll of 350,000 Xbox users and calibration with the auxiliary information from the exit poll data. The cross-tabulation of covariates about voting behavior-sex, race, age, education, state, party ID, political ideology, and reported 2008 vote-results in 176,256 poststratification cells, across which a multilevel logistic regression is fit. Next, the exit poll data with 101,638 respondents from the 2008 presidential election are used as the reference sample for poststratification to adjust for the discrepancies in the sample decompositions of the Xbox data to match the population control information. Their findings show that MRP estimates are in line with the forecasts from leading poll analysts, which were based on aggregating hundreds of traditional polls conducted during the election cycle, and non-representative samples have the potential for population-based inferences. The success of MRP mostly comes from the adjustment of selection bias to which vote swings are mostly attributed (Gelman et al., 2016) .
Current approaches for statistical inferences with nonprobability survey samples assume the relevant auxiliary information is available from a reference probability survey sample (Elliott and Valliant, 2017; Chen et al., 2019; Kim and Tam, 2020) . Incorporating auxiliary information as covariates in response propensity modeling, inverse propensity score weighting (IPW) estimators are usually used but are sensitive to misspecification, especially with a wide range of values, and can result in highly variable estimates, for example, discussed by Tan (2007) . Flexible predictive models are fit to construct stabilized pseudo-weights for nonprobability samples, such as Bayesian Additive Regression Trees (BART, Rafei et al. (2020) ) and kernel weighting approaches . In contrast to survey weighting, MRP constructs poststratification cells based on the auxiliary information and groups individuals to reduce the variability that is also stabilized by the multilevel outcome modeling.
Flexible modeling strategies: DR estimators improve the efficiency and the robustness of IPW estimators by combining a prediction model for the survey outcome (Bang and Robins, 2005) . Kang and Schafer (2007) show that the regression prediction estimator outperforms DR estimators with a predictive model. Moreover, most DR estimators do not apply to subgroup estimation. With a linear regression model, the well-known general regression (GREG) estimator improves estimation accuracy in the combination with poststratification adjustments (Deville and Särndal, 1992) . MRP replaces the linear regression in GREG with a multilevel model with smoothing effects across poststratification cells. Lauderdale et al. (2020) consider three pre-election polling applications and conclude that careful model specification is essential in applying MRP. Flexible models of survey outcomes are essential to yield robust model-based inferences. The election forecasting model applied to the Xbox survey is a Bayesian hierarchical model with only main effects. The existing examples of flexible models under MRP include hierarchical models with high-order interactions and global-local shrinkage prior specifications (Ghitza and Gelman, 2013; , Gaussian process (GP) regression models (Si et al., 2015) , BART (Bisbee, 2019) , and stacked regression (Ornstein, 2020) . The flexible models under MRP improve small area estimation and facilitate subdomain and overall population inferences. In addition to a careful selection of predictive auxiliary variables and their high-order interaction terms, we recommend incorporation a flexible function form of inclusion probabilities into the outcome modeling. Next we develop a systematic data integration framework under MRP.
Even with numerous successful applications, MRP is faced with challenges in practical settings. Most MRP applications integrate the probability/nonprobability sample of interest with large probability samples (e.g., ACS, CPS, or census records) and ignore the sampling uncertainty of the latter. That is, MRP treats the population cell counts N j 's from external data as fixed. In practice, the reference probability sample may be small, shown in Figure 1 . We have to estimate population information based on the reference probability sample and then use the estimates for calibration or imputation of outcome values for the nonsampled units. Often the joint population distributions of all auxiliary information are incomplete, e.g., with the joint distribution of only a subset of variables or only marginal distributions being available.
With a large number of auxiliary variables, the number of poststratification cells J could be large, for example, J = 2 × 4 × 4 × 5 × 50 = 8000, in the adjustment of sex (2 levels), race/ethnicity (4 levels), education (4 levels), age (5 levels) and state (50 levels). The resulting sample poststratification cells could be sparse and even empty, as the computational bottleneck for MRP. Little (1993) recommends collapsing small cells to reduce the variance, with a payoff of increased bias. Either 1) when the inclusion probabilities of the individuals are the same within cells; or 2) when the included individuals are similar to those excluded within cells, i.e., conditional MAR inside cells, the poststratification and MRP estimators will be unbiased. Combing cells may violate such assumptions, however, the data sparsity requires shrinkage to stabilize inferences. The Bayesian paradigm of MRP allows the data to determine the pooling effects. Figure 1 : Illustrative systematic data integration framework of combing a small probability (prob) sample and a nonprobability (Non-prob) sample, where they share a common set of covariates X and the survey outcome Y is only collected in Non-prob samples. We create an inclusion indicator I for the nonprobability samples. The probability sample will be used to estimate the population information for the nonsampled units.
We would like to generate synthetic population information if unknown and fit MRP with a flexbile model across cells and Bayesian shrinkage prior distributions, as a systematic data integration framework in Figure 1 on the use of auxiliary information. The general MRP estimator for the population mean would beθ
where bothN j andθ j are model-based estimates.
The MRP under the systematic data integration framework combines the quasi-randomization approach and the superpopulation modeling approach by constructing three main models.
1. The cell sizes: assume that the respondents within poststratification cell j are independently and quasi-randomly sampled from the population cell cases, (n 1 , . . . , n J | ψ) ∼ Multinomial ((cN 1 ψ 1 , . . . , cN J ψ J ), n),
where c is a normalizing constant, c = 1/ J j=1 N j ψ j . This can be approximated by J Poisson distributions when J is sufficiently large or the inclusion probability ψ j is sufficiently small.
Here the population cell counts N j 's are estimated from the reference sample. Often the relevant auxiliary information is partially complete in the reference sample and requires modeling to generate the synthetic population distribution (Reilly et al., 2001; Kastellec et al., 2015) . Ghitza and Gelman (2020) turn to large-scale voter registration databases to accommodate the census or exit polls for a more accurate picture of the target population with powerful and stable political covariates. A variety of nonparametric and parametric approaches are available to estimate N j . With survey weights available from the reference probability sample, weighted Bayesian bootstrap approaches have been developed to obtain the population distribution of X j (Zangeneh and Little, 2012; Dong et al., 2014; Zhou et al., 2016; Makela et al., 2018) . Si et al. (2015) have used survey weights in probability samples by assuming ψ j = 1/w j , where w j is the unique weight value in cell j and model (8) to yield smoothed estimates of N j 's with sparse cells even only containing one unit. Leemann and Wasserfallen (2017) combine MRP with raking when only population margins of the auxiliary variables are known. The raking model assumes that the auxiliary variables have independent effects on the inclusion, equivalent to the prediction model of inclusion probabilities below but with only the main effects of auxiliary models (Little and Wu, 1991; Si and Zhou, 2020) .
2. The inclusion probabilities: we can assume that they are concentrated around the average inclusion probability of the sample, for example, a Beta distribution with a mean of n/N . Alternatively, with covariates Z j that affect inclusion propensities and a link function g(·) (e.g., logit), we model the inclusion probabilities,
where Z j 's denote cell-wise covariates and γ represents the model parameters.
3. The outcome: within poststratification cells the units are included with the same probability and independently distributed. We assume that the outcome depends on the cell inclusion probabilities and follows a normal distribution with cell-specific mean and variance values,
where f (ψ j ) is a function of ψ j , and X j denotes the cell-wise covariates that predict the outcome variable, may include both main effects and high-order interaction terms, and overlap with Z j . For notational simplicity, we ignore individual-level covariates, which can be predictive, even continuous covariates, and added in the outcome model, though used in the poststratification. Flexible prior distributions can be introduced on f (ψ j ), such as penalized spline functions and nonparametric Bayesian distributions, to account for potential dependency structure and smoothing effects across cells. Si et al. (2015) induce a GP prior distribution to the mean function, such as f (ψ j ) ∼ GP (log(ψ j )β, Σ( ψ)), where the mean log(ψ j )β is a linear function of the logarithm of the inclusion probability, and the covariance matrix Σ( ψ) depends on the differences between ψ j 's, e.g., Cov(ψ j , ψ j ) = (1 − g)τ 2 exp(− (ψ j −ψ j ) 2 l 2 ) + gτ 2 I j=j with weakly informative or noninformative prior specifications for hyperparameters (g, τ, l). The smoothing function f (ψ j ) can be pre-specified or approximate with basis expansion functions if the number of cells is large. With rich information X j and correct specification of the mean structure X j β, the role of the function f (ψ j ) could be minimal. We recommend to include the flexible specification of f (ψ j ) that offers protection against model misspecification, similar to the doubly robust penalized spline of propensity prediction methods (Little and An, 2004; Breidt et al., 2005; Zhang and Little, 2009 ).
The computation for the models specified in (8), (9), and (10) can be implemented in a fully Bayesian procedure that integrates all sources of uncertainty. Stan makes the computation of MRP generally accessible. As a state-of-the-art platform for statistical modeling and high-performance statistical computation, Stan can perform fully Bayesian analyses with Markov chain Monte Carlo (MCMC) computation. The implementation of MRP is straightforward with the publicly available R packages such as Rstan (Stan Development Team, 2020), and Rstanarm (Goodrich and Gabry, 2020). For example, in Rstanarm, the function stan glmer fits a multilevel model and the function posterior predict imputes the outcome for all nonsampled units in the population. The imputation step draws posterior predictive samples and takes into account all sources of uncertainty, which generates multiple synthetic populations and is implicitly poststratification. We can obtain subgroup estimates by extracting the imputed values for the corresponding domains. We implement the proposed systematic data integration framework in Stan, which has not been available in the public R packages yet.
We conduct simulations to first evaluate the bias and variance trade-off of MRP with varying associations of the poststratification cell structure with the survey outcome and inclusion mechanism. Then we implement the systematic data integration framework with cell-wise covariates and illustrate the improvement of MRP on robust inferences. We compare the UnW estimatorȳ s , the PS estimatorȳ ps and the MRP estimatorθ mrp in inferences for the overall finite population and subdomain mean estimates.
To evaluate the randomness properties, we perform repeated sampling of 100 samples to assess the bias, standard error (SE), RMSE and nominal coverage rates (CR) of 95% confidence intervals. We simulate a population of 1000000 units with J = 10 cells of equal sizes. We are interested in the overall population and five subgroup mean estimates, where the gth subgroup composites g + 1 number of cells randomly selected from the total 10 cells, for g = 1, . . . , 5, the goal of which is to generate subgroups with various sample sizes. These six quantities of interest are denoted as All, G1, G2, G3, G4, G5 in the outputs shown in Figure 2 . For each repetition, we randomly draw n = 400 samples based on the pre-specified outcome modeling and inclusion mechanisms in Table 1 . For the MRP specification we assume θ j = β 0 + β 1 j with β 1 ∼ Cauchy(0, 3), a Cauchy prior distribution, and assign noninformative prior distributions to other parameters. This facilitates computation performances under repeated sampling, and can be replaced by the flexible outcome model in (10). 
Outcome Inclusion High (Y, X); high (I, X) Y ∼ N (−5 + 5j, 7 2 ) ψ j = 0.1 + 0.8(j − 1)/9 High (Y, X); low (I, X) Y ∼ N (−5 + 5j, 7 2 ) ψ j = 0.5 Low (Y, X); high (I, X) Y ∼ N (−5, 25 2 ) ψ j = 0.1 + 0.8(j − 1)/9 Low (Y, X); low (I, X) Y ∼ N (−5, 25 2 ) ψ j = 0.5 We consider four different simulation scenarios by varying the strength of association of the auxiliary variables with inclusion probabilities (I, X) and outcome (Y, X) from high to low. The four cases are described in Table 1 . Generally, across all cases, MRP is the most accurate approach with the smallest SEs and RMSEs, especially for the subgroups with small sample sizes. The top two panels in Figure 2 show that when the association between the outcome and the auxiliary variables used for the poststratification cell construction is high, the unweighted estimator is biased and has low coverage rates, and the MRP and PS estimators yield valid inferences. The findings hold the same when we change both of the associations to be medium.
When the outcome is weakly related to the cell structure, as the lower two panels in Figure 2 , the PS estimator has large variability and potential low coverage, while the MRP estimator yields competitive coverage rates with small SEs. When both associations are low, among three estimators, MRP yields the smallest SE and RMSE with comparable bias values. Interestingly, the UnW and PS estimators have similar performances when the cell structure is weakly related with both the outcome and inclusion mechanism. This may be due to substantial group sizes. In all, with a weakly informative prior distribution, the MRP estimator outperforms the unweighted and poststratification estimators with valid inferences and a bias-variance tradeoff, especially for small subgroups.
We use the ACS 2011-2015 data as auxiliary covariates and simulate the outcome variable to conduct the simulation study. Consider three ACS variables: sex (2 levels: female, male), race/ethnicity (race, 8 levels: white, black, Hispanic, Asian, AIAN, NHPI, Other, Multiple), and family income (inc, 6 levels: 2 2 ) , and the inclusion indicator I i ∼ Bernoulli(p i ) with p i = logit −1 (α 0 + α sex sex i + α race race i + α inc inc i ), where the coefficients β var 's and α var 's are vectors if the covariate var has more than 2 levels, β var 's are random draws between 1 and 10, and α var 's are randomly selected from {−1, 0, 1}, var ∈ {sex, race, inc}. The simulated outcome depends on main effects of the three covariates, and the inclusion mechanism depends on only partial main effects, fewer than those in the outcome model. The resulting inclusion probabilities range from 0.27 to 0.98 with a mean value of 0.75.
For inference of each simulated sample, we implement the systematic data integration framework given by models (8), (9), and (10). We assume f (ψ j ) ∼ GP (0, Σ( ψ)), where the covariance matrix Σ( ψ) has a pre-specified kernel function Cov(ψ j , ψ j ) = (1 − g)exp(−(ψ j − ψ j ) 2 ) + gI j=j and a hyperparameter g ∈ (0, 1) that is given a noninformative prior distribution. The hierarchical prior specification of the coefficients is: β var ∼ N (0, σ 2 y:var ), σ y:var ∼ N + (0, 3 2 ), α var ∼ N (0, σ 2 i:var ), σ i:var ∼ N + (0, 1 2 ), where N + (·) denotes a half-normal distribution restricted to the positive values.
We perform repeated sampling 100 times and set the average sample size as 10000. Even with this large sample size, not all 96 cells will be available in the sample. MRP predicts cell-wise estimates for all the population cells. As inferences with large subgroups, Table 2 presents the results of the systematic MRP integration framework for the overall population and five subgroup mean estimates, where the subgroups compose different numbers of randomly selected cells with varying sizes. Comparing the three estimators, MRP and poststratification generate unbiased and valid estimates, while the unweighted estimators are biased with low nominal coverage rates. Since the outcome depends on all three variables, the PS estimator is appropriate, and the MRP framework induces a flexible model that captures the true structure with slightly increased variability and conservative confidence coverages, mainly due to the nonparametric Bayesian GP prior specification.
In the assessment of the cell-wise mean estimates, we randomly draw one sample, and the selected sample covers 80 cells, 11 of which only have 1 unit, and 23 cells have fewer than 5 cases. MRP generates smaller bias and RMSE values than the PS estimator. Figure 4 in the Appendix summarizes the relative bias and RMSE of the sampled 80 cell mean estimates.
To evaluate the performances when the auxiliary information becomes redundant with a sparse cell structure, we consider another simulation scenario with the ACS auxiliary variables. Suppose that the outcome only depends on sex and race: Y i = β 0 + β sex sex i + β race race i + i , i ∼ N (0, 2 2 ), and the inclusion indicator depends on three variables with probabilities p i = logit −1 (α 0 +α sex sex i + α race race i + α inc inc i ), where the coefficients β var 's are random draws between 1 and 10, and α var 's are randomly selected from the range (−2.5, 1.5) with a break of 0.5. The resulting inclusion probabilities are generally low, ranging from 0 to 0.38 with a mean value of 0.04. The randomly drawn sample has 9939 cases and 72 cells, 19 of which have fewer than 5 units. We implement MRP with the same setup as described above. Figure 3 depicts the output for the 72 cell mean estimates under this sparse setting. MRP has substantially smaller biases and RMSEs than the PS estimator. This provides further supporting evidence that MRP outperforms the PS estimator and improves small area estimation with efficiency and accuracy, especially with a sparse data structure. 
We apply MRP to make inferences of our motivating nonprobability survey, the ABCD study. The ABCD study is a prospective cohort study and has collected social, health, imaging and genetics measures of 11,875 children aged 9-10 for environmental exposure, neuroimaging and substance use analysis from 21 U.S. research sites between 2016 and 2018. The ABCD sampling and recruitment process is designated to yield an overall sample that closely approximates national sociodemographics of the targeted U.S. children aged 9-10 ( Garavan et al., 2018) . However, the 21 research sites are not randomly chosen based on reasons of convenience such as neuroimaging resource allocation and accessibility, a form of nonprobability-based selection, and the child enrollment is conditional on the school and parental consents, resulting in potential selection and nonresponse bias.
Comparing with the ACS 2011-2015 data, Table 3 shows that the ABCD study oversamples 9-year old children, males, high-income families, and certain race/ethnicity groups, with slightly more representation of children from families with married couples and employment. By design, the race/ethnicity composition for major classes (White, Black, Hispanic) matches the ACS targets fairly closely with children of Asian ancestry being underrepresented and children with self-reported other (AIAN, NHPI, Multiple) race /ethnicity being overrepresented relative to the U.S. population of 9-and 10-year olds. This illustrates sociodemographic discrepancies from the U.S. population and potential selection bias analyzing data (Compton et al., 2019) . Heeringa and Berglund (2019) have constructed weights to match the ABCD sample to the ACS population control information by predicting pseudo-probabilities of sample inclusion and performing raking adjustments. Specifically, a multiple logistic regression model was fit to the concatenated ACS and ABCD data with covariates: age, sex, race/ethnicity, family income, family type, household size, patients' labor force status, and census region, and predict the pseudo-probabilities of inclusion in the ABCD sample, the inverse of which are treated as the initial weights. The initial weights are trimmed at the 2% and 98% quantiles of the distribution and calibrated with raking adjustments of age, sex and race/ethnicity, whose marginal distributions are matched to the ACS. However, the final weights are still widely spread and can result in unstable variances.
In contrast to weighting, we apply MRP to adjust for the sample discrepancies in the estimation of average cognition test scores for the overall U.S. population of 9/10-year olds and diverse sociodemographic groups of interest. The score is a total composite score of cognition based on the NIH Toolbox (ABCD, 2019). We use seven auxiliary variables in Table 3 to construct poststratification cells. The cross-tabulation results in 4800 (= 2 * 2 * 5 * 2 * 5 * 6 * 4) cells, 1517 of which are available in the ABCD data and 3128 cells are available in the ACS data. There are 3 cells in the ABCD that are not available in the ACS. In the ABCD data, 962 cells have fewer than or equal to 5 units, and 333 cells only include 1 unit.
We integrate these two datasets by fitting a model to the ABCD and then predicting the potential outcomes of the ACS dataset (in the adjustment of the ACS weights), similar to Figure 1 , but with a large probability sample of the ACS. The outcome model incorporates the main effects of seven auxiliary variables
Here we use dummy indicators for age (9, 10), sex (fem: female, male), and family type (mar: married, other), and multiple terms indicating levels for race/ethnicity, family income (inc), household size (hhsize), and family labor force status (lf), which are assigned with hierarchical prior distributions: α var j ∼ N (0, σ 2 var ), σ var ∼ N + (0, 1 2 ), var ∈ {race, inc, hhsize, lf }. Assume the error term inside cells follows a Cauchy distribution i ∼ Cauchy(0, 1). Our analysis results are not sensitive to the hyperparameter specification and generate the same findings when we increase the values (e.g., N + (0, 3 2 ), Cauchy(0, 3)).
We perform the MCMC computation in Stan with two chains of 2000 iterations each. The diagnostic measuresR have values around 1 and indicate convergence. The covariates age, sex, race/ethnicity, family income, household size and labor force status are significant predictors of the test score, while the marital status is not significant (α married = −0.85, with 95% CI (-1.96, 0.67)). The 10-year olds tend to have a higher score than the 9-year olds (α age = 4.61, with 95% CI (4.54, 4.67)), and the girls have higher scores than boys (α f em = 0.96, with 95% CI (0.91, 1.01)). The posterior mean estimates of the variance parameters are significantly different from 0,σ race = 2.27 (95% CI: 1.57, 3.33),σ inc = 2.16 (95% CI: 1.50, 3.17),σ hhsize = 0.78 (95% CI: 0.45, 1.87), and σ lf = 0.69 (95% CI: 0.31, 1.80), showing the covariates are predictive for the cognition outcome.
Since these covariates have different distributions between ABCD and ACS, we expect that the poststratification adjustment to the ACS will generate different results from the unadjusted ABCD sample analysis. Table 4 presents the mean estimates of cognition test scores and 95% CIs for the overall population and four sociodemographic subgroups of U.S. 9/10-year-old children. We can predict the outcome values for the 3128 cells in the ACS and also all the 4800 cells, and we have the population counts for the ACS cells. For the overall mean score estimates, the MRP estimates based on the 3128 cells are 86 (95% CI: 85.94, 86.07), which are different from the MRP estimate based on the 1517 cells: 85.32 (95% CI: 85.23, 85.4). For comparison, we present the MRP estimates based on the cells that are also available in the ABCD in Table 4 . We also apply the pseudo-weights of the ABCD baseline survey and compute the weighted estimates. Based on Table 4 , for the overall mean score estimates, the MRP estimate (85.32) is different from the weighted estimate (85.73), and both are significantly lower than the sample estimate without adjustments (86.2). MRP generates lower scores for girls from married families, low-income white children, and those whose parents are married and employed with larger households, but slightly higher scores for black children, compared to both the weighted and sample estimates. The variance of MRP estimators is lower than that of the weighted and unweighted sample estimators, an illustration of efficiency gains.
We have shown that cognitive performances vary across child groups with diverse sociodemographic and familial characteristics, and the adjustment of auxiliary variables with the ACS study yields substantially different results. MRP uses the highly predictive auxiliary information of the cognitive assessments and adjusts for the differential response propensities. The MRP estimates compose model-based predictions weighted by population cell counts and have efficiency gains. The weighted estimates use the raw values that are subject to measurement error. The poststratification of MRP borrows the joint distribution of the auxiliary variables from ACS, however, the weighting adjustment only uses their marginal distributions. It is possible that the joint distributions of the auxiliary variables are substantially different between the population and the sample, and the high-order interaction terms affect the sample inclusion propensities, leading to different estimates. Here, we incorporate auxiliary information into MRP and stress the potentially substantial impacts on inferences from the sample non-representation. The external validity assessment of different findings requires additional information and substantive knowledge, which will be discussed below.
Large-scale nonprobability surveys quickly emerge and demand qualified auxiliary information and robust statistical adjustments to achieve representative inferences. Integrating nonprobability samples with probability surveys, MRP adjusts for selection/nonresponse bias and data sparsity to improve survey estimates for population inferences, especially small area estimation. MRP constructs poststratification cells, fits flexible hierarchical models, and pools cell estimates weighted by the population cell counts as weights. The inference validity relies on the poststratification structure and model specification. The key to success is the availability of highly predictive auxiliary variables.
MRP assumes quasi-randomization given the auxiliary variables, where the inclusion probabilities are treated as equal inside cells. MRP fits a superpopulation model to predict outcome values of nonsampled units. This article develops guidelines on the use of auxiliary information in MRP and a systematic framework for data integration and inference. The systematic framework propagates all sources of uncertainty, including the estimation of unknown population control information. We have demonstrated that MRP can achieve a balance between robustness and efficiency.
This comprehensive study of MRP opens up several interesting future extensions. First, a sensitivity analysis framework will offer insights into potential bias if the sample inclusion is nonignorable, still correlated with the survey outcomes conditional on the available auxiliary information. Little et al. (2019) have proposed a measure of the degree of departure from ignorable sample selection. Depending on the magnitude of the correlation structure, a sensitivity index can be developed to reflect the range of potential bias values. Second, integrating multiple datasets, beyond two samples, can supplement the list of predictive auxiliary variables that can be sequentially imputed, and also help access the inferential validity. Data integration techniques learn about the common target population structure while accounting for the heterogeneity of different studies. Existing organic databases and administrative records could provide unprecedented information. Linking auxiliary information across multiple data sources has become a research priority for most statistical agencies. Third, MRP focuses on the modeling of a single survey outcome and can extend to multivariate outcomes. The selection of auxiliary variables will become the union of those predictive for any outcome. Weighting cells can be constructed based on predictive mean matching and propensity score subgrouping, to achieve double robustness. Fourth, MRP implicitly generates multiple synthetic populations according to the specified model and brings in methodological and practical challenges in model evaluation. The synthetic population prediction can utilize state-ofthe-art machine learning and deep learning algorithms to improve accuracy. The evaluation should include goodness of fit in the sample data and achieved representation of the target population. Model evaluation in terms of generalizability is challenging in practical application studies (Keiding and Louis, 2016) . Election forecasting can be checked with actual results, however, a gold standard does not exist in most areas. This is an area that needs further developments and collaborations across multidisciplines.
Qualified research emphasizes both sound design and analysis approaches. While soliciting the predictive auxiliary information for analysis, improving the data collection process is essential. MRP is a useful remedy after data collection. Building an integrated database combining multiple data sources to provide big data for a systematic design and analysis process and promote open science, is another promising direction.
In this application, we focus on the cognition assessment at baseline in the comparison of sociodemographic groups in the ABCD study and have demonstrated the use of auxiliary information in MRP to make representative inferences. The ABCD study is a longitudinal study and advocates "population neuroscience" for child development studies. Our future work also includes nonresponse bias adjustment due to attrition and population-based inferences with the brain and neurological outcome measures and methodology developments combining MRP with image modeling.
Title: Figure 