Data Integration by combining big data and survey sample data for finite population inference
Jae-kwang Kim Siu-Ming Tam
stat.ME





#1 1




1
  
  
  Data integration by combining big data and survey sample  data for finite population inference  
    Jae Kwang Kim
    

    Department of Statistics, Iowa State University

    and 

    Siu-Ming Tam 

    University of Wollongong and  Australian Bureau of Statistics
    December 30, 2023
===========================================================================================================================================================================

 

0
  

    Data integration by combining big data and survey sample  data for finite population inference  


   


The statistical challenges in using big data for making valid statistical inference in the finite population  have been well documented in  literature.  These challenges are due primarily to statistical bias arising from under-coverage in the big data source to represent the population of interest and measurement errors in the variables available in the data set.  By stratifying the population into a big data stratum and a missing data stratum,   we can estimate the missing data stratum by using a fully responding probability sample, and hence the population as a whole by using a data integration  estimator.
By expressing the data integration estimator as a regression estimator, we can handle  measurement errors in the variables in big data and also in the probability sample. We also propose a fully nonparametric classification method for identifying the overlapping units and develop a bias-corrected data integration estimator under misclassification errors. Finally, we  develop a two-step regression data integration estimator to deal with measurement errors in the probability sample.  An advantage of the approach advocated in this paper is that we do not have to make unrealistic missing-at-random assumptions for the methods to work. The proposed method is applied to the real data example using 2015-16 Australian Agricultural Census data.




Keywords:
Calibration weighting; 
Measurement error; Non-response;  Regression estimation;   Selection bias.






1.5 




§ INTRODUCTION



Suppose we are interested in estimating some finite population parameters, e.g. the finite
population mean, of a target population based on a data set.  If the data set comes from a probability sample,  parameter estimation is
straightforward,  and we can draw on the extensive literature on survey sampling over the past
century, e.g. <cit.>, <cit.>, <cit.>. However, if the data set comes from a
non-probability sample, e.g. from a big data
source, the estimation is less straightforward, and  the theory for making inference with non-probability
samples is not fully developed. <cit.> and
<cit.> addressed methodological uses and challenges of big data in the production of official statistics.








The perils and pitfalls in using big data are primarily under and over coverage, and self selection.  Bias from under coverage is akin to bias from non-random samples for inference, and the bias from self-selection is akin to nonresponse bias in surveys.  These biases have been discussed extensively in the statistics literature (see for example, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>).  Specific discussion of these biases can be found in <cit.>  for web data; <cit.> on data from smart phones and wearable devices; and <cit.> for social media data. The weighting methods considered in <cit.> and  <cit.>  are based on 
a missing-at-random  assumption (MAR)  of  <cit.>. 
The MAR assumption is a strong assumption and there is no way to verify this assumption from the data only. 





Survey data integration, which is developed to combine information for two independent surveys from the same target population, can be used to handle the selection bias of non-probability samples by incorporating a probability sample. <cit.>  proposed a mass imputation approach for survey integration. In <cit.>, the nearest neighbor matching imputation is used to identify the imputed value for each element in the probability sample. <cit.> developed a statistical theory for register-based statistics and data integration.
<cit.> discussed practical issues in sample matching for solving the selection bias in the non-probability sample.
While  matching-based imputation is promising and potentially useful in practice, it is still based on the missing-at-random  assumption.  <cit.> also considered a weighting  adjustment method based on parametric model assumptions on the  selection mechanism for the non-probability sample, but the MAR  assumption is still required.  <cit.> provided comprehensive reviews of statistical methods of data integration  for finite population inference.

In this paper, we propose a novel method of data integration for  handling big data by incorporating survey sample data. The sampling mechanism for big data  is not necessarily MAR. That is, there can be some systematic difference between the big data sample and the survey sample even after adjusting for the auxiliary variables.  We assume that the survey variables are observed in both samples, but allow them  to be inaccurately measured in one sample. Our approach is to treat the big data sample as a finite population of incomplete (or inaccurate) observations. Furthermore, the incomplete observations in the population can be treated as auxiliary information for calibration weighting <cit.>. Thus, standard  techniques such as calibration weighting for incorporating  auxiliary information from the finite population can be used directly. To conduct calibration estimation in the survey data, we need to identify the subset of the probability sample that also belongs to the big data sample. This  is somewhat similar in spirit to dual frame estimation <cit.>. In our application, the big data sample is subject to coverage errors, but the survey sample is not.   The proposed method is particularly useful for government statistical agencies which can effectively apply such matching.
 When the accurate matching is not possible, we propose a novel  classification method to identify the overlapping units using the matching variables observed from two data sources. Fully nonparametric propensity scores are obtained from the proposed classification procedure and  they  can be used to correct for the bias in applying the data integration estimator with inaccurate matching. 
 

 




The paper is organized as follows. In Section 2, basic setup is introduced. In Section 3, the basic idea for data integration is introduced.  In Section 4, a semi-supervised classification method is introduced to identify the overlapping units with big data. In Section 5, an efficient method for data integration is introduced. In Section 6, the proposed method is extended to the case of measurement errors in the sample observation.

Two  limited simulation studies are presented in Section 7 and an application of the proposed method to an official statistics is presented in Section 8. Some concluding remarks are made in Section 9.




§ BASIC SETUP



Consider a finite population U={1, ⋯, N} of size N.
From the finite population, we have two samples, denoted by A and B,
where A is  a probability sample and B is a big data sample obtained by an unknown selection mechanism. From both samples, we measure the study variable Y. 
Initially, we assume  that Y is measured without measurement error in sample A, but we shall relax this assumption in Section 6.
However, in  sample B, Y is not necessarily measured accurately. Thus, instead of observing y_i, we observe y_i^*, which is a contaminated version of y_i, from sample B.  For simplicity, we assume that

    y_i^* = β_0 + β_1 y_i + e_i,

where (β_0, β_1) is an unknown parameter and e_i ∼ (0, σ^2). Model (<ref>) implies that y^*_i can be systematically different from y_i. In the special case of (β_0, β_1)=(0,1),  there is no measurement bias in y^*_i. In addition, since the selection mechanism for the big data sample is unknown, it is subject to selection bias. Generally speaking, the selection bias of big data cannot be ignored, and adjusting for the selection bias is critical <cit.>.

To correct for the selection bias and measurement errors in the big data, we assume that we have a gold standard survey sample. Obtaining survey sample data is often expensive, but the gold standard can be used to improve the quality of the big data sample. Furthermore, optimal allocation of the resources can   make the final analysis more cost-effective.

To make sample A  a gold standard sample, a probability sampling design for selecting sample A is employed, and y_i are accurately observed from the sample. From sample A, we can compute T̂_a = ∑_i ∈ A d_i y_i, a design-unbiased estimator of T=∑_i=1^N y_i, where d_i = π_i^-1 is the design weight of unit i, and π_i is the first-order inclusion probability of unit i  in sample A.

Table 1 presents the data structure of our setup.
We also assume that it is possible to identify  elements in sample A  also belonging  to sample B.  That is, we can create δ_i for i ∈ A, where

    δ_i = {[     1 i ∈ B;     0       ].

Thus, we can observe δ_i in sample A if the individual-level matching is possible. We shall relax this assumption in Section 5. 




Our goal is to combine the observations in the two data sets  to find an improved estimator of T. By making a proper use of big data through weighting, we can obtain an improved estimator of T over T̂_a, which  completely ignores the information in the big data sample. Combining two data sources is called data integration,  and we will consider data integration as a general tool for making a proper use of big data for finite population inference.
Challenges in data integration are outlined in  <cit.> and <cit.>. 

<cit.> provided methods for adjusting such bias by using data integration. 
This paper extends the work  of <cit.>  to non-binary variables, and also addresses situations when there are measurement errors or matching errors  in the data sets.







§ DATA INTEGRATION  FOR HANDLING SELECTION BIAS



We first consider the simple case of no measurement errors in Y, i.e.,  y_i^* = y_i. Now, we can conceptually define δ_i in (<ref>) throughout the finite population. Thus, the set of elements with δ_i=1 is the big data sample. We can decompose

    T= ∑_i=1^N y_i = T_b + T_c,

where T_b = ∑_i =1^N δ_i y_i and T_c= ∑_i=1^N (1-δ_i) y_i. Since T_b can be obtained from sample B, we  only have to estimate T_c from  sample A. Thus, we can use

    T̂_DI = T_b + ∑_i ∈ A d_i (1-δ_i) y_i

as a design-based estimator of T obtained from two samples.  If the population size  N is known,  a better estimator is

    T̂_PDI = T_b + (N-N_b) ∑_i ∈ A d_i (1-δ_i) y_i/∑_i ∈ A d_i (1-δ_i)   ,

where N_b= ∑_i=1^N δ_i is the size of sample B.
Estimator T̂_PDI in (<ref>) is essentially a post-stratified estimator with the two post-strata defined by δ_i=1 and δ_i=0,  respectively.

The design variance of T̂_PDI in (<ref>) is 

    (  T̂_PDI ) = (N-N_b)^2 {∑_i ∈ A d_i (1-δ_i) y_i/∑_i ∈ A d_i (1-δ_i) }≈{∑_i ∈ A d_i (1-δ_i) (y_i - Y̅_c) } ,
 
where Y̅_c= ∑_i=1^N (1-δ_i) y_i/ (N-N_b). Here, the approximate equality  follows from Taylor  linearization applied to the ratio component in (<ref>).  
If the sampling design for sample  A is simple random sampling of size n with n/N≈ 0, we have

    (  T̂_PDI ) ≈ (1-W_b) N^2/n S_c^2,

where  W_b=N_b/N and  S_c^2 = (N-N_b)^-1∑_i=1^N (1-δ_i) (y_i - Y̅_c)^2.  Thus, the variance reduction of T̂_PDI compared with  T̂_a = ∑_i ∈ A d_i y_i is

    ( T̂_PDI ) / ( T̂_a ) =  (1- W_b)  S_c^2/ S^2  .

If S_c^2 ≈ S^2, the data integration estimator  is always more efficient than the design-based estimator using sample A only.
In fact, from (<ref>), the effective sample size  using the post-stratified data integration estimator is

    n^* = n  1/1-W_bS^2/S_c^2 .


Thus,  if we define c_a and c_b to be the per-unit cost of observing y_i in sample A and sample B, respectively,  the total cost function using post-stratified data integration estimation is
C_DI = c_a n + c_b N_b, while the total cost required to obtain the same efficiency of T̂_a  is
C_a = c_a n^*.
If S_c^2 ≈ S^2, we have

    C_DI - C_a = c_b N W_b - c_a n W_b/1-W_b .

Therefore, given the same efficiency, the cost for using post-stratified data integration estimator is lower than using sample A only  if

    c_b/c_a≤n/N 1/1-W_b.

Thus,  if the under-coverage rate of B is less than  (c_a/c_b)· (n/N) , the proposed data integration estimation  is cost-effective by  (<ref>).




§ EFFICIENT ESTIMATION


We now discuss how to further  improve the efficiency of the data integration estimator.  One approach is to  use the idea of ratio estimation for T by treating x_i = δ_i y_i as the auxiliary variable, which is observed throughout the finite population. Thus,

    R̂ = ∑_i=1^N x_i/∑_i ∈ A  d_i x_i

can be multiplied  to direct estimator to reduce the variance, that is, to improve efficiency. The resulting ratio estimator is

    T̂_RatDI  = T̂_a R̂ = T_b T̂_a/T̂_b ,

where T̂_b = ∑_i ∈ A d_i δ_i y_i and  T̂_a = ∑_i ∈ A d_i y_i.  Thus, T̂_RatDI in (<ref>) is called the ratio data integration estimator.  Note that we can express T̂_RatDI   as

    T̂_RatDI  = ∑_i ∈ A d_i ( T_b/T̂_b) y_i  = ∑_i ∈ A w_i  y_i ,

where w_i satisfies

    ∑_i ∈ A w_i x_i = ∑_i ∈ A d_i ( T_b/T̂_b) δ_i y_i = ∑_i=1^N δ_i y_i =    ∑_i=1^N x_i .

Thus, equality (<ref>) implies that the ratio data integration estimator satisfies the calibration property of the auxiliary variable in the sense that the estimator applied to x_i  matches  the known population total of x_i.


More generally, we can  apply the calibration estimation method  to  _i = (1, δ_i, δ_i y_i)^, since ∑_i=1^N (1, δ_i, δ_i y_i ) = (N, N_b, T_b) is known.  Specifically, we can find {w_i : i∈ A} that minimizes an  objective  function Q(d, w) subject to the calibration equation ∑_i ∈ A w_i _i = ∑_i =1^N _i. The  regression estimator is based on

    Q(d, w) = ∑_i ∈ A d_i ( w_i/d_i -1 )^2 .

The solution to the optimization problem is

    w_i = d_i  𝐗_N^( ∑_i ∈ A d_i _i _i^)^-1_i ,

where 𝐗_N = ∑_i=1^N _i.

To understand the solution in (<ref>), if we write _i =(1-δ_i, _1i^)^ with _1i = δ_i (1, y_i)^, the regression weight in (<ref>) reduces to

    w_i = {[ d_i 𝐗_1^_xx11^-1_1i             δ_i = 1;      d_i (N_c/N̂_c)           δ_i = 0 , ].

where 𝐗_1 = ∑_i=1^N _1i,
_xx11 = ∑_i ∈ A d_i _1i_1i^, N_c= N-N_b  and N̂_c = ∑_i ∈ A d_i (1-δ_i). The weights in (<ref>) satisfy

    ∑_i ∈ A w_i (δ_i, δ_i y_i)= (N_b, T_b),   ∑_i ∈ A w_i (1-δ_i)  = N_c .

The regression data integration estimator is then defined as

    T̂_RegDI = ∑_i ∈ A w_i y_i ,

where w_i is defined in (<ref>).  Inserting (<ref>) into (<ref>), we can write

    T̂_RegDI = ∑_i =1^N δ_i (1, y_i)^_1 + N_c T̂_c/N̂_c ,

where T̂_c = ∑_i∈ Ad_i(1-δ_i)y_i and

    _1 ={∑_i ∈ A d_i δ_i (1, y_i) (1, y_i)^}^-1∑_i ∈ A d_i δ_i  (1, y_i)^ y_i  =
    (0,1)^.

Therefore, the regression data integration estimator in (<ref>) is algebraically equivalent to the post-stratified data integration estimator in (<ref>). However, we can include other auxiliary variables observed throughout the finite population in the calibration equation; see   Remark 1 below for details.


For variance estimation, standard linearization methods or replication methods for regression estimator can be applied. For example, a linearization variance estimator for (<ref>) can be written as

    V̂ ( T̂_RegDI) = ∑_i ∈ A∑_j∈ Aπ_ij - π_i π_j/π_ijê_i /π_iê_j /π_j,

where π_ij is the joint inclusion probability of unit i and j, ê_i = y_i - _i^ and 

=( ∑_i ∈ A d_i _i _i^)^-1∑_i ∈ A d_i _i y_i. 




	In addition to y_i, if there is another variable z_i observed in both samples, we can incorporate this information into calibration estimation. That is, we use _i = (1-δ_i, δ_i, δ_i y_i, δ_i z_i)^ in the calibration estimation.  If z_i is observed throughout the finite population, we can use _i = (1-δ_i, δ_i, δ_i y_i,  z_i)^.
	


In some cases, the big data may have duplication and lead to over-coverage problems. In this case, we can still apply the idea of calibration estimation by modifying the definition of δ_i to be the  number of times that the unit  appears in sample B. In this case, we can use

    ∑_i ∈ A w_i ( 1, δ_i, δ_i y_i ) = ∑_i ∈ U (1, δ_i, δ_i y_i )
 
as the calibration equation.



 
The  proposed method is also applicable  when  measurement errors exist in addition to selection bias in big data sample. That is, instead of observing y_i, we observe y_i^*, an inaccurate measurement for y_i, in sample B. In sample A, in addition to  observing (y_i, δ_i), we assume that  it is possible to obtain y_i^* for units with δ_i=1 by matching. Thus, we observe (y_i, δ_i, δ_i y_i^*) in sample A.
In this case, we can still use δ_i y_i^* as a control for the calibration equation. Thus, instead of using _i= (1- δ_i, δ_i, δ_i y_i)^, we can use _i^*=(1-δ_i, δ_i, δ_i y_i^*)^ in (<ref>) to get the calibration weights satisfying ∑_i ∈ A w_i (1-δ_i) = N_c, ∑_i ∈ A w_i δ_i = N_b and ∑_i ∈ A w_i δ_i y_i^* =   ∑_i ∈ B y_i^*. 









§ SEMI-SUPERVISED CLASSIFICATION

 


The proposed method in Section 3 is based on the assumption that the big-data indicator function δ_i is observed for every element in sample A. If we have an access to the unique identifiers then it is possible to match the records accurately  and obtain δ_i. In other cases, we only have matching variables such as name, zip code, and date of birth, etc. In this case, we use these matching variables to obtain the best guess of δ_i, denoted by δ̂_i, based on the observed value of the matching variables 𝐳_i. Obtaining δ̂_i from the matching variables is a challenging classification problem. Furthermore, finding a bias-corrected estimator under misclassification error is not considered in the literature. In the context of multiple frame surveys, <cit.> developed a bias-adjustment method assuming that the misclassification probabilities are known. 



In our setup, note that δ_i is observed for sample B, as δ_i=1 if i ∈ B by definition. We do not observe δ_i  if  i ∈ A. Thus, this is a semi-supervised classification problem because the true label (δ_i) for classification is available only for sample B. Here, we shall propose a maximum likelihood method of semi-supervised classification under the setup of data integration. Note that unlike the probabilistic record linkage, we do not have to identify the pairs of matches and non-matches as in <cit.>. We have only to determine  whether each unit in sample A  belongs to the particular subpopulation B or not. 

To formally describe the idea of the proposed method, recall that 
the finite population U is  decomposed into two groups, U=B ∪ B^c.  We assume that π = P( δ = 1) is known and given by π = N_b/ N.  We have a probability sample A selected from U and observe _i instead of observing δ_i for all i ∈ A. 
 
 If the densities for two groups, p( |δ =1) and p( |δ = 0),  are known or estimated from the samples, then we can use 

    P( δ_i=1 |_i ) = π p( _i |δ_i  = 1) / (1-π) p( _i |δ_i = 0) + π p( _i |δ_i = 1)

to make classification for unit i ∈ A. We use δ̂_i=1 if we classify unit i as i ∈ B.  Otherwise, we use δ̂_i=0. The decision rule is 

    δ̂_i = 1 P̂ ( δ_i=1  |_i ) > 1/2 ,
 
where 

    P̂( δ_i=1 |_i ) = πp̂( _i |δ_i  = 1) / (1-π) p̂( _i |δ_i = 0) + πp̂( _i |δ_i = 1)  .

 Note that p( |δ =1) means the marginal density function of  among big data. Estimation of p( |δ =1) is straightforward as long as we have access to the big data. Thus, we have only to estimate  parameters in p( |δ=0). 
 


To discuss parameter estimation, suppose that =(z_1, ⋯, z_K) and each z_k can take one of D values among the set 𝒵_k= { z_k^(1) , ⋯, z_k^(D)} with unknown probabilities. We assume that 
    
    p( |δ=1) = ∏_k=1^K p_k ( z_k |δ = 1)
 
    where 
p_k (z_k |δ =1) = m_kd   z_k = z_k^(d)
    and     ∑_d=1^D m_kd = 1. Since 
     we can observe _i among δ_i=1, we can estimate m_kd using 
    
    m̂_kd = 1/N_B∑_i ∈ BI  ( z_ik = z_k^(d)) .

Now, for the model for p( |δ=0), we assume that 
    
    p( |δ=0) = ∏_k=1^K p_k ( z_k |δ = 0) ,

    where 
    p_k (z_k |δ =0) = u_kd   z_k = z_k^(d) 
    and     ∑_d=1^D u_kd = 1.
 If we define 
    
    γ_ik^(d) = 
    {[              1 z_ik = z_k^(d);              0                ].

    then we can express 
m_kd = P ( γ_ik^(d) = 1 |δ_i=1) and 
u_kd = P ( γ_ik^(d) = 1 |δ_i=0) . 



 To estimate u_kd, we use the following EM algorithm: 
	 
		
  * First note that,  if δ_i were observed, the complete-sample pseudo log-likelihood would be 
		
    l_com (  𝐮|,  )   =   ∑_i ∈ A d_i  δ_i log{π∏_k=1^K  m_ik}+  ∑_i∈ A d_i  (1- δ_i)  log{ (1-π )  ∏_k=1^K   u_ik}
 
		where 
    ( m_ik, u_ik) =  ∑_d=1^Dγ_ik^(d) (m_kd , u_kd) 
    .


		Note that there is no need to estimate m_kd again, because we have access to big data directly. Only u_kd are the parameters of interest. 
		
  * In the E-step, we need to evaluate the conditional expectation of l_com (  𝐮|, ) given the observed data. Thus, given the current parameters, we have only to compute 
		
    Q(  𝐮|𝐮^(t)  )    =    E {l_com (  𝐮|, ) |𝐮^(t)}
       =   ∑_i∈ A   d_i p̂_i^(t)log{π∏_k=1^K    m_ik} + ∑_i∈ A d_i  (1- p̂_i^(t) )  log{  (1-π )  ∏_k=1^K    u_ik} ,
 
		where 
		
    p̂_i^(t)   =    E ( δ_i|_i ;  𝐮̂^(t)    ) 
       =   π∏_k=1^K   m_ik/π∏_k=1^K   m_ik +(1- π)  ∏_k=1^K   û_ik^(t)
 
		and 
		(m_ik , û_ik^(t)  ) = ∑_d=1^Dγ_ik^(d) ( m_kd, û_kd^(t) ).  
		
  * The M-step is to maximize the Q over 𝐮 to update the parameters. The updating formula is 
		
    û_kd^( t+1)   =   ∑_i ∈ Ad_i (1- p̂_i^(t))  γ_ik^(d)/∑_i ∈ A d_i (1-p̂_i^(t) ) .
  
		
  * Set t=t+1 and go to Step 2. Continue until convergence. 
	 


Once δ̂_i are computed, we may want to  use, instead of (<ref>),

    ∑_i ∈ A w_i ( 1, δ̂_i, δ̂_i y_i ) = ∑_i ∈ U (1, δ̂_i, δ̂_i y_i )
 
as the calibration equation.  The calibration estimator using (<ref>) is equivalent to 

    T̂_PDI2 = T_b2 + (N-N_b2) ∑_i ∈ A d_i (1-δ̂_i) y_i/∑_i ∈ A d_i (1-δ̂_i)   ,

where 
(N_b2,  T_b2)  = ∑_i ∈ Uδ̂_i(1,  y_i). However, unless δ̂_i = δ_i, we do not observe N_b2 and T_b2 and cannot compute T̂_PDI2 above. 

To overcome this difficulty, note that p̂_i in (<ref>)  is a consistent estimator of E( δ_i |_i). Thus, as long as 

    P( δ=1 |, y) = P( δ=1 |)
 
holds, we can  estimate N_b2 and T_b2 consistently  by applying the standard propensity score method using p̂_i. That is, use 

    ( N̂_b2, 
    T̂_b2)  = ∑_i ∈ Uδ_i/p̂_i δ̂_i(1,  y_i)    = ∑_i ∈ Bδ̂_i/p̂_i (1, y_i )
 
as a propensity score  estimator of (N_b2,  T_b2)  = ∑_i ∈ Uδ̂_i(1,  y_i). Unlike <cit.>, the estimated propensity scores p̂_i are fully nonparametric.  
Ignoring estimation errors in p̂_i, we have 

    E_δ{ ( N̂_b2,  T̂_b2 ) }≅ E_δ{∑_i ∈ Uδ_i/p_i δ̂_i(1,  y_i)    } =  ∑_i ∈ Uδ̂_i(1,  y_i)  = (N_b2, T_b2 ),

where E_δ ( ·) denotes the expectation with respect to δ and the first equality holds because E( δ_i |_i, y_i ) = p_i.  
Thus, the resulting data integration estimator is 

    T̂_PDI2 = T̂_b2  + (N-N̂_b2) ∑_i ∈ A d_i (1-δ̂_i) y_i/∑_i ∈ A d_i (1-δ̂_i)   .
 

The data integration estimator in (<ref>) can be viewed as a calibration estimator with calibration equation 

    ∑_i ∈ A w_i ( 1, δ̂_i, δ̂_i y_i ) = ∑_i ∈ U (1, δ̂_iδ_i  /p̂_i, δ̂_i δ_i y_i  /p̂_i ) = (N,  ∑_i ∈ Bδ̂_i/ p̂_i,  ∑_i ∈ Bδ̂_i y_i/ p̂_i ),
 
 which requires computing p̂_i and δ̂_i for every unit in sample B. 
Condition (<ref>) can be understood as  the ignorability condition of the sampling mechanism for sample B. This condition is not as  strong  as it might look at first. If y is categorical, one can always include y into  and apply the proposed classification method. In this case, condition (<ref>) is always satisfied. For continuous y, we may categorize y first and include it into . See the simulation study in Section 7.2 for an example.  




§ HANDLING MEASUREMENT ERRORS IN SURVEY DATA



We now consider the case the measurement errors exist in the survey data. For example, survey data is collected annually, and the big data is available monthly. In this case, if we are interested in estimating parameters on a  monthly basis, we can treat the observed values  in the latest year from the survey data as an inaccurate measurement for y_i. Thus, we observe (δ_i,  y_i^*) from sample A and observe y_i from sample B. In this case, we can use the measurement error model (<ref>) to obtain a design-model based estimator of T=∑_i=1^N y_i.


To estimate T under measurement errors in sample A and selection bias in sample B, we consider the following two-step approach:


  * [Step 1] 	 Using the measurement model, estimate the parameters in E( y_i | y_i^*)= m(y_i^*; ) and obtain mass imputation for sample A. That is, create ŷ_i=m( y_i^*; β̂) for all elements in sample A. If the measurement error model is (<ref>), then we can use
	 ŷ_i = β̂_1^-1 ( y_i^* - β̂_0 ),
	 where (β̂_0, β̂_1) is the estimated parameter from the elements in A ∩ B.

  * [Step 2]
 Apply calibration estimation  using _i= (1-δ_i, δ_i, δ_i y_i)^. That is, the final estimator is
 
    T̂_RegDI = ∑_i ∈ A w_i ŷ_i,

 where w_i minimizes Q( d, w) subject to the calibration equation ∑_i ∈ A w_i _i = ∑_i ∈ U_i.



In Step 1, the bias-corrected predictor of y_i  is obtained from model (<ref>). In principle, since we observe (y_i, y_i^*) among those with δ_i=1 in sample A, we can treat this sample, A ∩ B, as the validation sample in the calibration study. If the mechanism for δ_i=1 depends on y only, then the measurement error model (<ref>) is non-informative in the sense of <cit.>.  In this case, we can estimate model parameters in (<ref>) consistently by the complete-case analysis. 
That is, we can use 
    ∑_i ∈ A d_i δ_i ( y_i^* - β_0 - β_1 y_i) (1, y_i) = (0,0)
 as an estimating equation for (β_0, β_1). 


For variance estimation of  T̂_RegDI in  (<ref>), we can use, similarly to (<ref>), 

    V̂ ( T̂_RegDI) = ∑_i ∈ A∑_j∈ Aπ_ij - π_i π_j/π_ijê_i /π_iê_j /π_j,

where  ê_i = ŷ_i - _i^𝐁̂ and  
𝐁̂ =( ∑_i ∈ A d_i _i _i^)^-1∑_i ∈ A d_i _i ŷ_i. Thus, we can safely ignore the effect of uncertainty of   in ŷ_i= m( y_i^*; )  for  variance estimation. See  Appendix A for a sketched  justification.







§ SIMULATION STUDY



 §.§ Simulation study one


In the first simulation, continuous Y variable is considered from the following model:

    y_i  =   3+ 0.7    (x_i - 2) + e_i ,

where x_i ∼ N(2, 1), e_i ∼ N(0, 0.51),  and e_i is independent of x_i. We generate a finite population of size N=1,000,000 from this model. 
 Also, we generate

    y_i^* =  2 + 0.9 ( y_i-3) + u_i

where u_i ∼ N(0, 0.5^2), and u_i is independent of y_i.





In this simulation, we repeatedly obtain two samples, denoted by A and B, by simple random sampling of size n=1,000 and by an unequal probability sampling of size N_B=500,000, respectively.
In selecting sample B, we create two strata, where stratum 1 consists of elements with x_i≤ 2, and stratum 2 consists of those with x_i >2.
Within each stratum, we select n_h elements by simple random sampling independently, where n_1=300,000 and n_2=200,000. Under this sampling mechanism,
the sample mean of B is smaller than the population mean. We assume that the stratum information is  not available at the time of data analysis.



We consider the following three scenarios:
	

  * [Scenario I] No measurement errors in both samples. Thus, we observe y_i in both samples.

  * [Scenario II]
	Measurement errors in sample B. Thus, we observe y_i in sample A and  y_i^* in sample B.

  * [Scenario III] 	Measurement errors in sample A. Thus, we observe y_i^* in sample A and y_i in sample B.

In addition, assume that we observe the matching indicator δ_i in sample A. If δ_i=1 in sample A, we observe (y_i, y_i^*).

We consider the following four estimators for the population mean of Y:
	

  * Mean A. Mean of sample A observations.

  * Mean B. Mean of sample B observations.

  * Post-stratified data integration estimator of the form (<ref>).

  * Regression data integration estimator of the form (<ref>).







In Scenario II, the post-stratified data integration estimator is computed using

    θ̂_PDI =  1/N{∑_i=1^N δ_i y_i^*  + (N-N_b) ∑_i ∈ A d_i (1-δ_i) y_i/∑_i ∈ A d_i (1-δ_i) }.

In Scenario III, the post-stratified data integration estimator is computed using

    θ̂_PDI = 1/N{∑_i=1^N δ_i y_i  + (N-N_b) ∑_i ∈ A d_i (1-δ_i) y_i^*/∑_i ∈ A d_i (1-δ_i) },

and the regression data integration estimator is computed from the two-step approach in  (<ref>).





Table 2 presents the result of the  simulation study based on 1 000 Monte Carlo samples.

From  Table 2,  mean A estimator is unbiased except for Scenario III, where systematic measurement errors exist in sample A.  Mean B estimator is always biased due to the selection bias in sample B. The bias is the largest in absolute values for Scenario II, where measurement errors exist in addition to the  selection bias. Variance of mean B estimator is the smallest because of the large sample size of sample B (N_B=500,000). The post-stratified data integration estimator is unbiased in Scenario I, which is consistent with our theory in Section  3.
The variance of the post-stratified estimator  is about half of the variance of the mean A estimator because N_B/N=0.5. If the rate W_B=N_B/N is larger, then the variance estimator  post-stratified estimator will be smaller as equation (4) may suggest. 
However, in Scenario II, the post-stratified data integration  estimator is biased because T_b=∑_i=1^N δ_i y_i is estimated  without correcting for the measurement errors. In Scenario III, it is biased because T_c=∑_i=1^N (1-δ_i) y_i is estimated from sample A without correcting for the measurement errors.           The regression data integration estimator is unbiased for all scenarios. It is the same as the post-stratified data integration estimator under Scenario I, as   discussed in (<ref>).


In addition, we also compute variance estimators of the regression data integration estimator using formula (<ref>). For example, in Scenario 2, we use 
    
    ê_i = {[ y_i- (b̂_0 + b̂_1 y_i^*)                  δ_i = 1;               y_i - y̅_c                δ_i = 0 , ].

    where (b̂_0, b̂_1) is the solution to 
∑_i ∈ A d_i δ_i ( y_i - b_0 - b_1 y_i^* ) (1, y_i^*) = (0,0). Based on 1,000 Monte Carlo samples, we compute the relative biases of the variance estimators. The relative biases are -0.0037, 0.028, and 0.019 for Scenarios 1, 2, and 3, respectively. Thus, we  conclude that the proposed variance estimators are nearly unbiased.  



 §.§ Simulation study two


In the second simulation study, we study the performance of the data integration estimator using matching variables. 
 In the simulation study, we first generate a finite population with (z_i1, z_i2, δ_i, y_i) as follows. First generate 

    z_1i∼{1, ⋯, 20 }

independently. Given z_1i, we generate δ_i from Bernoulli distribution with the probability 

    P ( δ_i=1 | z_1i  )  
    = {[   c   z_i1≤ 10; 2 c  z_i1 > 10 ].

where c is chosen such that the sum of the probabilities over the finite population  is equal to N_B. We set N=10,000 and N_B=5,000 in this simulation. 
 We also generate 
 
    y_i = {[ 4+ 0.5 (z_i2+e_i)          z_1i≤ 10; 6+ 0.3 (z_i2+e_i)         z_1i > 10;                   ].

 where
z_2i∼{1, ⋯, 10}, e_i ∼(0,1) and z_2i and e_i are mutually independent. Thus, we can treat z_2i as a categorization of  continuous variable y_i. 

From the finite population, we select sample A by simple random sampling of size n_A. 
Two values of n_A=|A| are considered: n_A=1,000 versus n_A=2,000. From sample A, we observe (z_i1, z_i2, y_i) but not δ_i. Thus, we apply the semi-supervised classification method using (z_1i, z_i2) as the matching variable. Note that as z_i2 is included in the matching to satisfy the ignorability condition  (<ref>) approximately. 


From 
each sample, we consider five  estimators of Y̅_N =N^-1∑_i=1^N y_i. 
	

  * Mean A. Mean of sample A observations.

  * Mean B. Mean of sample B observations.

  * Naive data integration (DI)  estimator: Treat δ̂_i  as if accurate and apply the data integration estimator using δ̂_i to get 

    T̂_PDI = T_B + ( N- N_b ) ∑_i ∈ A d_i (1- δ̂_i) y_i/∑_i ∈ A   d_i (1- δ̂_i) .


  * The proposed data integration estimator: 

    T̂_PDI2 = T̂_b2  + (N-N̂_b2) ∑_i ∈ A d_i (1-δ̂_i) y_i/∑_i ∈ A d_i (1-δ̂_i)   ,

where T̂_b2 and N̂_b2 are defined in (<ref>). 



  * The original data integration estimator  using the true indicator function δ_i. This estimator is computed as a benchmark for comparison. 







Table <ref> presents the performance of the five estimators. Mean B is seriously biased due to its selection bias. Naive DI estimator is also biased seriously due to the misclassification errors in δ̂_i. The proposed DI estimator is unbiased and is more efficient than the sample mean (Mean A) of the sample A, although the efficiency gain is not as significant as in the original DI estimator. The efficiency gain will increase with W_b=N_b/N. 








§ AN APPLICATION IN OFFICIAL STATISTICS


We now consider an application of the proposed method to a real data problem  using  2015-16 Australian Agricultural Census as the big data, which has 85% response rate. In addition, we use the 2014-15 Rural Environment and Agricultural Commodities Survey (REACS) as the probability sample (sample A)  for calibration. Our interest is to combine the Agricultural Census data with the REACS data to estimate  the total area of holdings (AOH), the total number of dairy cattle (DAIRY), the number of beef cattle (BEEF), and the number of tonnes of wheat for grain or seed produced (WHEET) for 2015-16. Thus, we observe y_i from the Agricultural Census data and observe y_i^* from REACS.

To apply the proposed method, define δ_i=1 if unit i participated at the Census and δ_i=0 otherwise. Thus, in REACS sample, we observe y_i in addition to y_i^* for δ_i=1. Using the matched sample in sample A, we can fit a measurement error model

    y_i^* = β_0 + β_1 y_i + u_i

and obtain ŷ_i =   β̂_1^-1 ( y_i^* - β̂_0 ) for all i ∈ A. Here, y_i is the true value of the study variable 
from 2015-2016 Census 
and y_i^* is its proxy value obtained from  2014-2015 REAC survey data.  


For each parameter, we compute the following three estimators:


  * Survey estimate (from REACS sample): θ̂_HT = ∑_i ∈ A w_i ŷ_i

  * Big data estimate (from Census):  θ̂_B = ∑_i ∈ B y_i

  * The two-step Data integration estimate using calibration weighting:

    θ̂_DI = ∑_i ∈ A w_i, calŷ_i

where w_i, cal satisfies
∑_i ∈ Aw_i, cal ( 1- δ_i  , δ_i, δ_i x_i ) = ∑_i ∈ U ( 1- δ_i  , δ_i, δ_i x_i )
and x_i includes major study variables.

The estimates are compared with the official numbers of the Australian Bureau of Statistics (ABS), which is obtained by applying imputation for item nonresponse in the Census.


    <  >



    <  >



Figure 1 and Figure 2 present the estimation results for AOH and DAIRY, respectively,  by eight states in Australia.
We do not report the results for other commodities to save space. The confidence intervals are constructed using the asymptotic normality with 90% nominal coverage rates.  The results in Figure 1 and Figure 2 can be summarized as follows: (1) The Big data estimates show serious negative biases due to the undercoverage of the big data (nonresponse in the Census), (2) The proposed data integration estimator shows narrower confidence intervals than the survey estimate,  (3) The effect of calibration weighting is reduced because of the measurement errors in sample A observations.
 Overall, the confidence intervals obtained from the proposed data integration estimators cover the official ABS estimates.




§ DISCUSSION





The proposed data integration methods  feature an independent probability sample for estimating the missing data stratum of the finite population, which can correct for the under-coverage bias of the big  data sample.  By treating big data as an incomplete sampling frame for the finite population, we can apply the calibration weighting method.
 In addition,  these methods are extended in this paper to handle measurement errors in either the Big Data source or the probability sample source.   Also, a fully nonparametric approach to propensity score estimation for big data sample participation is developed using a new  semi-supervised classification method. 
 



In practice, our methods are useful provided the following conditions apply:


  * Existence of a probability sample A which also measures y or provides a proxy measure y^*. Whilst the coincidental existence of such a sample is rare, where one, e.g. a national statistical offices, determines the benefits in using big data for inference outweighs the costs, one can design, develop and implement such a random sample to collect the measure of interest.  Where this occurs, the population count of the sample units, N,  is by definition known.


        
  * The calibration method is useful only if the coverage of B is substantial, which is not an unreasonable assumption if B is a big data set.  Also, when B is big, it can be assumed that A ∩ B is not empty for measurement error  adjustment, where warranted;
       




The nonparametric propensity scores  obtained from the semi-supervised classification method can be used to correct for the  coverage bias in big data sample. How to make valid statistical inference, including variance estimation, under the nonparametric propensity score adjustment is not pursued here and will be covered elsewhere.  Extensions to  small domain estimation <cit.> and analytic inferences using big data  will also be  future research topics.



§ ACKNOWLEDGEMENTS


The authors are grateful  to two anonymous referees and the co-editor for the very constructive comments. 
The research of the first author was partially supported  by a grant from US National Science Foundation (MMS-1733572). 



 





§ APPENDIX



 §.§ A. Justification for (<ref>) 



Let θ=N^-1 T, the finite population mean of Y,  be the parameter of interest.
We first consider variance estimation of the mass imputation estimator of the form

    θ̂_DI = 1/N∑_i ∈ A d_i ŷ_i,

where ŷ_i is a predictor of y_i using y_i^*. We use ŷ_i = m̂^-1 ( y_i^*) where m̂ ( y_i )= m( y_i ; ) = E( y_i^* | y_i ; ). The estimating equation for  can be written as 

    Û_β ( ) = N^-1∑_i ∈ A d_i δ_i { y_i^* -   m( y_i ; ) }( y_i ; ) = 0

for some ( y; ) such that Û_β ( ) is linearly independent. 
Writing
θ̂_DI = θ̂_DI ( ),  we can use Taylor linearization   to estimate the variance of θ̂_DI.
Using the standard argument <cit.>,  we can obtain

    θ̂_DI =  θ̂_DI ( _N ) - E{∂/∂^θ̂_DI ( _N )  }[ E{∂/∂^Û_β ( _N )  }]^-1Û_β( _N ) + o_p (n^-1/2) ,

where   _N is the probability limit of . 

After some algebra,  we can express (<ref>) as

    θ̂_DI = 1/N∑_i ∈ A d_i { q_i  + δ_i ( y_i^* -m (y_i; )   )  κ^_i    } + o_p( n^-1/2 ) .

where q_i = q_i (_N) is the solution to y_i^*  = m(q_i ; _N), _i =  (y_i; _N) and κ satisfies 

    ∑_i=1^N δ_i ṁ_i _i^κ = ∑_i=1^N  q̇_i

with  ṁ_i = ∂ m(y_i; ) / ∂ and q̇_i = ∂ q_i( )/ ∂. 
Using (<ref>),  we can express

    θ̂_DI - θ =( q̅_N - θ ) + (u̅_HT - u̅_N)  + o_p (n^-1/2),

where q̅_N = N^-1∑_i=1^N q_i,  u_i =  q_i+ δ_i {  y_i^* -  m(y_i; _N )  }  (   κ^_i  ), u̅_HT = N^-1∑_i ∈ A d_i u_i and u̅_N = N^-1∑_i=1^N u_i.

From (<ref>), we can obtain

    ( θ̂_DI - θ)
    =  ( q̅_N - θ ) + (u̅_HT - u̅_N) = V_1 + V_2 .

The first term is of order O(N^-1), and the second term is O(n^-1). The first term is negligible if n/N=o(1).
To estimate the second term of (<ref>), we can use

    V̂_2 = 1/N^2∑_i ∈ A∑_j ∈ Aπ_ij - π_i π_j/π_ijû_i /π_i û_j /π_j ,

where

    û_i = q̂_i  + δ_i {  y_i^* - m ( y_i; )  } ( κ̂^ _i )

and

    κ̂ = {∑_i ∈ A d_i δ_i ṁ_i _i^}^-1∑_i ∈ A d_i   q̇_i .


Next, we consider variance estimation for the calibration estimator 

    θ̂_RegDI=1/N∑_i ∈ A w_i ŷ_i ,

where w_i are the calibration weights satisfying the calibration equation ∑_i ∈ A w_i _i = ∑_i =1^N _i.    In this case, the linearization in (<ref>) changes to 

    θ̂_RegDI = 1/N∑_i ∈ A d_i { e_i  + δ_i ( y_i^* -m (y_i; _N )   )  κ_2^_i    } + o_p( n^-1/2 ) ,

where e_i = q_i - _i^𝐁, 𝐁  = ( ∑_i =1^N _i _i^)^-1∑_i =1^N _i  q_i and 
κ_2 satisfies 

    ∑_i=1^N δ_i ṁ_i _i^κ_2 = ∑_i=1^N e_i .

Since _i includes an intercept term, we have ∑_i=1^N e_i =0, which implies κ_2=0. 
 Therefore, for variance estimation of θ̂_RegDI, we can use (<ref>), where
ê_i =q̂_i - _i^𝐁̂
and 𝐁̂ = ( ∑_i  ∈ A d_i _i _i^)^-1∑_i ∈ A d_i _i q̂_i . 



chicago








