Combining Non-probability and Probability Survey Samples Through Mass Imputation
Jae Kwang Kim Seho Park Yilin Chen Changbao Wu




#1 1



1
  
  

  Combining Non-probability and Probability Survey Samples Through Mass Imputation
    Jae Kwang Kim^1,
  	 Seho Park^2, Yilin Chen^3, and Changbao Wu^3

    ^1 Department of Statistics, Iowa State University,

     ^2 Department of Biostatistics, Indiana University School of Medicine, 

    ^3 Department of Statistics and Actuarial Science,
University of Waterloo
    
===========================================================================================================================================================================================================================================================================================

 

0
  

    Combining Non-probability and Probability Survey Samples Through Mass Imputation


   




Summary. Analysis of non-probability survey samples requires auxiliary information at the population level. Such information may also be obtained from an existing probability survey sample from the same finite population. 
Mass imputation has been used in practice for combining non-probability and probability survey samples and making inferences on the parameters of interest  using the information collected only in the non-probability sample for the study variables.   

Under the assumption that the conditional mean function from the non-probability sample can be transported to the probability sample, we establish the consistency of the mass imputation estimator  and derive its asymptotic variance formula. Variance estimators are developed using either linearization or bootstrap. Finite sample performances of the mass imputation estimator are investigated through simulation studies. We also address important practical issues of the method through the analysis of a real world non-probability survey sample collected by the Pew Research Centre.


Keywords:  Auxiliary variables; Bootstrap variance estimator; Data integration; Ignorable sample selection; Model transportability; Selection bias



1 




§ INTRODUCTION



Probability sampling is a classical tool for obtaining a representative sample from a target population. Because the first-order inclusion probabilities are known under the given survey design,  probability sampling can provide design unbiased estimators and valid statistical inferences for finite population parameters. Although probability samples have been a tremendously successful story over the past 60 years, the approach has encountered many challenges in recent years, including high cost, low response rate, and inability to provide up-to-date information on variables of specific studies.
On the other hand, non-probability samples, such as those collected from web panels, have become increasingly popular due to the efficient recruitment process, quick responses and low maintenance expenses. See, for instance, <cit.>, for examples of web based non-probability survey samples. 

Statistical analysis of non-probability survey samples, however, faces many challenges as documented by <cit.>. The first major challenge is that the selection/inclusion mechanism for non-probability samples is typically unknown. Treating non-probability samples as if they are a simple random sample often leads to biased results. The second major challenge is that the estimation of the propensity scores (i.e., the probabilities of participating in the non-probability sample) under an assumed model requires auxiliary information at the population level. A popular framework used in recent years is to assume that the required auxiliary information on the target population is available from an existing probability survey sample. This framework was first used by <cit.> and followed by a number of other authors, including <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>, among others. 


Generally speaking, there are three possible approaches to analyzing non-probability survey samples: (i) inverse probability weighting (IPW) using estimated propensity scores; (ii) model-based prediction using an outcome regression model; and (iii) doubly robust procedures involving both the propensity scores and the outcome regression model. A central point to all approaches is the assumption that 
the selection/inclusion mechanism for non-probability survey samples is ignorable. It is essentially the missing at random (MAR) assumption of <cit.>. It is well known in the literature on missing data analysis that the MAR assumption cannot be formally tested using data from the sample itself. Practical applications of the propensity score based methods require a careful examination of the auxiliary variables included in the sample to see whether participation in the non-probability sample could be characterized by those variables. Calibration weighting methods for non-probability samples are also discussed in the literature. See, for instance, <cit.> and  <cit.>. The proposed methodologies, however, are closely related to the propensity scores based methods and require the same ignorability assumption on the selection mechanism <cit.>. 

Mass imputation, on the other hand,  is a model-based prediction method that includes both the non-probability and the probability samples. 

The existing probability sample is assumed to have no measurement on the study variable of interest and can be viewed as having 100% missing values for the study variable.
The non-probability sample contains  values observed on both the study variable and the auxiliary variables. 
The non-probability sample is then used as  training data to develop a prediction model that will be used to provide imputations  for the probability sample.
The key assumption for the validity of the mass imputation method is that the prediction model built based on the non-probability sample can be transported to the probability sample for imputation. See Section <ref> for detailed discussions. Mass imputation has also been developed in the context of two-phase sampling (), but it is not fully investigated in the context of survey integration for combining the non-probability sample with a probability sample. One notable exception is the sample matching method of <cit.>, but he did not provide any theoretical justifications for the proposed method. Furthermore, the sample matching method of  <cit.> is based on the nearest neighbor imputation method, which suffers from the curse of dimensionality when the number of auxiliary variables is large.
<cit.>
discussed composite estimation when one of the surveys is 
mass-imputed.
<cit.> discussed practical issues in sample
matching for mass imputation.

In this paper, we aim to fill the important research gap in survey sampling on mass imputation for analyzing non-probability survey samples. From a theoretical point of view, if the training data for building the imputation model were a probability sample, then the theory of <cit.> could be  directly  applicable. Under the assumptions that the support of auxiliary variables from the non-probability sample is the same as the sample space and the conditional mean function from the non-probability sample can be transported to the probability sample, we show that the method of <cit.> can be applied even when  non-probability samples are used as the training data. We develop rigorous asymptotic theory for the mass imputation estimator along with a linearization variance estimator and a proposed bootstrap method for variance estimation. Our investigation on practical aspects of the method was motivated by the analysis of a real non-probability survey sample collected in 2015 by the Pew Research Centre (PRC) from the United States of America. It turned out that there existed at least two probability survey samples taken in the same year from the same finite population with some common auxiliary variables included in both samples. In addition to demonstrations of the mass imputation method using the PRC non-probability sample, we also address major practical questions through the example, which includes (i) the assessment of auxiliary variables with regard to the ignorability assumption and the imputation model; (ii) the impact of relative sample sizes between the non-probability and the probability samples; and (iii) factors to be considered in choosing a probability sample when multiple samples are available. 


The paper is organized as follows. 
The basic setting is described in Section <ref> and the mass imputation is presented in Section <ref>. Main theoretic results on consistency and the asymptotic variance formula are presented in Section <ref>. A practically useful bootstrap variance estimator is proposed in Section <ref>. Results from a limited simulation study on the finite sample performances of the mass imputation estimator are reported in Section <ref>.
Important practical aspects of the method are discussed in Section <ref> through an application to the PRC non-probability survey sample. Some concluding remarks are given in Section <ref>. Proofs are presented in the online Supplementary Material document.








§ PROBLEM SETUP
 


Let (_i,y_i) be the measurements of the auxiliary variables  and the study variable Y associated with unit i, i=1,⋯,N, where N is the population size. Let θ_N = N^-1∑_i=1^Ny_i be the finite population mean of the study variable, which is the parameter of interest. 
Let A be a probability sample with observations on the auxiliary variables ; let B be the non-probability sample with information on both the study variable Y and the auxiliary variables . 
Let n_A=|A| and n_B=|B| be the respective sample sizes. Table <ref> presents the general setup of the two sample structure for data integration. The non-probability sample B  is not representative of the target population due to the unknown sample selection/inclusion mechanism. 



The two sample datasets can also be represented by {(_i,w_i), i ∈ A} and {(_i,y_i),i∈ B}, where w_i=π_i^-1 are the survey weights for the probability sample A and π_i = P(i∈ A) are the first-order inclusion probabilities.
If the study variable Y were observed for the probability sample A, the Horvitz-Thompson estimator of θ_N would be given by θ̂_HT = N^-1∑_i∈ Aw_iy_i. The mass imputation estimator to be introduced in the next section replaces the unobserved y_i by an imputed ŷ_i using data from both samples A and B. 


Our discussions on important practical issues related to the mass imputation estimator are motivated by the analysis of the survey sample collected by the Pew Research Centre in 2015. 

The dataset, denoted by PRC, contains  56 variables which aim to reveal the relations between people and their 
communities. There is a total of 9,301 cases in the PRC dataset which  are provided by eight different vendors with unknown sampling and data collection strategies. We treat the PRC dataset as a non-probability survey sample with the sample size n_B=9,301.
Four variables of the PRC dataset are of particular interest to our analysis, which are pertinent to the research questions on “people and their communities”, and are treated as study variables. They are listed at the bottom of Table <ref>, among them three are binary variables: Talk with neighbors frequently (Y_1), Participated in school groups (Y_2), Participated in service organizations (Y_3), and one is treated as a continuous  variable: Days had at least one drink last month (Y_4).






The general framework for mass imputation requires the availability of an existing probability sample taken from the same target population and containing a set of auxiliary variables which are common to the non-probability sample. It turns out that there are (at least) two such probability survey samples available, both were taken in 2015 from the general US population. This raises an important practical question on how mass imputation can be implemented under such scenarios, which will be further discussed in Sections <ref> and <ref>. 


The first probability sample is the volunteer supplement survey data from the Current Population Survey (CPS), which is one of the most  recognized surveys in the United States. The CPS dataset contains n_A = 80,075 cases with measurements on volunteerism, which is highly relevant to the study variables considered in the PRC dataset. The second probability  sample is the Behavioral Risk Factor Surveillance System (BRFSS) survey data. It  is  designed to measure behavioral risk factors for US residents and has a large sample size n_A = 441,456.  
Neither of the two probability  samples contains  measurements of the study variables, but both share a rich set of common survey items with the PRC dataset as shown in Table <ref>.

As a preliminary analysis, we examine marginal distributions of the variables contained in the three datasets. 
Table <ref> contains the estimated population mean using each of the three datasets. For the PRC dataset, the sampling strategy is unknown and no survey weights are  available, so the estimates presented are the unadjusted simple sample means. For the BRFSS and the CPS datasets where survey weights are available, survey weighed estimates are computed.
The “NA” in the table indicates that the variable is not available from the dataset. It can be seen that there are noticeable differences between the naive sample means from the PRC sample and the survey weighted estimates from the two probability samples for variables such as Origin (Hispanic/Latino), Education (High school or less), Household (with children), Health (Smoking) and Volunteer works. It is clear evidence that the PRC dataset is not a representative sample for the US population. 







§ THE MASS IMPUTATION ESTIMATOR



We now formally define the mass imputation estimator of the population mean θ_N under the two-sample setup presented in Table <ref>.
Let δ be the indicator variable for the unit being included in the non-probability sample B. Note that δ = 1 for all i∈ B and δ = 0 for all i ∉ B. This is fundamentally different from the standard setting for missing data analysis where the response indicator variable equals to 1 for respondents and 0 for nonrespondents, all part of the sample. It is also similar to missing data problems if we treat the entire finite population as a sample, which intuitively justifies the need for auxiliary information at the population level in analyzing non-probability survey samples.  


We first assume that each unit in the population has a non-zero probability to be included in the sample B, i.e.,

    P ( δ =1 |=) >0

for all  in the support of . 
Condition (<ref>) is often referred as as the Positivity  Assumption, which means that for any possible value , there is a positive probability for this type of units to be selected for sample B. It implies that the support of  in sample B coincides with the sample space of  in the population. The condition can easily be verified if all components of  are discrete. 

If certain components of  are continuous, it is still plausible
to check the condition based on the distributions 
of  between the weighted distribution in the probability sample A and the empirical distribution in the non-probability sample B. If the two sample supports overlap to each other, then (<ref>) is satisfied. <cit.> proposed a systematic way of  subsampling to handle the cases with limited overlaps in the context of causal inference.   


The most crucial part of the mass imputation approach is the prediction model f(y|), which can be estimated by using the observed (Y,) from the non-probability sample B if 

    f( y|, δ=1) = f ( y |)  .

The prediction model f(y|) can then be used to create mass-imputed values of y using observed  for all the units in the probability sample A. 

Equation (<ref>) can be termed as the transportability condition in the sense that the imputation model built from sample B is transportable to sample A. 
A sufficient condition for (<ref>) to hold is the ignorability condition for the non-probability sample B, which is similar to the missing at random (MAR) assumption  <cit.> widely used in the literature on missing data analysis:

    P( δ = 1 |, Y) = P( δ =1 |)  .

In other words, transportability is a weaker condition than ignorability. Unfortunately, the ignorability condition cannot be formally tested using the sample itself. 









Under assumptions (<ref>)-(<ref>), it is possible to
consider a mass imputation estimator based on nearest neighbor imputation as suggested by <cit.>. Nearest neighbor imputation is a nonparametric method that does not require any parametric model assumptions. While nonparametric imputation methods can provide robust estimation, it suffers from curse of dimensionality, and the asymptotic bias of the nearest neighbor imputation is not negligible if  the dimension of  is greater than one <cit.>.
In this paper, we consider a semi-parametric model for sample B with the first conditional moment specified as

    E(Y | = )=m(;_0)

for some unknown p× 1 vector _0 and a known function m(·;·).  
The specification of m(;) typically follows the mean function for generalized linear models. If Y is continuous, we can use a linear regression model and let m(;) = '; if Y is binary, we can use a logistic regression model and let m(;) =  {1+exp(-') }^-1; if Y represents count data, we can use a log-linear model and let m(;) = exp('). 
Under the assumed transportability condition, the model parameters  can be estimated by the non-probability sample, B. 

We assume that  is the unique solution to the estimating equations 

    Û ( ) = 1/n_B∑_i ∈ B{ y_i-m(_i;) }(_i;)=0

for some p-dimensional vector of functions (;). 

The estimating equations given in (<ref>) may be constructed as the quasi-score equations from generalized linear models.  

Note that E{Û ( _0) }=0 under the assumption (<ref>) and the model (<ref>).
Under suitable moment conditions on (Y, ) and smoothness conditions on m(;) with respect to  similar to those discussed in Section 3.2 of <cit.>, it can be shown that the solution to (<ref>), denoted as , is consistent for _0 under the ignorability assumption.


Once   is obtained from the non-probability sample B  where both Y and  are observed, we can obtain predicted values ŷ_i = m( 𝐱_i ; ) for all i ∈ A, which is the so-called “mass imputation”.
The mass imputation estimator for the finite population mean θ_N=N^-1∑_i=1^N y_i is computed as

    θ̂_I=1/N∑_i ∈ Aw_iŷ_i

using the survey weights w_i for the probability sample A. If the population size N is unknown, it can be estimated by N̂ = ∑_i∈ A w_i. Some asymptotic properties of the mass imputation estimator given in (<ref>) under the assumed semiparametric model (<ref>)  are presented in the next section.










§ MAIN THEORY



We now discuss asymptotic properties of the mass imputation estimator given in (<ref>). <cit.> investigated the asymptotic properties of the mass imputation estimator in the context of combining two independent probability samples. We now consider an extension to mass imputation using a non-probability sample as a training sample. 
For the asymptotic framework, we assume there is a sequence of finite populations and a sequence of samples as discussed in <cit.>, which provides a framework to allow the sample sizes go to infinity. Let  be the solution to (<ref>) and define ^* = plim, where the reference distribution is the sampling mechanism for sample B. 
Roughly speaking, if (<ref>) holds, then ^*= _0, where   _0 are the true parameters in the conditional model  (<ref>).
Under certain regularity conditions including   
V {Û ( ^*) } =O_p (n_B^-1 ) , 
we can establish that 

    - ^* = O_p (n_B^-1/2 ).
 
  The following theorem presents an asymptotic expression of the mass imputation estimator. 








 Suppose that (,y) has bounded fourth moments over the sequence of finite populations  and that (<ref>) holds. Under the regularity conditions stated in the Supplementary Material,  
the mass imputation estimator (<ref>) satisfies
θ̂_I=θ̃_I+o_p(n_B^-1/2),
where

    θ̃_I=N^-1∑_i ∈ Aw_im(_i;^*)+N^-1∑_i ∈ B{ y_i-m(_i;^*) } g(_i; ^*)
 
with 
g(_i; ^*) =  (_i;^*)',  
and   is the   solution to 

    ∑_i ∈ B(_i;^*) (_i;^*)'  =  ∑_i=1^N(_i;^*) ,

where  (;) = ∂ m(;)/∂. 
 


Proof of Theorem <ref> is presented in the online Supplementary Material.  
The expression (<ref>) is essentially the first-order Taylor linearization of θ̂_I in (<ref>) and the Taylor expansion is made around = ^*, not around the true parameter value _0. Note that we do not use  the ignorability  assumption in (<ref>) to establish (<ref>). 


Following (<ref>), we can express 

    θ̃_I - θ_N = N^-1( ∑_i ∈ Aw_im_i^* -  ∑_i=1^N m_i^* ) +  N^-1(   ∑_i =1^N  δ_i g_i^*e_i^*   -  ∑_i=1^N e_i^*  )  ,

where m_i^* = m(_i; ^*),  e_i^* = y_i - m(_i; ^*) and g_i^* = g(_i; ^*). 
Thus, ignoring the smaller order terms, the asymptotic bias 
is 

    E( θ̃_I - θ_N)      =    E{ N^-1(   ∑_i =1^N  δ_i g_i^*e_i^*   -  ∑_i=1^N e_i^*  ) }
       =    -   E(  N^-1∑_i=1^N e_i^*  ) 
       =     E{ Cov_N ( δ_i π_B^-1 , e_i^* )  } ,
 
where π_B = n_B/N and the second equality follows from E( δ_i e_i^* |_i) =0  by the definition of ^* and e_i^*.  Here, we used the notation  Cov_N (x_i, y_i) = N^-1∑_i=1^N (x_i - x̅_N) (y_i - y̅_N). 
If  ^* = _0, then e_i^*=y_i - m(_i; _0)= e_i and 
the mass imputation estimator in (<ref>) is unbiased. 
Otherwise, the bias is non-zero. 





On the other hand, the asymptotic bias of the naive estimator, the simple sample mean y̅_B = n_B^-1∑_i ∈ B y_i,  is given by 

    E( y̅_B - θ_N )     =    E{ n_B^-1∑_i =1^N δ_i (y_i - Y̅_N)  }
       =    E{ Cov_N (  δ_i π_B^-1  , y_i) } .
  
 Thus, comparing (<ref>) with (<ref>), we find that the absolute value of the bias of θ̃_I is smaller than that of y̅_B as the variance of e_i^* will be smaller than the variance of y_i. Thus, the mass imputation estimator reduces the bias even when the sampling mechanism for sample B is non-ignorable. 
 
 We now derive the asymptotic variance formula for the mass imputation estimator. By (<ref>), 
 we have 
V(θ̃_I-θ_N) = V_A + V_B, where

    V_A =  V( N^-1∑_i ∈ A w_i m_i^* - N^-1∑_i=1^N m_i^* )

 is the variance component for sample A and
 
    V_B   =     V{  
    N^-1 (   ∑_i =1^N  δ_i g_i e_i^*   -  ∑_i=1^N e_i^*  ) }
       =    V[  E { 
    N^-1∑_i =1^N  ( δ_i g_i -1) e_i^*   |,  }   ] + E [  V { 
    N^-1∑_i =1^N  (δ_i g_i-1) e_i^* |, }  ]
 
is the variance component for sample B. In (<ref>), only the design-based variance under the probability sampling design for sample A remains. 
The reference distribution in (<ref>) is the joint distribution of the sampling mechanism for sample B and the  superpopulation model in (<ref>). The conditional expectation E{·|,  } is with respect  to the model for  e_i^* = y_i - m( _i; ^*) conditional on . The first term of (<ref>) is unknown in general because the sampling mechanism for sample B is unknown. 


If the sampling mechanism for sample B is ignorable as defined in (<ref>), we have ^* = _0 and e_i^* = e_i. In this case, the first term of (<ref>) disappears  because E ( e_i |_i, δ_i ) = E{ y_i - m( _i; _0) |_i, δ_i }=E{ y_i - m( _i; _0) |_i }=0. Thus, under the ignorability assumption, we can simplify V_B as 
 
    V_B   =      E [  V { 
    N^-1∑_i =1^N  (δ_i g_i-1) e_i |, }] 
       =    E{ N^-2∑_i=1^N (δ_i g_i-1)^2 e_i^2 } .
 
The second equality follows by the independence of e_i's in the superpopulation model. 
Thus, the first term of (<ref>) can be understood as the additional variance increase under  model misspecification. 
Furthermore, if n_B/N=o(1) then we can simply use

    V_B= E [ N^-2∑_i ∈ B(  e_i g_i )^2 ].



Note that we can easily estimate V_B in (<ref>) even though the sampling mechanism for sample B is unknown.

The asymptotic variance V(θ̃_I-θ_N) consists of two parts. The first term V_A is of order O(n_A^-1) and the second term V_B  is of order O(n_B^-1). If n_A/n_B=o(1), i.e., the sample size n_B is much larger than n_A, the term V_B is of smaller order and the leading term of the total variance is V_A. Otherwise the two variance components both contribute to the total variance. 




Under the linear regression model Y_i = 𝐱_i'  + e_i with e_i ∼ (0, σ_e^2), independent among all i, we have  ŷ_i = 𝐱_i' where
= ( ∑_i ∈ B𝐱_i 𝐱_i' )^-1∑_i ∈ B𝐱_i y_i.
The mass imputation estimator of (<ref>) under the regression model is given by θ̂_I,reg = N^-1∑_i∈ Aw_i_i'.
If the probability sample A is selected by simple random sampling, the asymptotic variance of θ̂_I,reg is given by

    V ( θ̂_I, reg - θ_N ) ≈ V ( n_A^-1∑_i ∈ A𝐱_i' )   + V ( N^-1∑_i ∈ B e_i _i' )  ,

where
= ( N^-1∑_i ∈ B𝐱_i 𝐱_i' )^-1_N and _N = N^-1∑_i=1^N 𝐱_i.
If 𝐱_i = (1, x_i)' and = (β_0,β_1)' and assuming π_B = n_B/N is negligible,  the asymptotic variance reduces to

    V ( θ̂_I, reg - θ_N )    ≈   1/n_A (β_1 )^2σ_x^2+ 1/n_Bσ_e^2  +  E[ (x̅_N -x̅_B)^2 /∑_i ∈ B(x_i-x̅_B)^2]  σ_e^2 ,

where σ_x^2 = V(X)  and x̅_B = n_B^-1∑_i∈ Bx_i.  If sample B is a random sample from the population, then the third term is of order O(n_B^-2) and becomes negligible. However, since sample B is a non-probability sample, the third term might not be negligible.




Variance estimation for the mass imputation estimator (<ref>) requires estimation of the two components V_A and V_B. The first component can be estimated by


    V̂_A=1/N^2∑_i ∈ A∑_j ∈ Aπ_ij-π_iπ_j/π_ij w_im(_i;)w_jm(_j;)  ,

where  π_ij = P(i,j∈ A) are the joint inclusion probabilities and are assumed to be positive.

The second component can be estimated by

    V̂_B= 1/N^2∑_i ∈ Bê_i^2{' (_i;)}^2  ,

where ê_i=y_i-m(_i;) and
=  ( N^-1∑_i ∈ B𝐱_i 𝐱_i' )^-1 N^-1∑_i ∈ A w_i  𝐱_i.
The total variance of θ̂_I can be estimated by V̂ ( θ̂_I - θ_N ) = V̂_A+V̂_B.








§ BOOTSTRAP VARIANCE ESTIMATION




The variance estimator presented in Section 4 is based on the linearization method. The closed-form formula for the asymptotic variance is simple to implement. However, to compute V̂_B in (<ref>), we need to use individual observations of (_i, y_i) in sample B, which is not necessarily available when only the sample A with mass imputed responses is released to the public data users.
Note that one of the goals of mass imputation is to produce a representative sample A with synthetic observations on the response variable using sample B as a training dataset. Once the mass imputation is performed, the training data is no longer necessary in computing point estimators. It is therefore desirable to develop a variance estimation method that does not require access to observations in sample B.


To achieve this goal, we propose a bootstrap method <cit.> for variance estimation that creates a replicated set of synthetic data {ŷ_i^(k),
 i  ∈ A } corresponding to each set of bootstrap weights { w^(k)_i , i  ∈ A }, k=1,⋯,L associated with sample A only. The method enables users to correctly estimate the variance of the mass imputation estimator θ̂_I without access to the training data {(y_i , _i ) : i ∈ B } from sample B. The data file will contain additional columns of {ŷ^(k)_
i : i ∈ A} associated with the columns of  bootstrap weights {w^(k)_i ;  i  ∈ A }, k =1, ⋯ , L, where L is the number of replicates
created from sample A. 
<cit.> also considered a similar method in the context of survey integration from non-nested two-phase sampling.


In order to develop a valid bootstrap method for the mass imputation estimator θ̂_I in (<ref>), it is critical to develop a valid bootstrap method for estimating V(  ) when  is computed from (<ref>).  Note that, under assumptions (<ref>) and (<ref>), we can obtain

    V (  ) ≐  J^-1Ω J^-1'

where J =  E{ n_B^-1∑_i ∈ B_i _i'  }, Ω = E{ n_B^-2∑_i ∈ B E( e_i^2 |)  _i _i' }
with _i =  ( 𝐱_i; _0 ) and _i =  ( _i; _0). The reference distribution for (<ref>) is the joint distribution of the superpopulation model (<ref>) and the unknown sampling mechanism for the non-probability sample B. Interestingly, the variance formula in (<ref>) equals  exactly to the variance of  when  sample B is selected by simple random sampling (SRS). That is, even though the sampling design for sample B is not SRS, its effect on the variance of  is essentially the same with SRS. This is due to the MAR assumption in (<ref>) which makes the effect of the sampling design for estimating  ignorable even though it is still not ignorable for θ= E(Y). Therefore, we can safely ignore the sampling design for sample B when estimating  and develop a valid bootstrap method for variance estimation of  using the bootstrap method for SRS.

Our proposed bootstrap method can be described as the following four steps:


  Step 1. Create the kth set of replication weights {w_i^(k), i∈ A} based on the sampling design for the probability sample A.

  Step 2. Generate the kth bootstrap sample of size n_B from sample B using simple random sampling with replacement and compute ^(k) using the same estimation equations (<ref>) applied to the bootstrap sample.

  Step 3. Use ^(k) obtained from Step 2 to compute ŷ_i^(k) = m( _i; ^(k)) for each i ∈ A and create a new column {ŷ_i^(k),i∈ A} alongside {w_i^(k), i∈ A} for the sample A dataset.

  Step 4. Repeat Steps 1-3, independently, for k=1,⋯, L for a pre-chosen L.
 

 Step 1 needs to follow standard practice in survey sampling on creating replication weights for design-based variance estimation. 
 See, for instance, Chapter 10 of <cit.>. 
 
 The final sample A dataset contains additional columns for the replicate versions of mass imputed values {ŷ_i^(k),i∈ A} and the sets of replication weights {w_i^(k), i∈ A}, k=1,⋯,L.  The replicate versions of the mass imputation estimator θ̂_I are computed as

    θ̂_I^(k) =1/N∑_i ∈ A w_i^(k)ŷ_i^(k) ,        k=1,⋯,L ,

and the resulting bootstrap variance estimator of θ̂_I  is computed as

    V̂_b ( θ̂_I ) = 1/L∑_k=1^L ( θ̂_I^(k) - θ̂_I )^2.

Note that once { (w^(k)_i ,   ŷ^(k)_
i ): i ∈ A} are obtained for k=1, ⋯, L, in addition to the original mass imputation data { (w_i ,   ŷ_
i ): i ∈ A},  the users do not need to have access to sample B. 


The following theorem establishes the consistency of the proposed bootstrap variance estimator.
Proof of the theorem is presented in the online Supplementary  Material.




Suppose that the assumptions of Theorem 1 hold. Assume that the sampling mechanism for sample B is ignorable as defined in (<ref>). 
Under the additional assumption stated in the Supplementary Material document, the bootstrap variance estimator given by (<ref>)
satisfies

    V̂_b(θ̂_I)=V(θ̂_I - θ_N)+o_p(n_B^-1)  .









§ SIMULATION STUDY



A limited simulation study is performed to (i) evaluate the departure of the imputation model to the performance of mass imputation estimator, (ii) check the validity of the proposed variance estimators, and (iii) compare the mass imputation estimator with the inverse probability weighted (IPW) estimator of <cit.> using a logistic regression model for estimating the propensity scores. 



 The setup for the simulation study employed  a 3 ×  2 factorial structure with two factors. The first factor is the superpopulation model that generates the finite population. The second factor is  the sample size  for sample B.  We used two levels for n_B, where  n_B=500 or n_B=1,000. 
 We generated
 the following three models for finite  populations of size N=100,000. The variables, x_i and e_i, are independently generated from N(2,1) and N(0,1), respectively. 
 
 The study variable y_i are constructed differently for each model: 
 
 
   Model I:    y_i = 1 +2x_i + e_i. 
 
   Model II:   y_i = 3 + x_i + 2e_i. 
 
   Model III:   y_i = 2.5 + 0.5x_i^2 + e_i. 
 
 Model I generates a finite population with a high correlation between x and y (r^2 = 0.8), Model II generates a finite population with a low correlation (r^2 = 0.2), and Model III generates a finite population where the linear relationship fails. Model III is included to check the effect of model mis-specification for the imputation model.


From each of the three populations, we generated two independent samples. We use  simple random sampling of size n_A=500 to obtain  sample A. In selecting sample B of size n_B, where n_B ∈{500, 1000 }, we create two strata where Stratum 1 consists of elements with x_i≤ 2 and Stratum 2 consists of elements with x_i>2. Within each stratum, we select n_h elements by simple random sampling, independent between the two strata,  where n_1=0.7 n_B and n_2=0.3 n_B. We assume that the stratum information is unavailable at the time of data analysis.
Using the two samples A and B, we compute four estimators of θ_N = N^-1∑_i=1^N y_i:





1) The sample mean  from sample A:   θ̂_A=n_A^-1∑_i ∈ A y_i.



2) The naive estimator (sample mean) from sample B:   θ̂_B = n_B^-1∑_i ∈ B y_i.



3) The mass imputation estimator from sample A  given in (<ref>) using ŷ_i = β̂_0 + β̂_1 x_i where (β̂_0, β̂_1) are  the estimated regression coefficients obtained from sample B.



4) The inverse probability weighted (IPW) estimator proposed by <cit.>:
θ̂_IPW=N^-1∑_i ∈ Bπ̂_i^-1y_i,
where the propensity scores, π_i=π(_i;) = {1+exp(-ϕ_0 - ϕ_1 x_i ) }^-1 with _i = (1,x_i)' and = (ϕ_0,ϕ_1)', are estimated by using  which solves the following score equations:

    U()=∑_i ∈ B_i-∑_i ∈ Aw_iπ(_i;)_i=0.



The sample mean of sample A, not available in practice but computed  in the simulation, serves as a gold standard estimator. Results are based on 1,000 repeated simulation runs.
Table <ref> presents the Monte Carlo bias, the Monte Carlo variance, and the relative mean squared error of the four point estimators. The relative mean squared error  of estimator θ̂ is defined as

    ReMSE =MSE(θ̂)/MSE(θ̂_A) .

Note that the sample mean θ̂_A from sample A is  not available in practice but is computed here as the gold standard. Table <ref> shows that, with models I and  II where the linear regression model holds, the mass imputation estimator is unbiased for the population mean.  The naive mean estimator of sample B underestimates the population mean for all scenarios considered in the simulation. When the size of sample B for training data is larger than the size of sample A (n_B=1,000), it is possible that the mass imputation estimator has a smaller MSE than the gold standard. Under the simple regression model,
the asymptotic  variance of the mass imputation estimator is

    V(θ̂_I - θ_N) ≈1/n_Aσ_x^2 β_1^2 + 1/n_Bσ_e^2 E  { 1 +  (x̅_N - x̅_B)^2/ s_x,B^2 } ,

where s_x,B^2 = (n_B-1)^-1∑_i∈ B(x_i - x̅_B)^2, while the variance of sample mean of sample A  is V(θ̂_A)=σ_y^2/n_A=(β_1^2σ_x^2+σ_e^2)/n_A. Thus, if n_B is much larger than n_A, the mass imputation estimator can be more efficient than the sample mean of A. The IPW estimator is less efficient than the mass imputation estimator for all cases considered in the current simulation setup.

The imputation model using the simple linear regression is incorrectly specified for model III. The mass imputation estimator is modestly biased. Nonetheless, the performance in terms of MSE is better than the IPW estimator because the mass imputation estimator has much smaller variance than the IPW estimator, even though the absolute bias is larger, when the linear relationship fails.












Table <ref> presents Monte Carlo mean and relative bias of the two variance estimators of the mass imputation estimator using linearization and bootstrap.  Both variance estimators show negligible relative biases. In particular, the bootstrap variance estimator shows good performance even under model III.





§ APPLICATION TO THE PRC DATASET



We now apply the proposed mass imputation techniques described in Section <ref> and compute point and variance estimates of population means  for 
the four study variables in the PRC dataset, which were introduced in Section 2:

    
  Y_1: Talk with neighbors frequently, 
    
  Y_2: Participated in school groups, 
    
  Y_3: Participated in service organizations, 
    
  Y_4: Days had at least one drink last month.

Note that PRC is the non-probability sample and CPS and BRFSS are the two existing probability samples to be used to provide information on  auxiliary variables. The PRC sample was collected through web panels, the BRFSS survey was conducted by telephone interviews, and the CPS survey used a mixed mode with 11 percent of units surveyed by telephone and the remaining part of the survey by personal interviews. Further details are provided in the online  Supplementary Material document. 

Prior to any formal analysis of the datasets, we would like to check whether the positivity and the transportability assumptions are reasonable. The support of the auxiliary variables in the PRC dataset can be checked and compared to the support in the CPS sample or the BRFSS sample, which confirms the positivity assumption. The transportability assumption cannot be formally tested. We first looked into the auxiliary variables to see if they are meaningful predictors for the study variables, and the answer is clearly affirmative. We then conducted a sensitivity analysis by (i) choosing a “test study variable” Z which is available from the non-probability sample and one of the probability samples; (ii) fitting the test model f(z |) separately by using the non-probability sample and the probability sample; (iii) comparing the two fitted models for similarities and differences. Details of the sensitivity analysis are reported in the Supplementary Material document. Note that results from such analysis do not formally prove or disprove the validity of the assumption. For instance, the transportability of the model f(x_2 | x_1) does not imply the transportability of the model f(y| x_1,x_2).


Our analysis focuses on three important practical aspects of the method. We first investigate how the relative sizes of sample A and sample B impact the mass imputation estimator and  its associated variance estimator. 
We then discuss the covariates selection for the prediction model, and apply different strategies based on the availability of auxiliary variables from one or both probability samples. Lastly, we make comparisons between the mass imputation estimators and the IPW estimators.



 §.§ Impact of relative sample sizes


Recall from Section 2 that the two probability samples CPS and BRFSS contain respectively n_A = 80,075 and 441,456 cases while the non-probability PRC  dataset has n_B=9,301 units. The sample size n_A  is much larger than n_B for both CPS and BRFSS. To investigate other possible scenarios where n_A  and n_B has different ratios, we draw three subsamples from the original BRFSS dataset by the simple random sampling without replacement method. The resulting subsamples, denoted by  BRFSS(1), BRFSS(2) and BRFSS(3), have sample size  n_A^*=80,000, 8,000 and 800, respectively. 
Survey weights for each of the subsamples are computed as w_i n_A/n_A^*, where w_i is the weight  of unit i in the original  BRFSS sample with n_A=441,456.

The choices of covariates for the prediction model are constrained by the set of covariates observed in both the probability sample and the non-probability sample. The common set of covariates between  BRFSS and PRC,  however, differs from the set between CPS and PRC. We use the set of covariates which are available in all three datasets for the current investigation. We use a logistic regression model for  binary responses and a linear regression model for the continuous response. 

We first compute the point estimates of the population means by the proposed mass imputation method using five different probability samples, and the results are presented in Table <ref>. The first row specifies which probability sample is used as sample A except for the column under “PRC” which represents the naive estimates of simple sample means.
It can be seen that the three larger probability samples BRFSS, BRFSS(1) and BRFSS(2) with n_A=441,456 and n_A^*=80,000  and 8,000 produce almost identical results. The smallest probability sample BRFSS(3) with n_A^*=800 leads to noticeably different results, indicating potential inconsistent estimates when the size of the probability sample is too small. We also notice that, despite the differences between the CPS and the BRFSS samples, their corresponding  mass imputation  estimators  produce similar results 
with the use of the same set of auxiliary variables.

We further look at  variance estimators with different probability samples, using the linearization and bootstrap methods described in Sections 4 and 5. The linearization variance estimator  is based on the formula V(  θ̃_I - θ_N) = V_A + V_B given by Theorem 1, where V_A is the design-based variance component under the probability sampling design for sample A, and  V_B is the variance component for sample B under the prediction model. Unfortunately, detailed design information other than the survey weights is not available for either the BRFSS or the CPS sample. We use an approximate variance formula for V_A by assuming that the survey design is single-stage PPS sampling with replacement, a strategy often used by survey data analyst for the purpose of variance estimation. The bootstrap variance estimator  is computed based on L=5,000 bootstrap samples.

Results for variance estimators with decomposition of the variance components (V̂_A, V̂_B) for the linearization method are reported in Table <ref>.  The variance estimates  have been multiplied by 10^5 for binary response variables and  by 10^2 for the continuous variable. Variance decomposition cannot be done for the bootstrap method. We have the following major observations from Table <ref>: (i) the two variance estimation methods, linearization and bootstrap, give similar results in total variance for all cases; (ii) the two original large probability samples CPS and BRFSS produce similar total variances for all cases, with the design variance component V̂_A making a negligible contribution; (iii) when the size of the probability sample becomes smaller, from BRFSS(1) to BRFSS(2) and then BRFSS(3), the total variance becomes larger, and the variance component V̂_A from the probability sampling design for sample A becomes dominant. 




 §.§ Impact of the prediction model
 

The second part of the analysis is on covariate selection for the prediction model. We consider following selection strategies: (a) Use the common set of variables available in all three datasets. This is what has been used in Section 7.1, and the resulting model is denoted as ξ( Partial). (b) Use the common set of variables available in the two datasets PRC and CPS or PRC and BRFSS, leading to two larger but different common sets of variables for the two probability samples. The resulting models are denoted as ξ( All). (c) Use the set of covariates which are available in both PRC and CPS (or both PRC and BRFSS) but are selected through the backward variable selection algorithm. The resulting models are denoted as ξ( Select). P-values for covariates in  model ξ( All) and ξ( Select) are listed in Table <ref> and Table <ref>.

Mass imputation based point estimates θ̂_I for the population means of the four response variables using different sets of covariates with the probability Sample A as either CPS or BRFSS are presented in Table <ref>, with the linearization variance estimates displayed in the parentheses (multiplied by 10^5 for binary variables and by 10^2 for the continuous variable). Major observations from Table <ref> can be summarized as follows. (i) With a chosen sample A, the point estimates under model ξ( All) are very similar to the estimates under model ξ( Select) which drops some covariates using the backward variable selection algorithm; (ii) the point estimates under model ξ( All) behave quite differently to the estimates under model ξ( Partial) which excludes some covariates based on the availability from a second probability sample; (iii) the linearization variance estimates under model ξ( Select) are always smaller than the variance estimates under model ξ( All), showing some efficiency gain by eliminating non-significant factors from the model; (iv) the point estimates using CPS are very similar to the estimates using BRFSS under the same model ξ( Partial) with the same set of covariates; (v) the point estimates using CPS under model ξ( All) are quite different to the estimates using BRFSS under model ξ( All). Note that the two models use two different sets of covariates. Some additional covariates available with one particular probability sample may be highly correlated to some of the response variables. 
For instance,  the  covariate “Volunteer works”, which is  available in CPS but not in BRFSS, is related to the response variables Y_1, Y_2 and Y_3. Similarly, health related covariates, which are available in BRFSS but not in CPS, likely explain the evident difference of θ̂_I for response  Y_4 when the BRFSS probability sample is used. 




 §.§ Comparing mass imputation with inverse probability weighting
 

In the final part of the analysis of the PRC non-probability sample, we make comparisons between mass imputation estimator θ̂_I and IPW estimator θ̂_IPW of <cit.>, which is also examined in the simulation study reported in Section 6. The mass imputation estimates are computed based on model ξ( Select). The IPW estimators are based on the logistic regression model for the propensity scores using the same set of covariates as in the model ξ( Select). 
The point estimates along with the associated linearization variance estimates (multiplied by 10^5 for binary variables and by 10^2 for the continuous variable) for the four response variables are reported in Table <ref>. The two estimators θ̂_I and θ̂_IPW provide   comparable results in both bias and variance, regardless whether CPS or BRFSS is used as the probability sample A. 







§ CONCLUDING REMARKS



The use of non-probability survey samples as an efficient and cost-effective data source has become increasingly popular in recent years. Theoretical developments on analysis of non-probability samples, however, severely lag behind the need of making valid inference from such datasets. Non-probability survey samples are biased and do not represent the target population. Valid inferences require supplementary information on the population. The mass imputation approach proposed here relies on the availability of a probability survey sample from the same target population with information on covariates. The covariates are also measured for the non-probability sample and need to possess two crucial features for the framework discussed in this paper: 
(a) they characterize the inclusion/exclusion mechanism for units in the non-probability sample so that the ignorability or the transportability assumption could be justified; and (b) they are relevant to the response variable in terms of prediction power, a condition required for the mass imputation estimator.

Our analyses of the PRC non-probability sample shed light on two important practical issues. The first is the choice of an existing probability sample from the same target population if two or more such samples are available. The decision rests on the set of auxiliary variables which is available in both the non-probability and the probability samples. The size of the probability sample is less crucial as long as it is not too small. Two different probability samples tend to produce similar results if the same set of auxiliary variables is used. The second is the interpretation of the mass imputation estimator when the non-probability survey and the probability survey use different modes of data collection. The estimator is computed using the probability sample under mass imputation but the final estimate carries the mode from the non-probability sample since the model is built based on the observed study variable in the non-probability sample. 


One of the main advantages of the mass imputation approach over the propensity score approach is that we can use  subject matter knowledge on the study variables  to build a good prediction model that can lead to  better estimation.
Also, as discussed in Section 4, the bias of the mass imputation estimator can be small even in misspecified models. Thus,  mass imputation can provide  rather robust estimates. The price we pay is the extra modeling work for each survey item. 



The theoretical results on mass imputation presented in this paper focus on the estimation of finite population means. This is in line with traditional approaches in survey sampling where the basic theory is developed for the Horvitz-Thompson estimator. Inferences on other finite population parameters under mass imputation remain to be a research topic for future development. The main theoretical results are also based on deterministic imputation under the semiparametric model (<ref>). Another research topic is to explore the use of random imputation, including multiple imputation <cit.> and fractional imputation <cit.>, for analyzing non-probability survey samples.  


Statistics Canada has been implementing the modernization initiatives in recent years, which call for a culture switch from the traditional survey-centric approach by the agency to using data from multiple sources. One of the questions arising from the discussions is the role of traditional probability-based surveys, and there are even questions on the necessity of their existence in the future. Our theoretical results presented in this paper call for probability survey samples with rich information on auxiliary variables. A few large scale high quality probability surveys representing the target population can play significant roles in analyzing data from non-probability survey samples. Also,  statistical analysis of  combined or fused data is an important problem <cit.>. Further research on this direction will be presented  elsewhere. 



§ SUPPLEMENTARY MATERIAL
 

The online supplementary material provides proofs of the Theorems 1 and 2, results from a sensitive analysis, and additional information about the three samples used in the application. 



§ ACKNOWLEDGMENTS


The authors are grateful to  the two anonymous referees, the AE, and the Editor for their very constructive comments. The research of J.K. Kim was partially supported by a grant from US National Science Foundation (MMS-1733572). The research of C. Wu was supported by a grant from the Natural Science and Engineering Research Council of Canada. 



































dcu





