Estimating General Parameters from Non-Probability Surveys Using Propensity Score Adjustment

mathematics Article Estimating General Parameters from Non-Probability Surveys Using Propensity Score Adjustment Luis Castro-Martín , María del Mar Rueda * and Ramón Ferri-García Department of Statistics and Operational Research, University of Granada, 18071 Granada, Spain; [email protected] (L.C.-M.); [email protected] (R.F.-G.) *Correspondence: [email protected] Received: 20 October 2020; Accepted: 21 November 2020; Published: 23 November 2020   Abstract: This study introduces a general framework on inference for a general parameter using nonprobability survey data when a probability sample with auxiliary variables, common to both samples, is available. The proposed framework covers parameters from inequality measures and distribution function estimates but the scope of the paper is broader. We develop a rigorous framework for general parameter estimation by solving survey weighted estimating equations which involve propensity score estimation for units in the non-probability sample. This development includes the expression of the variance estimator, as well as some alternatives which are discussed under the proposed framework. We carried a simulation study using data from a real-world survey, on which the application of the estimation methods showed the effectiveness of the proposed design-based inference on several general parameters. Keywords: nonprobability surveys; propensity score adjustment; survey sampling MSC: 62D05 1. Introduction Nonprobability samples are increasingly common in empirical sciences. The rise of online and smartphone surveys, along with the decrease of response rates in traditional survey modes, have contributed to the popularization of volunteer surveys where sampling is non-probabilistic. Moreover, the development of Big Data involves the analysis of large scale datasets whose obtention is conditioned by data availability and not by a probabilistic selection, and therefore they can be considered large nonprobability samples of a population [1]. The lack of a probability sampling scheme can be responsible for selection bias. Following the description from [ 1 , 2 ], we can distinguish the target population, UT , the subpopulation that a given selection method can potentially cover, Upc , and the fraction of the subpopulation that is finally covered, Uf c , and whose individuals might participate in the survey. Selection bias occurs when the characteristics of the individuals in Uf c differ significantly from those in UT in a way that could affect final estimates. Typically, differences between individuals in UT and individuals in Upc are caused by a lack of coverage induced by the survey administration mode (for example, an online questionnaire cannot be administered to the population without internet access), while differences between Upc and Uf c are caused by the variability in the propensities to participate between social-demographic groups (for example, an online questionnaire accesible in a thematic website might only be fulfilled by visitors of the website who have a specific interests that could influence the results). Following the rise of nonprobability samples, a class of methods for reducing selection bias have been proposed in the last decades. These methods were developed from different perspectives according to the availability of auxiliary information. We can mention calibration, Propensity Score Mathematics 2020,8, 2096; doi:10.3390/math8112096 www.mdpi.com/journal/mathematics Mathematics 2020,8, 2096 2 of 14 Adjustment (PSA), Statistical Matching and superpopulation modelling as the most relevant techniques to mitigate selection bias produced by coverage and self-selection errors. Calibration weighting was originally developed by [ 3 ] as a method to correct representation issues in samples with coverage or non-response errors. It only requires a vector of auxiliary variables available for each individual of the sample and the population totals of those variables. Calibration is able to remove selection bias in nonprobability samples if the selection mechanism is ignorable [ 4 ], and despite being originally developed for parametric estimation, further work [ 5 – 7 ] has extended calibration to distribution function, quantile and poverty measures estimation. Propensity Score Adjustment (PSA) and Statistical Matching require, apart from the nonprobability sample, a probability sample to do the adjustments. PSA was originally developed for balancing groups in non-randomized clinical trials [ 8 ] and it was adapted for non-response adjustments shortly after [ 9 , 10 ]. The application of PSA for removing bias in nonprobability surveys was theoretically developed in [ 11 , 12 ]. Statistical Matching was firstly proposed in [ 13 ] and extended in [ 14 ] for non-response adjustments. The difference between both methods is the sample used in the estimators: PSA estimates the propensity of each individual of the nonprobability sample to participate in the survey and then this propensity is used to construct the weights of the estimators, while Statistical Matching adjusts a prediction model using data from the nonprobability sample, applies it in the probability sample to predict their values for the target variable y and uses them in the parametric estimators. To the best of our knowledge, PSA and Statistical Matching has not been developed for nonparametric estimation. Superpopulation modelling requires data from the complete census of the target population for the covariates used in the adjustment, which is assumed to be a realization (sample) of a superpopulation where the (unknown) target values follow a model. It is based on the works by [ 15 , 16 ], where the main idea is to fit a regression model on the target variable with data from the nonprobability sample, and use the model to predict the values of the target variable for each individual in the population. The prediction can be used for estimation using a model-based approach or some alternative versions such as model-assisted and model-calibrated. LASSO models [ 17 ] and Machine Learning predictors [ 18 , 19 ] have been studied as alternatives to ordinary least squares regression in superpopulation modelling. The interest of society on poverty and inequality has increased in the last decades given the successive economic cycles and crisis. In such a context, official poverty rates and the percentage of people in poverty (or under a poverty threshold) are some important measures of a country’s wealth. The common characteristic of many poverty measures is their complexity. The literature on survey sampling is usually focused on the goal of estimating linear parameters. However, it is usual that the variable of interest in poverty studies is a measure of wages or income, where the distribution function becomes a relevant tool because it is required to calculate the proportion of people with low income, the poverty gap and other measures. Estimators for the cumulative distribution function, quantiles [ 20 , 21 ] and poverty measures [ 22 ] can be found in literature regarding probability samples, but there is hardly any work on the estimation of these parameters when the samples are obtained from volunteers. In this paper, we aim to develop a framework for statistical inference on a general parameter with non probability survey samples when a reference probability sample is available. After introducing the problem of the mean estimation for volunteer samples in Section 2, in Section 3, we consider the problem of the estimation for a general parameter through general estimating equations. Section 4 presents a new estimator for a general parameter through the use of PSA to estimate the propensity score of each individual in the survey weighted estimating equation and major theoretical results are presented. Results from simulation studies are reported in Sections 5and 6presents the concluding remarks. Mathematics 2020,8, 2096 3 of 14 2. Approaches to Estimation of a Mean for Volunteer Online Samples Let UT be the target population with N elements and sv a nonprobability sample drawn from a subset of UT , Uv , with a size of nv≤N . Let y be the target variable of the survey, whose mean in the population UT is denoted as Y . The sample estimation of Y , ˆ Y , is done using the Horvitz-Thompson estimator: ˆ YHT =∑i∈svwiyi ∑i∈svwi (1) where w is a vector of weights that accounts for the lack of representativity of sv caused by selection bias. If no auxiliary information is given, the weight would be the same for every unit, wi=N/nv , which requires to assume that the sample was drawn under a simple random sampling scheme. This is a naïve assumption given that sv is not probabilistic, this is, the probability of being in the sample is unknown and/or null for any of the units in UT. Let x be a matrix of covariates measured in sv along with y . If the population totals of the covariates, X , are available, it is possible to estimate the mean using a vector of weights obtained with calibration, wCAL . The calibration weights aim to minimize the distance between the original and the new weights min wCAL i E"∑ i∈sv G(wi,wCAL i)#(2) while respecting the calibration equations ∑ i∈sv wCAL ixi=X. (3) Some choices for the distance G( ., . ) were listed in [ 3 ], along with the resulting estimators. Calibration weighting for selection bias treatment was studied in [ 4 ], where post-stratification, which is a special case of calibration [ 23 ], was used to mitigate the bias caused by different selection mechanisms, showing its efficacy when the selection of the units of svis Missing At Random (MAR). If a reference sample, sr , drawn from the population UT is available and a number of covariates x have been measured both in sv and sr , two procedures can be done to reduce selection bias present in sv. Let Ivbe an indicator variable of an element being in sv, this is Ivi =(1i∈sv 0i/∈sv (4) Propensity Score Adjustment (PSA) assumes that each element of UT has a probability (propensity) of being selected for svwhich can be formulated as πv i=Pr(Ivi =1|xi,yi)(5) where πv i is the propensity of the i -th individual to participate in sv . The random mechanism behind this probability is the selection mechanism that governs the nonprobability sample. If the selection is Missing Completely At Random (MCAR), then πv i=Pr(Ivi = 1 ) and the selection bias is null, while if the selection is MAR then πv i=Pr(Ivi = 1 |xi) and the selection mechanism is considered ignorable. This does not mean that the selection bias should be ignored but rather it can be treated with the right techniques. In PSA, we consider the situation where true propensities are not known and therefore have to be estimated; we do it by combining sv and sr into a sample. The probability that Iv= 1 is then estimated using a prediction model, traditionally a logistic regression one: ˆ πv i=1 1+exp{−βxi}(6) Mathematics 2020,8, 2096 4 of 14 Alternative models, such as non-linear regression and Machine Learning classification algorithms, have been studied in literature as a substitute of logistic regression (see [ 24 ] for a review). The resulting propensities can be used to adjust new weights, wPSA, with different alternatives: •A simple inverse probability weighting is proposed by [25] wPSA1 i=wi ˆ πv i (7) which is a similar approach to the formula used in [26] wPSA2 i=1−ˆ πv i ˆ πv i (8) • Alternatively, individuals of the combined sample ( sv∪sr ) can be grouped in g equally-sized strata of similar propensity scores from which an average propensity is calculated for each group. Let πg be the mean propensities of the g -th strata. [ 2 ] use the means as in (7) to calculate the new weights: wPSA3 i=wi πgi (9) where girefers to the strata to which the i-th individual of svbelongs. • A similar approach can be found in [ 12 ], but instead of using the means, a factor is calculated for each strata: fg= ∑k∈srg˜ wk/∑k∈sr˜ wk ∑i∈svgwi/∑i∈svwi (10) where srg and svg are respectively the individuals from the probability and nonprobability sample that belong to the g -th strata, and ˜ w is the vector of design weigths of the reference sample. The final weights are obtained by multiplying the original weights and the correction factor: wPSA4 i=wi·fgi(11) PSA has been proven to successfully remove selection bias when prognostic covariates are chosen [ 11 ] and further adjustments, such as calibration, are applied in the estimations [ 2 , 12 , 27 ]. A recent paper [ 28 ] shows a real application of PSA in web panel surveys where the reductions in bias, although present, were not large enough to consider the estimates as unbiased. As an alternative to PSA, Statistical Matching is another method to mitigate selection bias when a reference sample is available. For the matter, a prediction model for y using x as the dependent variables is built using data from sv . The model is subsequently applied on the reference sample to obtain the estimates from the predicted values of yin sr,ˆ y: ˆ Y=∑ k∈sr wkˆ yk(12) The choice of prediction models has been studied in literature; the usual method is linear regression but other approaches such as donor imputation [ 13 ] or Machine Learning algorithms [ 19 , 29 ] have been listed as alternatives. Under certain conditions, Statistical Matching can reduce bias and mean square error to a greater extent than PSA [29]. When a complete census of the entire target population is available, with information on the covariates present in sv , superpopulation modelling can be applied to remove selection bias [ 19 ]. In this paper we consider the case when auxiliary information is available only from a reference probability survey. Mathematics 2020,8, 2096 5 of 14 3. Estimation of a General Parameter by Using PSA Let y be the variable of interest in a survey and yi be the value of the i -th unit in that variable, i= 1, . . . , N . Suppose we want to estimate a finite population parameter θN of dimension p≥ 1 defined as the solution of the census estimating equations: U(θN) = 1 N∑ U ui(yi,θN) = 0 (13) where ui(yi,θN)is be a function of θN. Some unidimensional parameters of interest can be: •the population total Tyfor ui= (yi−θN/N), •the population mean ¯ Yfor ui= (yi−θN), • the population distribution function Fy(t) for ui= ( 1 (yi≤t)−θN) with 1 (·) being the indicator function, •the finite population quantile of order j,Qjfor ui= (1(yi≤θN)−j, where 0 <j<1, We denote by ˆ θthe solution of the equation: ˆ U(θN) = ∑ U Iviui(yi,θN)/πv i=∑ sv ui(yi,θN)/πv i=0. (14) It is clear the Er(ˆ U(θN)) = U(θN) where r stands for the model of the selection mechanism for the sample sv , this is, the true model that fits propensity scores. If πv i are known we can get the consistent estimator of θN by solving the equation above. For the study of the properties of this estimator we consider a quasi-probability approach or pseudo-design-based inference ([ 19 ]) and we treat the volunteer sample as a realization of a Poisson sampling with probabilities πv i. For any sample design that verifies certain regularity conditions, the solution to b U(θ) = 0 provides a consistent estimator for the parameter θN (see [ 30 ]). Poisson sampling verifies these conditions, so that the consistency of the estimator is obtained immediately from the result of [ 30 ]. The normality of the estimator is demonstrated by [ 31 ], who also obtains the asymptotic variance of the estimator. From said expression and taking into account that in Poisson sampling the extractions are independent and therefore the probability of second order is given by πv ij =πv iπv jwe can obtain the variance of ˆ θ: V(ˆ θ) = J(ˆ θ)−1var(ˆ U(θ))J0(ˆ θ)−1(15) being J(θ) = 1 N∑U∂ui/∂θ and var(ˆ U(θ)) = ∑U(1−πv i)u2 i/(πv i)2 4. Estimation of a General Parameter with Estimated Propensities The propensity scores πv i are not known are impossible to estimate using the nonprobability sample sv alone, so additional information must be included. Let sr be a reference probability sample, of size nr , selected from UT under a sampling design (sd , pd) where the first order inclusion probabilities, πp i=∑sr3ipd(sr),i=1, . . . , nr, are known and non-null. The covariates of the propensity model x have been measured both in sv and sr , while the variable of interest yis only available for those individuals in sv. Suppose that the propensity scores can be modelled parametrically as πv i=P(Ivi =1/xi) = m(λo,xi)i=1, . . . , N(16) for some known function m(·) with second continuous derivatives with respect to an unknown parameter λo. We estimate the propensity scores by using data of both the volunteer and the probability sample. The maximum likelihood estimator (MLE) of πv i is m(ˆ λ , xi) where ˆ λ corresponds to the value of lambda that maximizes the log-likelihood function: Mathematics 2020,8, 2096 6 of 14 l(λ) = ∑ U (Ivilog(m(λ,xi)) + (1−Ivi)log(1−m(λ,xi)) = ∑ sv log m(λ,xi) 1−m(λ,xi)+∑ U log(1−m(λ,xi)). (17) As it is usual in survey sampling, we consider the pseudo-likelihood given that some units of the population have not been sampled: ˜ l(λ) = ∑ sv log m(λ,xi) 1−m(λ,xi)+∑ sp 1 πp i log(1−m(λ,xi)). (18) We propose thus a two phase procedure in this manner: Step 1: Calculate ˆ λpl by solving the score equations: ∂˜ l(xi,λ)/∂λ =0 Step 2: Calculate ˆ θvas the solution of the estimating function: ˆ UV(θ) = ∑ U Iviui(yi,θ)1 m(ˆ λpl,xi)=0 (19) We consider the following asymptotic framework for theoretical development, which is equivalent to the framework in [ 32 ]. Let UTν be a sequence of finite populations of size Nν . Each UTν has an associated non-probability sample svν of size nvν and an associated probability sample spν of size npν . We consider that the population size Nν→∞ , the nonprobability sample size nvν→∞ and the probability sample size npν→∞ as ν→∞ . For notational simplicity the index ν is suppressed for the rest of the paper. The properties of the estimator ˆ θv are developed under both the model for the propensity scores and the survey design for the probability sample. We make the following assumptions: •A.1. The estimating function ui(yi,θ,λ)is twice differentiable with respect to θand λ. • A.2. The propensities and the sampling design ensure that ˆ UV(θ)−U(θ) = Op(n−1/2) for any θ∈Θ. • A.3. The propensities and the sampling design ensure that ˆ UV(θ) is asymptotically Normal with mean U(θ)and entries of the variance at the order O(n−1)for any fixed θ∈Θ. Theorem 1. Under the conditions A.1, A.2 and A.3, ˆ θv is a consistent and asymptotically normal estimator for θ . Proof. Under assumed conditions, ˆ UV(θ) = U(θ) + Op(n−1/2) , thus by using the mean value theorem, ˆ θv has the same asymptotic behaviour that ˆ θwhich is consistent for θand asymptotically normal distributed (see Section 3). Variance estimation for ˆ θv can be handled by combining the two estimating equations, ˜ l and ˆ Uv , into a single system as it is done in [33]. The MLE of λ,ˆ λpl is the solution to the equations: U2(λ) = ∑ sv ∂log m(λ,xi) 1−m(λ,xi)/∂λ +∑ sp ∂1 πp i log(1−m(λ,xi))/∂λ =0 Mathematics 2020,8, 2096 7 of 14 and the PSA estimator of θNis the solution to the estimating equations U1(θ,λ) = ∑ sv ui(yi,θ)1 m(λpl,xi)=∑ sv g1(yi,xi,θ,λ) = 0. Let U(θ , λ) = (U0 1(θ , λ) , U0 2(λ))0 . Let ψ= (θ0 N , λo)0 be the true parameter values defined through the census estimating equations and ˆ ψ= ( ˆ θ0 N,ˆ λo 0)0the solutions to U(θ,λ) = 0. We need an additional assumption: • A.4. The propensities, the sampling design and the estimating function satisfy ∂b U/∂ψ=Op( 1 ) and ∂2b U/∂ψ∂ψ0=Op(1). Theorem 2. Under the conditions A.1, A.2, A.3 and A.4, the asymptotic variance-covariance matrix of ˆ ψ is given by the expression: V(ˆ ψ) = H−1V(ˆ U(θ,λ))H0−1(20) with H= H11 H12 0H22! H11 =E{∂ ∂θ U1(θN,λ)} H21 =E{∂ ∂λU1(θN,λ)} H22 =E{∂ ∂λU2(λ)} Proof. Since ˆ θv and ˆ λ are consistent estimator of respective parameters, we can write ˆ ψ=ψ+Op( 1 ) and the Taylor series expansion gives: ˆ ψ=ψ−H−1b U(θ,λ) + Op(kˆ ψ−ψk2), Thus the asymptotic variance of ˆ ψis given by: V(ˆ ψ) = H−1V(b U(θ,λ))H0−1. Taking into account the two random mechanisms, and the probabilities of the conditional expectation, we have V(b U(θ , λ)) = VpEr(b U(θ , λ)) + EpVr(b U(θ , λ)) where r stands for the model of the selection mechanism for the sample sv and p refers to the probability sampling design for sp . The asymptotic variance of ˆ ψ depends on the probability of selecting the sample sp under the given sampling design and the selection mechanism described by the propensity model. Plug-in estimators can be used to construct variance estimators for all the required components but it is not a simple issue. In practice, and as described in [ 7 ], the use of jackknife [ 34 ] and bootstrap techniques [ 35 ] in the variance estimation for nonlinear parameters should be more advantageous because of their wide applicability for different cases and conditions. Direct applications of bootstrap methods for estimating the variance-covariance matrix of ˆ ψ involve solving the equation U(θ , λ) = 0 repeatedly for each bootstrap sample. Multiplier Bootstrap with Estimating Functions was proposed by [36]. 5. Simulation Study 5.1. Data Data for the simulation study come from a wave of the Spanish Living Conditions Survey collected between 2011 and 2012 [ 37 ], which contains an annual thematic module that, in 2012, was dedicated to Mathematics 2020,8, 2096 8 of 14 household conditions. The survey sampling follows a two-phase cluster sampling, where the primary units are the households and the secondary units are their members. In 2012, the final sample included 33,573 individuals. For this study, the dataset was filtered to rule out individuals and variables with high quantities of missing data. After this procedure, the dataset employed as pseudopopulation of the study had a size of N= 28,210 individuals and p= 60 available variables. From this pseudopopulation, two probability samples of size nr were drawn according to the following sampling strategies: • The first sample, sr1 , was drawn with a stratified cluster sampling, where the strata were defined by the Autonomous Communities (NUTS2 regions) and the clusters were the households, which were drawn with probabilities proportional to the household size. The number of households to be selected, m , was estimated dividing nr by the medium household size in order to reach the aforementioned size of nr= 2000, resulting in m= 902 households. The final sample size of sr1 was nr1=2003. • The second sample, sr2 , was drawn with an unequal probability sampling, where probabilities were proportional to the minimum income of the individual’s household to make ends meet (variable HS130 in [37]). The extraction of the nonprobability sample, sv , was done with unequal probability sampling from the full pseudopopulation, where the probability of selection for the i -th individual, pi , was given by the formula: pi=1 1+exp(−2x1 i+0.2x2 i+0.01x3 i+0.2x41 i+0.4x42 i)(21) where •x1 i=1 when the i-th sampled individual has a computer at home, and x1 i=0 otherwise. •x2 i=1 when the i-th sampled individual is a man, and x2 i=0 otherwise. •x3 iis the age (in years) of the i-th sampled individual. •x41 i= 1 when the i -th sampled individual lives in a medium population density area, and x41 i=0 otherwise. •x42 i= 1 when the i -th sampled individual lives in a low population density area, and x42 i= 0 otherwise. The reasoning behind this sampling procedure is to take into account more similar mechanisms to self-selection procedures that take place in real nonprobability surveys. We have considered three different sample sizes, nv = 2000, 4000, 6000. 1000 simulation runs were performed for each procedure and sample size, drawing a sample in each run. 5.2. Simulation In each simulation, the parameters to be estimated were the following: •The Gini coefficient [38], which measures the income inequality, estimated as b Gy=∑k∈sv 1 πk(2b Fy(yk)−1)yk ∑k∈svyk/πk • The proportion of individuals with a disposable income below the at-risk-at-poverty threshold. This measure can be referred to as poverty incidence, poverty proportion, poverty risk or HCI ([ 39 ] and is estimated as ˆ HCI =1 N∑ k∈sv 1 πk I(y<0.6Q0.5) •The interquartile range, estimated as d IQR =b Q0.75 b Q0.25 Mathematics 2020,8, 2096 9 of 14 •The interdecile range, estimated as d IDR =b Q0.9 b Q0.1 . Every parameter was estimated with and without applying PSA so we could evaluate its performance. In order to estimate the propensities, a logistic regression model was chosen: m(ˆ λ,xi) = exp(ˆ λTxi) 1+exp(ˆ λTxi) 1000 simulations were executed for each context. The resulting mean bias, standard deviation and Root Mean Square Error were measured in relative numbers to make them comparable across different scenarios. The formulas used for their calculation can be found below: RBias (%) =  ∑1000 i=1ˆ θ(i) 1000 −θN·100 θN (22) RStandard deviation (%) = s∑1000 i=1(ˆ θ(i)−ˆ ¯ θ)2 999 ·100 θN (23) RMSE (%) = pRBias2+RSD2(24) with ˆ θ(i)the estimation in the i-th simulation and ˆ ¯ θthe mean of the 1000 estimations. 5.3. Results The relative mean bias of the estimations can be observed in Tables 1–3. We can observe that PSA reduces the bias in all situations, specially in the estimation of HCI. PSA using the reference sample drawn with probabilities proportional to the income, sr2 , provided much less biased estimates overall. Table 1. Relative mean bias (%) of each parameter without applying PSA. Size Gini HCI IQR IDR 6000 6.7 80.4 7.9 9.4 2000 3.2 93 3.8 3.1 4000 3.1 86 3.7 3 6000 3 79 3.5 3 Table 2. Relative mean bias (%) of each parameter applying PSA with the stratified reference sample. Size Gini HCI IQR IDR 2000 1.7 3 2.5 1 4000 2.1 3.3 2.7 1 6000 2.2 3.1 2.7 0.9 Table 3. Relative mean bias (%) of each parameter applying PSA with the proportional reference sample. Size Gini HCI IQR IDR 2000 0.3 1.1 0 0.2 4000 0.1 1.3 0.1 0.3 6000 0 1.1 0.1 0.5 The relative standard deviation of the estimations can be observed in Tables 4–6. The standard deviation remained stable across estimates of Gini coefficient, IQR and IDR, even with small gains for