Improving Survey Inference Using Administrative Records Without Releasing Individual‐Level Continuous Data
Introduction Probability samples play an important role in survey research, facilitating inference about health measures in large finite populations using moderately sized samples. However, probability surveys have suffered increasingly high nonresponse rates in the past several decades [ , , , , ], which makes these probability surveys often nonrepresentative, challenging the validity of survey inference. Auxiliary information about the target population can be used to improve survey inference . Such data can be obtained from sources such as administrative records. When only discrete auxiliary variables are available, poststratification or raking can be used to weight the sample to reduce bias in the survey estimation [ , , , ]. Poststratification requires the joint population distributions of the auxiliary variables to be known, whereas raking utilizes the more commonly available marginal population distributions . Although poststratification and raking are easy to implement, the resultant weights can be highly variable, and consequently, the weighted estimators can be unstable [ , , , , ]. Predictive inference represents an alternative framework for correcting nonresponse bias, offering the advantage of improved efficiency in survey inference [ , , , , , ]. Multilevel regression and poststratification (MRP) is a widely used model‐based alternative to the described weighted estimation approaches with important applications in social and political sciences . MRP models survey outcomes using a multilevel regression and can generate accurate survey estimates even from nonprobability samples . Although MRP has an advantage in improving efficiency in survey estimation over poststratification and raking, it requires all the auxiliary variables to be discrete. When administrative records contain individual‐level continuous auxiliary variables, MRP models cannot be directly applied. Instead, prediction models that allow flexible associations between survey outcomes and these continuous auxiliary variables are desirable. Liu et al. proposed a regularized prediction approach using soft Bayesian Additive Regression Trees (BART) , which predicts survey outcomes in the population using detailed individual‐level data from administrative records. This approach has been shown to effectively reduce bias and improve efficiency in survey estimates. However, individual‐level continuous administrative records are often inaccessible due to confidentiality concerns . With discrete auxiliary variables, only frequency tables of the population data need to be released to the public to allow weighting or model‐based adjustments. In contrast, continuous auxiliary variables require individual‐level population information to facilitate prediction or weighting. When confidentiality concerns prevent the release of individual‐level data, continuous auxiliary variables are often discretized using percentiles or meaningful cutoffs before making the population data available. Using discretized versions of continuous variables may lead to efficiency loss in predictive survey inference, particularly when there is a strong relationship between the survey outcome and the continuous variables. Therefore, we seek a data analysis strategy that makes the best use of continuous information while avoiding the need to release individual‐level continuous variables. In this paper, we propose a two‐step strategy for predictive survey inference that maximizes the use of continuous data information in the population while controlling disclosure risk. In the first step, statistical agencies with access to confidential administrative records link survey data with administrative records and estimate inclusion propensities using both continuous and discrete auxiliary variables. They then create a modified population data set by replacing individual‐level continuous data with continuous inclusion propensity scores and discretized versions of the continuous variables. In the second step, data users, who do not have access to the confidential records, perform predictive survey inference using the survey sample and the modified population data. This division of tasks ensures that statistical agencies handle the confidential data while data users conduct survey inference with the modified population data. We conduct simulation studies to evaluate the performance of this two‐step strategy compared with alternative methods. We also illustrate this two‐step strategy using the Ohio Army National Guard Mental Health Initiative (OHARNG‐MHI) Survey, estimating the percentage of lifetime alcohol abuse among all service members in the Ohio National Guard in 2008. Our R package, AuxSurvey , is available on Github: https://github.com/zjg540066169/AuxSurvey and provides researchers with a user‐friendly interface for conducting analyses we discuss in this paper. 

Methods 2.1 Notation and Background We consider a target population consisting of N units with survey outcome Y . Let Z = ( Z 1 , … , Z p ) T denote p discrete auxiliary variables and X = ( X 1 , … , X q ) T be q continuous auxiliary variables, known for the population. For simplicity, we let q = 1 . The continuous variable X is usually not publicly available due to confidentiality concerns. Instead, before releasing to the public, X is discretized, denoted with X ∗ , using percentiles or meaningful cutoffs. The population can be partitioned into J disjoint and exhaustive cells or poststrata defined by the joint distribution of ( Z , X ∗ ) , with N j units in cell j , j = 1 … J and ∑ j = 1 J N j = N where N j > 0 . With Y ‾ j = ∑ i = 1 N j y i / N j , the population mean of Y within cell j , the overall population mean Y ‾ can be written as, (1) θ = Y ‾ = ∑ i = 1 N y i N = ∑ j = 1 J N j Y ‾ j ∑ j = 1 J N j Let s denote a probability survey sample of size n selected from the population with survey outcome values, y 1 , … , y n . We can then use ( Z , X ∗ ) to divide the sample into J cells; the corresponding sample size in cell j is n j with ∑ j = 1 J n j = n and the sample mean of Y in cell j is y ‾ j = ∑ i = 1 n j y i / n j . Assuming that y ‾ j is an unbiased estimate of Y ‾ j , the poststratification estimator of θ can be written as (2) θ ^ = ∑ j = 1 J N j y ‾ j ∑ j = 1 J N j = ∑ i = 1 n w i y i ∑ i = 1 n w i where w i = N j / n j for sample unit i in cell j , j = 1 , … , J . Alternatively, w i in formula ( ) can be created using the raking method when only population margins of ( Z , X ∗ ) are available. Raking weights are obtained via an iterative proportional fitting procedure that begins by adjusting design weights to the marginal distribution of the first auxiliary variable. These adjusted weights are then updated to conform to the marginal distribution of the second auxiliary variable. This process of updating the adjusted weights is carried out for each auxiliary variable and the first iteration ends when the weights are adjusted using the last auxiliary variable. Subsequent iterations are performed until the weights conform to the marginal distributions of all the auxiliary variables, that is, the algorithm converges . The Newton–Raphson method can also be used to obtain the same set of weights . When ( Z , X ∗ ) are related to survey outcome Y and the distributions of ( Z , X ∗ ) differ between sample and population, the poststratification and raking methods can effectively reduce bias in the estimation of θ , by assigning a weight w i for each sample unit. However, the weights can be highly variable, especially when p +q is large and n j is small for some j , j = 1 , … , J . Consequently, the weighted estimator can be unstable. Alternatively, model‐based approaches that predict the nonsampled units in the population based on a prediction model of the sample units can be used to improve efficiency in the θ estimation. The MRP model  has been widely used as a model‐based alternative to poststratification. MRP uses a hierarchical regression model to model the survey outcome of interest given the auxiliary variables, (3) g ( E ( Y i | Z i , X i ∗ ) ) = α 0 + ∑ k = 1 p 1 α k Z i k + ∑ l = p 1 + 1 p α m [ i ] Z l + α m [ i ] x ∗ where g ( · ) is a link function, such as using an identity link for continuous Y and a logit link for binary Y ; α 0 is the intercept; α k is the slope associated with the binary auxiliary variable Z k , k = 1 , … , p 1 and p 1 ≤ p ; and α m [ i ] Z l and α m [ i ] x ∗ correspond to varying coefficients associated with polytomous variables Z l , l = ( p 1 + 1 ) , … , p , and the discretized continuous variable X ∗ , respectively, where m [ i ] indicates the category of the polytomous variable to which the i th unit belongs. The varying coefficients α m [ i ] Z l and α m [ i ] x ∗ are given independent normal prior distributions, α m [ i ] Z l ∼ i i d N ( 0 , τ l 2 ) and α m [ i ] x ∗ ∼ i i d N ( 0 , τ x ∗ 2 ) , and the variance components ( τ ( p 1 + 1 ) 2 , … , τ p 2 , τ x ∗ 2 ) are assigned a hyper‐prior distribution each using ∼ inv‐ χ 2 ( ν , σ 0 2 ) with a weakly informative or noninformative prior for ν and σ 0 . Using Bayesian simulations, a posterior draw of θ is then obtained by replacing Y ‾ j in ( ) with E ^ ( Y j | Z , X ∗ ) ( d ) , the d th draw from the posterior distribution of E ( Y | Z , X ∗ ) for units in cell j , with d = 1 , … , D , (4) θ ^ ( d ) = ∑ j = 1 J N j E ^ ( Y j | Z , X ∗ ) ( d ) ∑ j = 1 J N j The posterior mean or median of θ ^ ( d ) serves as the estimate for θ . The hierarchical structure in MRP partially pools estimates in poststrata formed by the auxiliary variables [ , , ], and improves the estimation of Y ‾ j in cells with sparse data by borrowing data information from other poststrata. It also improves efficiency in the estimation of population means in the presence of dispersed weights when models are well constructed . The model in ( ) can also be extended to include two‐way or higher order interactions between auxiliary variables . 2.2 Generalized Additive Model of Inclusion Propensity When the administrative records contain continuous auxiliary variables, discretizing the continuous variables may result in the loss of important information, especially when there are strong smooth relationships between the continuous auxiliary variables and the survey outcome. In the missing data literature, Little and An  showed that a penalized spline of propensity prediction model can improve mean estimation when imputation models are misspecified, where the logit‐transformed estimated response propensity is included in the model using a spline. We borrow this modeling idea. In addition to including X ∗ as covariates in the prediction model of the survey outcome, if the sample inclusion propensity is available, we can also include the logit‐transformed inclusion propensity as a covariate in the prediction model. With the sample inclusion propensity in conjunction with discretized variables, we aim to best utilize the continuous auxiliary information when access to the individual‐level continuous data is not feasible. Thus, we propose a two‐step strategy. In the first step, statistical agencies who have access to the confidential administrative records link survey sample with administrative records and estimate the inclusion propensity for both sample and nonsampled units using the continuous and discrete auxiliary variables measured in the population. They then create a modified population data set by replacing the continuous variables X with X ∗ and the estimated inclusion propensity. In the second step, the data users, who do not have access to administrative records, utilize this modified population data together with the sample data for survey inference. The sample inclusion propensity can be estimated using a propensity model [ , , ]. Let the sample inclusion indicator I be coded as 1 for the units in the sample and 0 for the rest of the units in the population. The estimated inclusion propensity, denoted by π ^ i for unit i in the population, can be obtained via fitting a logistic or probit regression of I on the discrete auxiliary variables Z and the original continuous variables X by linking survey sample with administrative records. When the number of auxiliary variables is small, parametric or semi‐parametric regression models can be used to estimate the propensity score. When there exist high‐dimensional auxiliary data, machine learning methods such as binary BART are recommended . The estimated sample inclusion propensity can be included in the modified population data together with the discrete auxiliary variables, and serve as a covariate in the prediction model for Y . Despite the advantages of MRP, the varying coefficients of α m [ i ] x ∗ in the MRP model in ( ) assume an exchangeable prior distribution. Because X ∗ is discretized from the continuous X , the ordering of the categories of X ∗ may matter. When there is a smooth relationship between X ∗ and Y , an exchangeable prior distribution is not sufficient. Therefore, we extend the MRP model in ( ) to allow a smooth relationship of Y with both X ∗ and logit ( π ^ ) , (5) g ( E ( Y i | Z i , X i ∗ ) ) = α 0 + ∑ k = 1 p 1 α k Z i k + ∑ l = p 1 + 1 p α m [ i ] Z l + s 1 ( x i ∗ ) + s 2 ( logit ( π ^ i ) ) where s 1 ( x i ∗ ) and s 2 ( logit ( π ^ i ) ) are smooth functions of x ∗ and logit ( π ^ i ) , respectively. These smooth functions allow flexible associations between the specified covariates and the survey outcome and thus protect against potential model misspecification. When the association between X and Y varies between different levels of a categorical variable, say Z 1 with H levels, model ( ) can be extended to reflect the interaction effect by allowing different smooth functions of X ∗ and logit ( π ^ ) in each level of Z 1 , (6) g ( E ( Y i | Z i , X i ∗ ) ) = α 0 + ∑ k = 1 p 1 α k Z i k + ∑ l = p 1 + 1 p α m [ i ] Z l + ∑ h = 1 H s 1 h ( x i ∗ ) + ∑ h = 1 H s 2 h ( logit ( π ^ i ) ) where s 1 h ( x i ∗ ) and s 2 h ( logit ( π ^ i ) ) are smooth functions for x ∗ and logit ( π ^ i ) in category h of Z 1 . The smooth functions can be modeled using spline or kernel functions. In this article, we use a smooth spline to model each of these associations. We execute the Bayesian generalized additive model (GAM) in ( ) and ( ) using the stan_gamm4() function in rstanarm , an R package that estimates models using RStan . RStan is an R interface to Stan for obtaining Bayesian inference using the No‐U‐Turn Sampler, a variant of Hamiltonian Monte Carlo . The stan_gamm4() function fits the specified GAM by adding priors on the hyperparameters of smooth splines, which is different from performing (restricted) maximum likelihood estimation with the lme4 package in R. Bayesian estimation provides better estimates for the uncertainty in the parameter estimates. We monitor the convergence of our parameter estimates using the convergence measure R ^ that suggests the chains mix well if close to 1. Models ( ) and ( ) are used to yield predictions for Y among the nonsampled units in the population [ , , ]. Our model‐based predictive estimator of θ is written as (7) θ ˜ m = N − 1 ∑ i ∈ s y i + ∑ i ∉ s ŷ i = N − 1 ∑ i = 1 N ŷ i + ∑ i ∈ s ( y i − ŷ i ) where y i denotes the observed Y in the sample, and ŷ i denotes the predicted Y based on the prediction model. The posterior distribution of θ is simulated by generating a large number of draws using the predictive estimator, with the median of these draws defining the generalized additive model of propensity (GAMP) estimator. 

Notation and Background We consider a target population consisting of N units with survey outcome Y . Let Z = ( Z 1 , … , Z p ) T denote p discrete auxiliary variables and X = ( X 1 , … , X q ) T be q continuous auxiliary variables, known for the population. For simplicity, we let q = 1 . The continuous variable X is usually not publicly available due to confidentiality concerns. Instead, before releasing to the public, X is discretized, denoted with X ∗ , using percentiles or meaningful cutoffs. The population can be partitioned into J disjoint and exhaustive cells or poststrata defined by the joint distribution of ( Z , X ∗ ) , with N j units in cell j , j = 1 … J and ∑ j = 1 J N j = N where N j > 0 . With Y ‾ j = ∑ i = 1 N j y i / N j , the population mean of Y within cell j , the overall population mean Y ‾ can be written as, (1) θ = Y ‾ = ∑ i = 1 N y i N = ∑ j = 1 J N j Y ‾ j ∑ j = 1 J N j Let s denote a probability survey sample of size n selected from the population with survey outcome values, y 1 , … , y n . We can then use ( Z , X ∗ ) to divide the sample into J cells; the corresponding sample size in cell j is n j with ∑ j = 1 J n j = n and the sample mean of Y in cell j is y ‾ j = ∑ i = 1 n j y i / n j . Assuming that y ‾ j is an unbiased estimate of Y ‾ j , the poststratification estimator of θ can be written as (2) θ ^ = ∑ j = 1 J N j y ‾ j ∑ j = 1 J N j = ∑ i = 1 n w i y i ∑ i = 1 n w i where w i = N j / n j for sample unit i in cell j , j = 1 , … , J . Alternatively, w i in formula ( ) can be created using the raking method when only population margins of ( Z , X ∗ ) are available. Raking weights are obtained via an iterative proportional fitting procedure that begins by adjusting design weights to the marginal distribution of the first auxiliary variable. These adjusted weights are then updated to conform to the marginal distribution of the second auxiliary variable. This process of updating the adjusted weights is carried out for each auxiliary variable and the first iteration ends when the weights are adjusted using the last auxiliary variable. Subsequent iterations are performed until the weights conform to the marginal distributions of all the auxiliary variables, that is, the algorithm converges . The Newton–Raphson method can also be used to obtain the same set of weights . When ( Z , X ∗ ) are related to survey outcome Y and the distributions of ( Z , X ∗ ) differ between sample and population, the poststratification and raking methods can effectively reduce bias in the estimation of θ , by assigning a weight w i for each sample unit. However, the weights can be highly variable, especially when p +q is large and n j is small for some j , j = 1 , … , J . Consequently, the weighted estimator can be unstable. Alternatively, model‐based approaches that predict the nonsampled units in the population based on a prediction model of the sample units can be used to improve efficiency in the θ estimation. The MRP model  has been widely used as a model‐based alternative to poststratification. MRP uses a hierarchical regression model to model the survey outcome of interest given the auxiliary variables, (3) g ( E ( Y i | Z i , X i ∗ ) ) = α 0 + ∑ k = 1 p 1 α k Z i k + ∑ l = p 1 + 1 p α m [ i ] Z l + α m [ i ] x ∗ where g ( · ) is a link function, such as using an identity link for continuous Y and a logit link for binary Y ; α 0 is the intercept; α k is the slope associated with the binary auxiliary variable Z k , k = 1 , … , p 1 and p 1 ≤ p ; and α m [ i ] Z l and α m [ i ] x ∗ correspond to varying coefficients associated with polytomous variables Z l , l = ( p 1 + 1 ) , … , p , and the discretized continuous variable X ∗ , respectively, where m [ i ] indicates the category of the polytomous variable to which the i th unit belongs. The varying coefficients α m [ i ] Z l and α m [ i ] x ∗ are given independent normal prior distributions, α m [ i ] Z l ∼ i i d N ( 0 , τ l 2 ) and α m [ i ] x ∗ ∼ i i d N ( 0 , τ x ∗ 2 ) , and the variance components ( τ ( p 1 + 1 ) 2 , … , τ p 2 , τ x ∗ 2 ) are assigned a hyper‐prior distribution each using ∼ inv‐ χ 2 ( ν , σ 0 2 ) with a weakly informative or noninformative prior for ν and σ 0 . Using Bayesian simulations, a posterior draw of θ is then obtained by replacing Y ‾ j in ( ) with E ^ ( Y j | Z , X ∗ ) ( d ) , the d th draw from the posterior distribution of E ( Y | Z , X ∗ ) for units in cell j , with d = 1 , … , D , (4) θ ^ ( d ) = ∑ j = 1 J N j E ^ ( Y j | Z , X ∗ ) ( d ) ∑ j = 1 J N j The posterior mean or median of θ ^ ( d ) serves as the estimate for θ . The hierarchical structure in MRP partially pools estimates in poststrata formed by the auxiliary variables [ , , ], and improves the estimation of Y ‾ j in cells with sparse data by borrowing data information from other poststrata. It also improves efficiency in the estimation of population means in the presence of dispersed weights when models are well constructed . The model in ( ) can also be extended to include two‐way or higher order interactions between auxiliary variables . 

Generalized Additive Model of Inclusion Propensity When the administrative records contain continuous auxiliary variables, discretizing the continuous variables may result in the loss of important information, especially when there are strong smooth relationships between the continuous auxiliary variables and the survey outcome. In the missing data literature, Little and An  showed that a penalized spline of propensity prediction model can improve mean estimation when imputation models are misspecified, where the logit‐transformed estimated response propensity is included in the model using a spline. We borrow this modeling idea. In addition to including X ∗ as covariates in the prediction model of the survey outcome, if the sample inclusion propensity is available, we can also include the logit‐transformed inclusion propensity as a covariate in the prediction model. With the sample inclusion propensity in conjunction with discretized variables, we aim to best utilize the continuous auxiliary information when access to the individual‐level continuous data is not feasible. Thus, we propose a two‐step strategy. In the first step, statistical agencies who have access to the confidential administrative records link survey sample with administrative records and estimate the inclusion propensity for both sample and nonsampled units using the continuous and discrete auxiliary variables measured in the population. They then create a modified population data set by replacing the continuous variables X with X ∗ and the estimated inclusion propensity. In the second step, the data users, who do not have access to administrative records, utilize this modified population data together with the sample data for survey inference. The sample inclusion propensity can be estimated using a propensity model [ , , ]. Let the sample inclusion indicator I be coded as 1 for the units in the sample and 0 for the rest of the units in the population. The estimated inclusion propensity, denoted by π ^ i for unit i in the population, can be obtained via fitting a logistic or probit regression of I on the discrete auxiliary variables Z and the original continuous variables X by linking survey sample with administrative records. When the number of auxiliary variables is small, parametric or semi‐parametric regression models can be used to estimate the propensity score. When there exist high‐dimensional auxiliary data, machine learning methods such as binary BART are recommended . The estimated sample inclusion propensity can be included in the modified population data together with the discrete auxiliary variables, and serve as a covariate in the prediction model for Y . Despite the advantages of MRP, the varying coefficients of α m [ i ] x ∗ in the MRP model in ( ) assume an exchangeable prior distribution. Because X ∗ is discretized from the continuous X , the ordering of the categories of X ∗ may matter. When there is a smooth relationship between X ∗ and Y , an exchangeable prior distribution is not sufficient. Therefore, we extend the MRP model in ( ) to allow a smooth relationship of Y with both X ∗ and logit ( π ^ ) , (5) g ( E ( Y i | Z i , X i ∗ ) ) = α 0 + ∑ k = 1 p 1 α k Z i k + ∑ l = p 1 + 1 p α m [ i ] Z l + s 1 ( x i ∗ ) + s 2 ( logit ( π ^ i ) ) where s 1 ( x i ∗ ) and s 2 ( logit ( π ^ i ) ) are smooth functions of x ∗ and logit ( π ^ i ) , respectively. These smooth functions allow flexible associations between the specified covariates and the survey outcome and thus protect against potential model misspecification. When the association between X and Y varies between different levels of a categorical variable, say Z 1 with H levels, model ( ) can be extended to reflect the interaction effect by allowing different smooth functions of X ∗ and logit ( π ^ ) in each level of Z 1 , (6) g ( E ( Y i | Z i , X i ∗ ) ) = α 0 + ∑ k = 1 p 1 α k Z i k + ∑ l = p 1 + 1 p α m [ i ] Z l + ∑ h = 1 H s 1 h ( x i ∗ ) + ∑ h = 1 H s 2 h ( logit ( π ^ i ) ) where s 1 h ( x i ∗ ) and s 2 h ( logit ( π ^ i ) ) are smooth functions for x ∗ and logit ( π ^ i ) in category h of Z 1 . The smooth functions can be modeled using spline or kernel functions. In this article, we use a smooth spline to model each of these associations. We execute the Bayesian generalized additive model (GAM) in ( ) and ( ) using the stan_gamm4() function in rstanarm , an R package that estimates models using RStan . RStan is an R interface to Stan for obtaining Bayesian inference using the No‐U‐Turn Sampler, a variant of Hamiltonian Monte Carlo . The stan_gamm4() function fits the specified GAM by adding priors on the hyperparameters of smooth splines, which is different from performing (restricted) maximum likelihood estimation with the lme4 package in R. Bayesian estimation provides better estimates for the uncertainty in the parameter estimates. We monitor the convergence of our parameter estimates using the convergence measure R ^ that suggests the chains mix well if close to 1. Models ( ) and ( ) are used to yield predictions for Y among the nonsampled units in the population [ , , ]. Our model‐based predictive estimator of θ is written as (7) θ ˜ m = N − 1 ∑ i ∈ s y i + ∑ i ∉ s ŷ i = N − 1 ∑ i = 1 N ŷ i + ∑ i ∈ s ( y i − ŷ i ) where y i denotes the observed Y in the sample, and ŷ i denotes the predicted Y based on the prediction model. The posterior distribution of θ is simulated by generating a large number of draws using the predictive estimator, with the median of these draws defining the generalized additive model of propensity (GAMP) estimator. 

Simulation Study 3.1 Design We simulate a population of size N = 3,000, and generate three independent binary variables Z = ( Z 1 , Z 2 , Z 3 ) with the marginal probabilities of ( 0 . 7 , 0 . 5 , 0 . 4 ) and one continuous variable X with a standard normal distribution. We consider two survey outcomes generated using additive nonlinear models. The first is a continuous outcome with Y 1 ∼ N ( 15 + 2 . 5 Z 1 − Z 2 + Z 3 − 2 X + 3 . 75 X 2 , 3 ) . The second is a binary outcome Y 2 generated via logit P r ( Y 2 = 1 ) = − 2 . 5 + 0 . 75 Z 1 − 2 . 5 Z 2 + 1 . 5 Z 3 − 0 . 25 X + 1 . 5 X 2 . Both models show a smooth association between X and survey outcomes. We then select a sample with approximate n = 600 cases from the population. We repeat this simulation process 500 times. We compare the performance of various estimators under three settings using absolute empirical bias, root mean squared error (RMSE), and average width and coverage rate of the 95% confidence or credible interval (CI). Credible intervals for Bayesian methods are calculated using the highest probability density method. Let θ ^ t be an estimate of the population mean θ t in the t t h simulation, t = 1 , … , 500 . The absolute empirical bias and RMSE are defined as follows, Absolute bias = | 1 500 ∑ t = 1 500 ( θ ^ t − θ t ) | , RMSE = 1 500 ∑ t = 1 500 ( θ ^ t − θ t ) 2 Estimators with smaller absolute bias, smaller RMSE, shorter 95% CI, and coverage rate closer to the nominal level are desired. 3.2 Setting One: X is Associated With Both Outcome Y and Inclusion I In this setting, we consider samples selected via an inclusion model, logit ( π ) = − 1 . 25 − Z 1 + 1 . 25 Z 2 − 0 . 75 Z 3 + 0 . 75 X − 0 . 10 X 2 , so that X is associated with both the survey outcomes Y 1 and Y 2 and the sample inclusion probability π . The propensity scores are estimated by fitting a binary BART model with predictors Z and X by linking the survey sample with the corresponding population data. BART is a sum‐of‐trees machine learning model . It is less sensitive to model misspecification and allows for nonlinear effects and multi‐way interactions between auxiliary variables and outcomes of interest. The estimated propensity scores π ^ are the posterior mean of the predictive probability of inclusion. We estimate population means for the overall population and for the subset where Z 1 = 1 . We compare the performance of various estimators, including an unweighted estimator using sample mean, eight weighted estimators, and four prediction model‐based estimators. The eight weighted estimators are IPW and eIPW : the inverse propensity weighted estimators. IPW uses true propensity scores π , while eIPW uses the estimated propensity scores π ^ from BART. PostStrat 3 and ( eIPW + PostStrat 3 ) : poststratification using the population joint distributions of Z and X 3 ∗ (the discretized X using population tertiles). ( eIPW + PostStrat 3 ) applies poststratification on the top of the base weights constructed using the inverse of the estimated propensity score from BART. Raking 5 and ( eIPW + Raking 5 ) : raking using the population margins of Z and X 5 ∗ (the discretized X using population quintiles). ( eIPW + Raking 5 ) applies raking on the top of the base weights constructed using the inverse of the estimated propensity score from BART. Raking 10 and ( eIPW + Raking 10 ) : raking using the population margins of Z and X 10 ∗ (the discretized X using population deciles) and raking on the top of the base weights. Poststratification using Z and X 5 ∗ or X 10 ∗ can lead to sparse or zero poststratification cells in the sample and thus is not considered here. Both the unweighted and weighted estimates are obtained using the survey package in R , in which the finite population correction is incorporated in variance estimation. We also consider four prediction models implemented using the rstanarm package in R : True Model: the true outcome model. MRP 10 : MRP model in ( ) using Z and X 10 ∗ . GAMP 10 and eGAMP 10 : GAMP model in ( ) using Z , X 10 ∗ , and logit propensity scores. GAMP 10 uses logit ( π ) while eGAMP 10 uses logit ( π ^ ) . Table presents the simulation results for this setting. The sample mean, which ignores the unequal probability of inclusion, performs poorly, exhibiting large bias and RMSE, with the 95% CI yielding a zero coverage rate. The eight weighted estimators demonstrate much smaller bias and RMSE compared with the sample mean. Among the eight weighted estimators, IPW yields the smallest bias but the largest RMSE. Using the estimated propensity score increases bias but reduces RMSE, which is expected due to the lower variation in the estimated propensity score compared with the true propensity score. Applying poststratification and raking on top of the base weights constructed with estimated propensity scores results in smaller bias and RMSE compared with poststratification and raking without considering the base weights for PostStrat 3 and raking 5 , but the improvement is small for raking 10 . Additionally, raking estimators using X 10 ∗ perform better than those using X 5 ∗ and yield the smallest RMSE among all weighted estimators. As expected, the predictive estimator using the true outcome model performs best, with the smallest bias, RMSE and average width, and coverage rate close to the nominal level. The MRP 10 performs similarly to raking 10 . The proposed GAMP 10 yields smaller bias but a wider 95% CI compared with MRP 10 for Y 1 and achieves both smaller bias and a shorter 95% CI for Y 2 . When replacing the true response propensity with the estimated one, the bias increases but the RMSE decreases and coverage rate of 95% CI improves for Y 1 while the results for Y 2 remain largely unchanged. The 95% CI coverage rates are significantly below the nominal level for all methods except the true model for Y 1 . This occurs because the continuous covariate X has a strong smooth relationship with both Y 1 and inclusion I , but discretizing X into quintiles or deciles across methods leads to information loss. The results of the subgroup analysis are presented in Supporting Information eTable 1. The conclusions are similar to the overall population analysis, except that the model‐based estimators now show greater efficiency gains compared with the weighted estimators. 3.3 Setting Two: X is Associated With Y but Not I In the second setting, we simulate one additional independent continuous auxiliary variable from a standard normal distribution, W ∼ N ( 0 , 1 ) , in the population. Then, samples are selected using inclusion model, logit ( π ) = − 1 . 25 − Z 1 + 1 . 25 Z 2 − 0 . 75 Z 3 + 0 . 75 W − 0 . 1 W 2 . In this setting, predictors Z , W , X are used to fit BART to estimate propensity scores π ^ . Note that W is associated with the sample inclusion probability π but not the survey outcomes Y 1 or Y 2 , and X is associated with the survey outcomes but not the sample inclusion. We compare the same thirteen estimators as setting one. All the weighted and predictive estimators, except for the true model estimator, use all available auxiliary information of ( Z , X ∗ , W ∗ ) . Specifically, PostStrat 3 uses ( Z , X 3 ∗ , W 3 ∗ ), raking 5 uses ( Z , X 5 ∗ , W 5 ∗ ), and raking 10 uses ( Z , X 10 ∗ , W 10 ∗ ), where W 3 ∗ , W 5 ∗ and W 10 ∗ are the discretized W using population tertiles, quintiles, and deciles, respectively. The MRP 10 uses ( Z , W 10 ∗ , X 10 ∗ ), whereas GAMP 10 and eGAMP 10 also include logit ( π ) and logit ( π ^ ) , respectively. Table presents the results for estimators on the overall population. The results for the subgroup analysis are shown in Supporting Information eTable 2. Similar to the first setting, the sample mean performs poorly. Among the eight weighted estimators, IPW yields the smallest bias but the largest RMSE. Using estimated response propensity, eIPW results in larger bias but smaller RMSE compared with IPW. The raking 5 shows a smaller bias and RMSE and a 95% CI coverage rate closer to the nominal level than PostStrat 5 . However, unlike setting one, raking 10 does not lead to significant improvements than raking 5 , with improvements only seen in RMSE for Y 1 . Applying poststratification and raking on the top of the base weights constructed using the inverse of propensity score does not lead to clear improvements compared with those ignoring the base weights, resulting in reduced RMSE but larger bias for Y 1 and reduced bias but larger RMSE for Y 2 . All the four prediction model‐based estimators have lower bias and reduced RMSE than the weighted estimators for both Y 1 and Y 2 . The inclusion of W in this setting increases the variation in the sample weights and thus the uncertainty in the weighted estimates, but the impact of W is relatively small in the prediction model‐based estimators. Because only W (and not X ) is related to π , the estimated inclusion propensity π ^ does not contain useful information about X . As such, the GAMP 10 and eGAMP 10 that incorporate the inclusion propensities does not offer improvements over MRP 10 . All model‐based estimators yield 95% CIs with coverage rate close to the nominal level. 3.4 Setting Three: Nonadditive Association With Both Y and I In the third simulation setting, we consider two survey outcomes generated from models with nonadditive association between X and Y . The first is a continuous outcome Y 1 ∼ N ( 15 + 2 . 5 Z 1 − Z 2 + Z 3 − 2 X + X 2 + Z 1 ∗ X − 2 . 5 Z 1 ∗ X 2 , 2 ) . The second is a binary outcome generated via logit P r ( Y 2 = 1 ) = − 1 . 75 + 0 . 75 Z 1 − 1 . 5 Z 2 + 1 . 5 Z 3 − 1 . 5 X + X 2 + Z 1 ∗ X − 2 . 5 Z 1 ∗ X 2 . Then, samples are selected using inclusion model logit ( π ) = − 0 . 9 − 0 . 5 Z 1 + 0 . 75 Z 2 − Z 3 + 0 . 5 X − 0 . 05 X 2 + 0 . 5 Z 1 ∗ X − 0 . 75 Z 1 ∗ X 2 , which also includes nonadditive association between X and π . We consider the similar weighted and predictive estimators as before. For raking, we use one‐variable margins for Z 2 and Z 3 and two‐variable margins for Z 1 and X 5 ∗ to model the interaction between Z 1 and X . The raking estimators using X 10 ∗ are not included due to the potential for sparse samples in the two‐variable margins for Z 1 and X 10 ∗ . For MRP 10 , we include interactions between Z 1 and X 10 ∗ . For GAMP 10 , the interaction between Z 1 and logit ( π ) is also included. Similarly, interaction between Z 1 and logit ( π ^ ) is specified for eGAMP 10 . To further assess whether the predictive estimators are sensitive to misspecification in the prediction models, we also compare eGAMP 10 to two alternative eGAMP models that either omit the interaction between Z 1 and X 10 ∗ or omit both interactions and compare MRP 10 to the alternative that ignores the interaction between Z 1 and X 10 ∗ for the continuous outcome Y 1 . Table presents the results for this setting. The results for the subgroup analysis are available in Supporting Information eTable 3. The results on the sample mean and weighted estimators are similar to the first setting. However, unlike the first setting, the predictive estimators show significant improvements over the weighted estimators, with much smaller bias, lower RMSE, and 95% CI coverage rates closer to the nominal level. Compared with MRP 10 , GAMP 10 and eGAMP 10 result in lower bias, improved efficiency, and closer to nominal level coverage, especially for the continuous outcome Y 1 . Table shows the sensitivity analysis for raking 5 , MRP 10 , and eGAMP 10 . Omitting the interaction between Z 1 and X ∗ results in much large bias and RMSE for all three estimators. However, the impact on the eGAMP 10 ( omitting interaction Z 1 × X 10 ∗ ) is smallest among the three estimators, with smaller bias and RMSE than the other estimators. Omitting both the Z 1 × X 10 ∗ and Z 1 × logit ( π ^ i ) interactions leads to poor performance of eGAMP 10 . 

Design We simulate a population of size N = 3,000, and generate three independent binary variables Z = ( Z 1 , Z 2 , Z 3 ) with the marginal probabilities of ( 0 . 7 , 0 . 5 , 0 . 4 ) and one continuous variable X with a standard normal distribution. We consider two survey outcomes generated using additive nonlinear models. The first is a continuous outcome with Y 1 ∼ N ( 15 + 2 . 5 Z 1 − Z 2 + Z 3 − 2 X + 3 . 75 X 2 , 3 ) . The second is a binary outcome Y 2 generated via logit P r ( Y 2 = 1 ) = − 2 . 5 + 0 . 75 Z 1 − 2 . 5 Z 2 + 1 . 5 Z 3 − 0 . 25 X + 1 . 5 X 2 . Both models show a smooth association between X and survey outcomes. We then select a sample with approximate n = 600 cases from the population. We repeat this simulation process 500 times. We compare the performance of various estimators under three settings using absolute empirical bias, root mean squared error (RMSE), and average width and coverage rate of the 95% confidence or credible interval (CI). Credible intervals for Bayesian methods are calculated using the highest probability density method. Let θ ^ t be an estimate of the population mean θ t in the t t h simulation, t = 1 , … , 500 . The absolute empirical bias and RMSE are defined as follows, Absolute bias = | 1 500 ∑ t = 1 500 ( θ ^ t − θ t ) | , RMSE = 1 500 ∑ t = 1 500 ( θ ^ t − θ t ) 2 Estimators with smaller absolute bias, smaller RMSE, shorter 95% CI, and coverage rate closer to the nominal level are desired. 

Setting One: X is Associated With Both Outcome Y and Inclusion I In this setting, we consider samples selected via an inclusion model, logit ( π ) = − 1 . 25 − Z 1 + 1 . 25 Z 2 − 0 . 75 Z 3 + 0 . 75 X − 0 . 10 X 2 , so that X is associated with both the survey outcomes Y 1 and Y 2 and the sample inclusion probability π . The propensity scores are estimated by fitting a binary BART model with predictors Z and X by linking the survey sample with the corresponding population data. BART is a sum‐of‐trees machine learning model . It is less sensitive to model misspecification and allows for nonlinear effects and multi‐way interactions between auxiliary variables and outcomes of interest. The estimated propensity scores π ^ are the posterior mean of the predictive probability of inclusion. We estimate population means for the overall population and for the subset where Z 1 = 1 . We compare the performance of various estimators, including an unweighted estimator using sample mean, eight weighted estimators, and four prediction model‐based estimators. The eight weighted estimators are IPW and eIPW : the inverse propensity weighted estimators. IPW uses true propensity scores π , while eIPW uses the estimated propensity scores π ^ from BART. PostStrat 3 and ( eIPW + PostStrat 3 ) : poststratification using the population joint distributions of Z and X 3 ∗ (the discretized X using population tertiles). ( eIPW + PostStrat 3 ) applies poststratification on the top of the base weights constructed using the inverse of the estimated propensity score from BART. Raking 5 and ( eIPW + Raking 5 ) : raking using the population margins of Z and X 5 ∗ (the discretized X using population quintiles). ( eIPW + Raking 5 ) applies raking on the top of the base weights constructed using the inverse of the estimated propensity score from BART. Raking 10 and ( eIPW + Raking 10 ) : raking using the population margins of Z and X 10 ∗ (the discretized X using population deciles) and raking on the top of the base weights. Poststratification using Z and X 5 ∗ or X 10 ∗ can lead to sparse or zero poststratification cells in the sample and thus is not considered here. Both the unweighted and weighted estimates are obtained using the survey package in R , in which the finite population correction is incorporated in variance estimation. We also consider four prediction models implemented using the rstanarm package in R : True Model: the true outcome model. MRP 10 : MRP model in ( ) using Z and X 10 ∗ . GAMP 10 and eGAMP 10 : GAMP model in ( ) using Z , X 10 ∗ , and logit propensity scores. GAMP 10 uses logit ( π ) while eGAMP 10 uses logit ( π ^ ) . Table presents the simulation results for this setting. The sample mean, which ignores the unequal probability of inclusion, performs poorly, exhibiting large bias and RMSE, with the 95% CI yielding a zero coverage rate. The eight weighted estimators demonstrate much smaller bias and RMSE compared with the sample mean. Among the eight weighted estimators, IPW yields the smallest bias but the largest RMSE. Using the estimated propensity score increases bias but reduces RMSE, which is expected due to the lower variation in the estimated propensity score compared with the true propensity score. Applying poststratification and raking on top of the base weights constructed with estimated propensity scores results in smaller bias and RMSE compared with poststratification and raking without considering the base weights for PostStrat 3 and raking 5 , but the improvement is small for raking 10 . Additionally, raking estimators using X 10 ∗ perform better than those using X 5 ∗ and yield the smallest RMSE among all weighted estimators. As expected, the predictive estimator using the true outcome model performs best, with the smallest bias, RMSE and average width, and coverage rate close to the nominal level. The MRP 10 performs similarly to raking 10 . The proposed GAMP 10 yields smaller bias but a wider 95% CI compared with MRP 10 for Y 1 and achieves both smaller bias and a shorter 95% CI for Y 2 . When replacing the true response propensity with the estimated one, the bias increases but the RMSE decreases and coverage rate of 95% CI improves for Y 1 while the results for Y 2 remain largely unchanged. The 95% CI coverage rates are significantly below the nominal level for all methods except the true model for Y 1 . This occurs because the continuous covariate X has a strong smooth relationship with both Y 1 and inclusion I , but discretizing X into quintiles or deciles across methods leads to information loss. The results of the subgroup analysis are presented in Supporting Information eTable 1. The conclusions are similar to the overall population analysis, except that the model‐based estimators now show greater efficiency gains compared with the weighted estimators. 

Setting Two: X is Associated With Y but Not I In the second setting, we simulate one additional independent continuous auxiliary variable from a standard normal distribution, W ∼ N ( 0 , 1 ) , in the population. Then, samples are selected using inclusion model, logit ( π ) = − 1 . 25 − Z 1 + 1 . 25 Z 2 − 0 . 75 Z 3 + 0 . 75 W − 0 . 1 W 2 . In this setting, predictors Z , W , X are used to fit BART to estimate propensity scores π ^ . Note that W is associated with the sample inclusion probability π but not the survey outcomes Y 1 or Y 2 , and X is associated with the survey outcomes but not the sample inclusion. We compare the same thirteen estimators as setting one. All the weighted and predictive estimators, except for the true model estimator, use all available auxiliary information of ( Z , X ∗ , W ∗ ) . Specifically, PostStrat 3 uses ( Z , X 3 ∗ , W 3 ∗ ), raking 5 uses ( Z , X 5 ∗ , W 5 ∗ ), and raking 10 uses ( Z , X 10 ∗ , W 10 ∗ ), where W 3 ∗ , W 5 ∗ and W 10 ∗ are the discretized W using population tertiles, quintiles, and deciles, respectively. The MRP 10 uses ( Z , W 10 ∗ , X 10 ∗ ), whereas GAMP 10 and eGAMP 10 also include logit ( π ) and logit ( π ^ ) , respectively. Table presents the results for estimators on the overall population. The results for the subgroup analysis are shown in Supporting Information eTable 2. Similar to the first setting, the sample mean performs poorly. Among the eight weighted estimators, IPW yields the smallest bias but the largest RMSE. Using estimated response propensity, eIPW results in larger bias but smaller RMSE compared with IPW. The raking 5 shows a smaller bias and RMSE and a 95% CI coverage rate closer to the nominal level than PostStrat 5 . However, unlike setting one, raking 10 does not lead to significant improvements than raking 5 , with improvements only seen in RMSE for Y 1 . Applying poststratification and raking on the top of the base weights constructed using the inverse of propensity score does not lead to clear improvements compared with those ignoring the base weights, resulting in reduced RMSE but larger bias for Y 1 and reduced bias but larger RMSE for Y 2 . All the four prediction model‐based estimators have lower bias and reduced RMSE than the weighted estimators for both Y 1 and Y 2 . The inclusion of W in this setting increases the variation in the sample weights and thus the uncertainty in the weighted estimates, but the impact of W is relatively small in the prediction model‐based estimators. Because only W (and not X ) is related to π , the estimated inclusion propensity π ^ does not contain useful information about X . As such, the GAMP 10 and eGAMP 10 that incorporate the inclusion propensities does not offer improvements over MRP 10 . All model‐based estimators yield 95% CIs with coverage rate close to the nominal level. 

Setting Three: Nonadditive Association With Both Y and I In the third simulation setting, we consider two survey outcomes generated from models with nonadditive association between X and Y . The first is a continuous outcome Y 1 ∼ N ( 15 + 2 . 5 Z 1 − Z 2 + Z 3 − 2 X + X 2 + Z 1 ∗ X − 2 . 5 Z 1 ∗ X 2 , 2 ) . The second is a binary outcome generated via logit P r ( Y 2 = 1 ) = − 1 . 75 + 0 . 75 Z 1 − 1 . 5 Z 2 + 1 . 5 Z 3 − 1 . 5 X + X 2 + Z 1 ∗ X − 2 . 5 Z 1 ∗ X 2 . Then, samples are selected using inclusion model logit ( π ) = − 0 . 9 − 0 . 5 Z 1 + 0 . 75 Z 2 − Z 3 + 0 . 5 X − 0 . 05 X 2 + 0 . 5 Z 1 ∗ X − 0 . 75 Z 1 ∗ X 2 , which also includes nonadditive association between X and π . We consider the similar weighted and predictive estimators as before. For raking, we use one‐variable margins for Z 2 and Z 3 and two‐variable margins for Z 1 and X 5 ∗ to model the interaction between Z 1 and X . The raking estimators using X 10 ∗ are not included due to the potential for sparse samples in the two‐variable margins for Z 1 and X 10 ∗ . For MRP 10 , we include interactions between Z 1 and X 10 ∗ . For GAMP 10 , the interaction between Z 1 and logit ( π ) is also included. Similarly, interaction between Z 1 and logit ( π ^ ) is specified for eGAMP 10 . To further assess whether the predictive estimators are sensitive to misspecification in the prediction models, we also compare eGAMP 10 to two alternative eGAMP models that either omit the interaction between Z 1 and X 10 ∗ or omit both interactions and compare MRP 10 to the alternative that ignores the interaction between Z 1 and X 10 ∗ for the continuous outcome Y 1 . Table presents the results for this setting. The results for the subgroup analysis are available in Supporting Information eTable 3. The results on the sample mean and weighted estimators are similar to the first setting. However, unlike the first setting, the predictive estimators show significant improvements over the weighted estimators, with much smaller bias, lower RMSE, and 95% CI coverage rates closer to the nominal level. Compared with MRP 10 , GAMP 10 and eGAMP 10 result in lower bias, improved efficiency, and closer to nominal level coverage, especially for the continuous outcome Y 1 . Table shows the sensitivity analysis for raking 5 , MRP 10 , and eGAMP 10 . Omitting the interaction between Z 1 and X ∗ results in much large bias and RMSE for all three estimators. However, the impact on the eGAMP 10 ( omitting interaction Z 1 × X 10 ∗ ) is smallest among the three estimators, with smaller bias and RMSE than the other estimators. Omitting both the Z 1 × X 10 ∗ and Z 1 × logit ( π ^ i ) interactions leads to poor performance of eGAMP 10 . 

Application to the Ohio Army National Guard Study The Ohio Army National Guard Mental Health Initiative (OHARNG‐MHI) study provides information about the prevalence and risk factors of mental health‐related outcomes among National Guard service members with the aim of identifying areas of intervention that can be modified during deployment to improve the psychological well‐being of soldiers . We apply the methods described in this paper to the OHARNG‐MHI study to estimate the percentage of lifetime alcohol abuse among all the service members as a data illustration. The target population for the OHARNG‐MHI study included all active members of the OHARNG between June 2008 and February 2009. All members with address information listed with the Guard were notified of the study via mailed letter and opt‐out card. Although some members chose to opt out of the study, others refused participation when contacted or were not contacted before the cohort closed . Furthermore, service members with no or incorrect telephone numbers could not be contacted to complete the 60 min structured computer‐assisted telephone interview . As such, the statistical analysis of the OHARNG‐MHI study data is complicated by nonsampling errors due to survey nonresponse and sampling frame undercoverage. Information on age group (17–24 years, 25–34 years, 35–44 years, 45 years or older), gender (male, female), race (White, Black, Other), rank (enlisted, officer), marital status (single, married, other), and number of years in service for the target population ( N = 10,994) was available in OHARNG administrative files. A total of n = 2,600 service members completed the survey. Table shows the distribution of the auxiliary information in the population and survey sample. It indicates that older, other race, officer rank, and married service members were overrepresented in the sample, whereas younger, Black race, enlisted rank, and single service members were underrepresented. The outcome of interest, presence of lifetime alcohol abuse, was only measured among survey participants. Figure shows the association between number of years in service and the logit‐transformed proportion of lifetime alcohol abuse in the sample. Here, the proportion of lifetime alcohol abuse is computed as the number of participants with lifetime alcohol abuse divided by the total number of participants in each year of service. The plot of the logit‐transformed proportion by years of service is overlaid with the fitted loess curve. The proportion increases with the number of years in service. Figure is created similarly and shows a smooth association between the logit‐transformed sample inclusion propensity and years in service. The logit‐transformed sample inclusion probability has a linear association with the number of years in service, increasing from 0 to about 20 years, and then plateaus. The data analysis includes three main steps. First, statisticians who have access to the confidential OHARNG administrative files, which contain the individual‐level continuous years of service variable, estimate the sample inclusion propensities. For this analysis, a binary BART machine learning model is utilized, incorporating all available auxiliary variables from the administrative files as covariates, including both discrete and continuous variables listed in Table . In the second step, a modified population data set is created for sharing with data users. The continuous years of service variable in the administrative files is replaced with its discretized version, categorized using deciles of the population values, along with the estimated inclusion propensities, π ^ . Finally, data users who do not have access to the confidential administrative files use the survey sample and the modified population data to conduct survey inference. The proportion of lifetime alcohol abuse in the OHARNG population is estimated using several methods: unweighted sample mean, eIPW, raking, eIPW + raking, MRP, and eGAMP. For raking, all discrete covariates from Table and the discretized version of years in service using deciles are used for weighting. For eGAMP, a GAM is fitted for the binary lifetime alcohol abuse variable, including all discrete covariates from Table , a spline of the deciles of the years of service variable, and a spline of the estimated inclusion propensities. For comparison, we also provide the GAM(x) estimator by fitting a GAM model regressing on a spline of the original continuous years in service variable and the other discrete auxiliary variables. Figure shows estimates and 95% CIs for the proportion of lifetime alcohol abuse. The unweighted sample mean estimates a higher proportion of lifetime alcohol abuse compared with the weighted and predictive estimators. This is expected, given that service members with more years of service were more likely to be included in the sample and also had a higher proportion of lifetime alcohol abuse (Figure ). All the weighted and predictive estimators yield similar point estimates for the proportion of lifetime alcohol abuse. However, eGAMP produces a wider 95% CI than the others. Due to the limited variation in the estimated inclusion propensities among the sample units in this application, we do not observe the typically wide confidence intervals associated with the weighted estimators. 

Discussion We consider finite population inference from a nonrepresentative sample where a number of auxiliary variables, both continuous and discrete, are measured in both the sample and the population via administrative records. This auxiliary information can be used to improve survey inference of population quantities through weighting or predictive models. However, individual‐level continuous administrative records are often inaccessible due to confidentiality concerns. The common practice is to discretize the continuous auxiliary variables using percentiles or other meaningful cutoffs, which can lead to a loss of information. Motivated by this challenge, we develop a method for predictive survey inference that makes the best use of continuous auxiliary information in the administrative records without requiring the release of individual‐level continuous variables for the entire population. We propose a two‐step strategy. In the first step, statistical agencies with access to confidential population data estimate probabilities of inclusion in the sample for all population units using both continuous and discrete auxiliary variables. They then create a modified population data for data users by replacing the continuous variables with their discretized versions, along with the continuous estimated inclusion propensities. In the second step, data users, who do not have access to the protected population data, utilize the modified population dataset and the survey sample for statistical inference of population quantities using the proposed GAMP model. The GAMP model extends MRP by using a hierarchical structure of discrete variables and spline functions of continuous variables, incorporating discretized continuous variables and estimated inclusion propensities. This approach is inspired by Little and An , who included spline functions of the logit‐transformed response propensity score in the prediction models, which yields robust estimates of sample means even when the model is misspecified. Similarly, including the logit‐transformed estimated propensities in GAMP's predictive models can improve survey inference of population quantities when discretized continuous variables do not fully capture the associations between the original continuous variables and the survey outcome of interest. We provide an R package, AuxSurvey , for conducting analyses using our proposed GAMP estimator. We assess the performance of the proposed GAMP method and compare it to existing weighting and prediction model‐based approaches through a simulation study. Our simulations confirm established findings on the importance of appropriate statistical analyses to adjust for nonsampling errors such as nonresponse and undercoverage [ , , ]. Both weighting and model‐based approaches yield more accurate population estimates than unadjusted estimates across all simulation settings. Importantly, our proposed GAMP estimator outperforms the weighting approaches using IPW, poststratification, or raking, providing less biased and more efficient estimates with a 95% CI coverage rate closer to the nominal level. GAMP also outperforms MRP when continuous auxiliary variables have more complicated associations with both the outcome of interest and sample inclusion. GAMP has a smaller bias and improved efficiency compared with modeling the discretized continuous variables alone in MRP (e.g., simulation setting three). When interactions exist between continuous variables and a categorical variable in the outcome model, omitting these interaction effects can lead to large bias and RMSE in raking and MRP estimators. In contrast, the bias and RMSE of GAMP are relatively small as long as the interactions between categorical variables and the logit‐transformed propensity scores are included, even if the interactions between categorical variables and discretized continuous variables are omitted. In cases where continuous auxiliary variables are not associated with the outcome of interest but are associated with sample inclusion (e.g., simulation setting two), the GAMP model, including a spline of the estimated inclusion propensities, does not offer additional benefits and incurs a slight increase in bias and RMSE compared with MRP. Nevertheless, the GAMP estimator still outperforms the weighted estimators in this setting. In practice, the number of auxiliary variables available for the entire finite population is often limited, so flexible parametric models are usually preferred. In this paper, we consider a Bayesian GAM implemented using the stan_gamm4() function in rstanarm . Using the hierarchical structure similar to MRP, our GAM model partially pools estimates in multiple categories of discrete auxiliary covariates. Including spline functions of the discretized continuous variables and the estimated inclusion propensities, the model is flexible enough to catch their potential nonlinear associations with the survey outcome of interest. When there are interactions between the smooth splines and discrete auxiliary variables, the model can also be easily modified to reflect these interaction effects using stratified splines. Furthermore, the Bayesian framework is straightforward for quantification of uncertainty. In the settings where there exists a high dimension of auxiliary variables, machine learning techniques, such as BART, can be used instead. It is crucial to balance data utility and confidentiality protection. The two‐step strategy we propose aims to ensure that continuous auxiliary data in administrative records remain useful for survey inference while protecting confidentiality. Although we focus on continuous auxiliary variables, the proposed method can also be applied to settings with sensitive categorical covariates, such as fine‐level geographic identifiers. In such cases, the fine‐level geographic identifier is used to estimate the inclusion propensity, which is then released along with coarse‐level geographic information for predictive survey inference. To avoid assigning the exact same estimated propensity score to units with identical categorical predictors, a random draw from the posterior distribution of the propensity score—rather than the posterior mean—can be obtained from the binary BART model for inclusion propensity. The proposed GAMP method has several limitations. First, like other predictive estimators, its effectiveness depends on the inclusion of auxiliary variables that are predictive of the outcomes. Additionally, these auxiliary variables in the sample should have common support to those in the target population. If important predictors of outcomes are not available in the administrative records, or if the ranges of the auxiliary variables in the sample are narrower than those in the population, the model predictions may not perform well. Second, the sample inclusion propensity often requires estimation. Our simulation shows that the performance of GAMP using the true versus BART‐estimated propensity scores is comparable, although the eGAMP estimator (using the BART‐estimated inclusion propensity) results in larger bias but smaller RMSE compared with GAMP using the true inclusion propensity. In the eGAMP estimator, we ignore the uncertainty associated with π estimation. One way to address the uncertainty from estimating π is obtaining multiple sets of estimated propensity scores from their posterior distributions and repeating the predictive models multiple times with these sets of estimated π values. Third, estimating inclusion propensity often involves linking the survey sample with administrative records, which may be subject to record linkage errors and need further correction. In conclusion, our study advocates for a prediction model‐based approach that leverages continuous auxiliary information in administrative records to improve survey inference while controlling disclosure risk, thus eliminating the need to release individual‐level continuous auxiliary variables. When survey inference is conducted by those with access to confidential records, a predictive model using splines on the continuous variables can be directly applied. However, for data users without such access, our two‐step strategy offers an effective solution. Although the use of predictive inferences, Bayesian hierarchical models, and splines for propensity scores is not new, our major contribution lies in combining these components to leverage record‐level continuous variables for improving survey inference in nonrepresentative probability samples or nonprobability samples while maintaining confidentiality. This paper focuses on model‐based predictive inference that includes the propensity score as a covariate. An alternative approach that combines weighting and prediction is model‐assisted estimation . This research also opens up several exciting directions for future investigation. For example, modeling two separate propensities—one for undercoverage (noncontact) and one for response (responded out of those contacted)—and including them as smoothed terms in the prediction model could be explored. Another potential extension is to model the unknown population distribution of auxiliary variables and incorporate the uncertainty in estimating this distribution for survey inference . 

Sharifa Z. Williams and Qixuan Chen conceived and designed the research and wrote the paper with contributions from Yajuan Si, Jungang Zou, and Sandro Galea. Sharifa Z. Williams, Jungang Zou, and Yutao Liu conducted the simulation study. Sharifa Z. Williams performed the real data analysis. Jungang Zou wrote the R package. Sharifa Z. Williams, Qixuan Chen, and Yajuan Si interpreted the results. Sandro Galea contributed the application data from the Ohio Army National Guard Mental Health Initiative study. 

Yutao Liu is now an employee of AstraZeneca and may or may not own stock options. 

The authors declare no conflicts of interest. 

Data S1. Supporting Information. 

