Doubly Robust Inference when Combining Probability and Non-probability Samples with High-dimensional Data

TL;DR

This paper develops a method for integrating high-dimensional probability and non-probability samples to make unbiased population inferences, using variable selection and doubly robust estimation techniques.

Contribution

It introduces a two-step approach combining penalized variable selection with a doubly robust estimator for finite population means in high-dimensional settings.

Findings

The method achieves selection consistency for general probability samples.

The doubly robust estimator is root-n consistent if either the sampling or outcome model is correct.

The approach effectively mitigates selection bias in high-dimensional survey data.

Abstract

Non-probability samples become increasingly popular in survey statistics but may suffer from selection biases that limit the generalizability of results to the target population. We consider integrating a non-probability sample with a probability sample which provides high-dimensional representative covariate information of the target population. We propose a two-step approach for variable selection and finite population inference. In the first step, we use penalized estimating equations with folded-concave penalties to select important variables for the sampling score of selection into the non-probability sample and the outcome model. We show that the penalized estimating equation approach enjoys the selection consistency property for general probability samples. The major technical hurdle is due to the possible dependence of the sample under the finite population framework. To overcome this challenge, we construct martingales which enable us to apply Bernstein concentration inequality for martingales. In the second step, we focus on a doubly robust estimator of the finite population mean and re-estimate the nuisance model parameters by minimizing the asymptotic squared bias of the doubly robust estimator. This estimating strategy mitigates the possible first-step selection error and renders the doubly robust estimator root-n consistent if either the sampling probability or the outcome model is correctly specified.

Tables

Table 1. Table 1: Two data sources. “ √ square-root \surd ” and “?” indicate observed and unobserved data, respectively.

		Sample weight ${\pi }^{-1}$	Covariate $X$	Study Variable $Y$
Probability	$1$	$\surd$	$\surd$	?
Sample	$⋮$	$⋮$	$⋮$	$⋮$
${𝒪}_A$	$n_A$	$\surd$	$\surd$	?
Non-probability	$n_A+1$	?	$\surd$	$\surd$
Sample	$⋮$	$⋮$	$⋮$	$⋮$
${𝒪}_B$	$n_A+n_B$	?	$\surd$	$\surd$

Table 2. Table 2: Simulation results for selection performance for the proposed double penalized estimating equation procedure under four scenarios: under OM I (II), the outcome model is correctly specified (misspecified), and under PSM I (II), the probability of sampling score model is correctly specified (misspecified)

	${\beta }^{\ast }$					${\alpha }^{\ast }$
	Under	Over	FN	FP		Under	Over	FN	FP
	$(\times {10}^2)$	$(\times {10}^2)$				$(\times {10}^2)$	$(\times {10}^2)$
Continuous outcome
(i) OM I and PSM I	0.0	31.8	0.0	1.4		0.0	0.0	0.0	0.0
(ii) OM II and PSM I	70.6	15.0	0.9	0.2		0.0	0.0	0.0	0.0
(iii) OM I and PSM II	0.0	32.8	0.0	1.4		100.0	100.0	4.0	1.0
(iv) OM II and PSM II	0.0	0.4	0.0	0.4		100.0	100.0	3.5	4.3
Binary outcome
(i) OM I and PSM I	0.0	0.0	0.0	0.0		0.0	0.0	0.0	0.0
(ii) OM II and PSM I	100.0	0.0	2.1	0.0		0.0	0.0	0.0	0.0
(iii) OM I and PSM II	0.0	0.0	0.0	0.0		100.0	100.0	4.0	1.0
(iv) OM II and PSM II	100.0	0.0	4.0	0.0		100.0	96.0	4.0	1.0

Table 3. Table 3: Simulation results for the coverage properties for the continuous and binary outcomes: empirical coverage rate and (empirical coverage rate ± 2 × \pm 2\times Monte Carlo standard error)

	Continuous outcome			Binary outcome
(i) OM I and PSM I		$95.2(93.3,97.1)$			$95.7(93.9,97.6)$
(ii) OM II and PSM I		$94.6(92.6,96.6)$			$95.5(93.6,97.4)$
(iii) OM I and PSM II		$96.2(94.2,97.8)$			$95.6(93.8,97.5)$
(iv) OM II and PSM II		$88.2(85.3,91.1)$			$42.9(38.3,47.6)$

Table 4. Table 4: Point estimate, standard error and 95 % percent 95 95\% Wald confidence interval

	$Y_1$ (days had at least one drink last month)				$Y_2$ (whether voted local elections)
	Est	SE	CI		Est$\times {10}^2$	SE$\times {10}^2$	CI$\times {10}^2$
Naive	5.36	0.90	(5.17,5.54)		75.3	0.5	(74.4,76.3)
Proposed method	4.84	0.15	(4.81,4.87)		71.8	0.2	(71.3,72.2)

Equations

$${\hat{\mu }}_{\mathrm{IPW}}={\hat{\mu }}_{\mathrm{IPW}}(\hat{\alpha })=\frac{1}{N}\sum _{i=1}^N\frac{I_{B,i}}{{\pi }_B(X_i^{\mathrm{T}}\hat{\alpha })}{Y}_i.$$

$$\sum _{i=1}^N\left\{\frac{I_{B,i}}{\pi (X_i^{\mathrm{T}}\alpha )}-\frac{I_{A,i}}{{\pi }_{A,i}}\right\}h(X_i;\alpha )=0,$$

$${\hat{\mu }}_{\mathrm{reg}}={\hat{\mu }}_{\mathrm{reg}}(\hat{\beta })=\frac{1}{N}\sum _{i=1}^N{I}_{A,i}{d}_{A,i}m(X_i^{\mathrm{T}}\hat{\beta }),$$

$${\hat{\mu }}_{\mathrm{cal}}={\hat{\mu }}_{\mathrm{cal}}=\frac{1}{N}\sum _{i=1}^N{\omega }_i{I}_{B,i}{Y}_i,$$

$${\hat{\mu }}_{\mathrm{dr}}={\hat{\mu }}_{\mathrm{dr}}(\hat{\alpha },\hat{\beta })=\frac{1}{N}\sum _{i=1}^N\left[\frac{I_{B,i}}{{\hat{\pi }}_B(X_i^{\mathrm{T}}\hat{\alpha })}\{Y_i-m(X_i;\hat{\beta })\}+I_{A,i}{d}_{A,i}m(X_i;\hat{\beta })\right].$$

$${\alpha }^{\ast }=\arg \underset{\alpha \in {\mathbb{R}}^p}{\min }E\left[{\pi }_B(X)\log \frac{{\pi }_B(X)}{{\pi }_B(X^{\mathrm{T}}\alpha )}+\{1-{\pi }_B(X)\}\log \frac{1-{\pi }_B(X)}{1-{\pi }_B(X^{\mathrm{T}}\alpha )}\right],$$

$$U_1(\alpha )=\frac{1}{N}\sum _{i=1}^N\left\{\frac{I_{B,i}}{{\pi }_B(X_i^{\mathrm{T}}\alpha )}-\frac{I_{A,i}}{{\pi }_{A,i}}\right\}{X}_i.$$

$$U_2(\beta )=\frac{1}{N}\sum _{i=1}^N{I}_{B,i}\left\{Y_i-m(X_i^{\mathrm{T}}\beta )\right\}{X}_i.$$

$$U^{\mathrm{p}}(\alpha,\beta)=U(\alpha,\beta)-\left(\begin{array}[]{c}q_{\lambda_{\alpha}}(|\alpha|)\mathrm{sign}(\alpha)\\ q_{\lambda_{\beta}}(|\beta|)\mathrm{sign}(\beta)\end{array}\right),$$

$$q_{\lambda }(|\theta |)=\lambda \left\{I(|\theta |<\lambda )+\frac{(a\lambda -|\theta |{)}_{+}}{(a-1)\lambda }I(|\theta |\ge \lambda )\right\},$$

$$\begin{gathered} \text{a.bias}({\alpha }^{\ast },{\beta }^{\ast })=E\left\{{\hat{\mu }}_{\mathrm{dr}}({\alpha }^{\ast },{\beta }^{\ast })-\mu \right\} \\ =E\left[\frac{1}{N}\sum _{i=1}^N\left\{\frac{I_{B,i}}{{\pi }_B(X_i^{\mathrm{T}}{\alpha }^{\ast })}-1\right\}\{Y_i-m(X_i^{\mathrm{T}}{\beta }^{\ast })\}\right]+E\left\{\frac{1}{N}\sum _{i=1}^N\left(I_{A,i}{d}_{A,i}-1\right)m(X_i^{\mathrm{T}}{\beta }^{\ast })\right\}. \end{gathered}$$

$$\frac{\partial\{\text{a.bias}(\alpha,\beta)\}^{2}}{\partial(\alpha_{\mathcal{C}}^{\mathrm{\scriptscriptstyle T}},\beta_{\mathcal{C}}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}}=2\times\text{a.bias}(\alpha,\beta)\times\left(\begin{array}[]{c}I_{B}\left\{\frac{1}{\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)}-1\right\}\{Y-m(X^{\mathrm{\scriptscriptstyle T}}\beta)\}X_{\mathcal{C}}\\ \left\{\frac{I_{B}}{\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)}-d_{A}I_{A}\right\}\partial m(X^{\mathrm{\scriptscriptstyle T}}\beta)/\partial\beta_{\mathcal{C}}\end{array}\right)$$

$$J(\alpha,\beta)=\left(\begin{array}[]{c}J_{1}(\alpha,\beta)\\ J_{2}(\alpha,\beta)\end{array}\right)=\left(\begin{array}[]{c}\frac{1}{N}\sum_{i=1}^{N}I_{B,i}\left\{\frac{1}{\pi_{B}(X_{i}^{\mathrm{\scriptscriptstyle T}}\alpha)}-1\right\}\{Y_{i}-m(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)\}X_{i\mathcal{C}}\\ \frac{1}{N}\sum_{i=1}^{N}\left\{\frac{I_{B,i}}{\pi_{B}(X_{i}^{\mathrm{\scriptscriptstyle T}}\alpha)}-d_{A,i}I_{A,i}\right\}\partial m(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)/\partial\beta_{\mathcal{C}}\end{array}\right)$$

$${\hat{\mu }}_{\text{p-dr}}={\hat{\mu }}_{\text{p-dr}}(\hat{\alpha },\hat{\beta })=\frac{1}{N}\sum _{i=1}^N\left\{I_{B,i}\frac{Y_i-m(X_i^{\mathrm{T}}\hat{\beta })}{{\pi }_B(X_i^{\mathrm{T}}\hat{\alpha })}+I_{A,i}{d}_{A,i}m(X_i^{\mathrm{T}}\hat{\beta })\right\},$$

$$U^{\mathrm{p}}(\widetilde{\theta})=U(\widetilde{\theta})-\left(\begin{array}[]{c}q_{\lambda_{\widetilde{\alpha}}}(|\widetilde{\alpha}|)\mathrm{sign}(\widetilde{\alpha})\frac{|\widetilde{\alpha}|}{\epsilon+|\widetilde{\alpha}|}\\ q_{\lambda_{\widetilde{\beta}}}(|\widetilde{\beta}|)\mathrm{sign}(\widetilde{\beta})\frac{|\widetilde{\beta}|}{\epsilon+|\widetilde{\beta}|}\end{array}\right)=0,$$

$$\nabla (\theta )=\partial U(\theta )/\partial {\theta }^{\mathrm{T}}$$

$$\frac{\partial U_1(\alpha )}{\partial {\alpha }^{\mathrm{T}}}$$

$$\frac{\partial U_2(\beta )}{\partial {\beta }^{\mathrm{T}}}$$

$$E(\theta )$$

$${\tilde{\theta }}_j^{[k]}={\tilde{\theta }}_j^{[k-1]}+{\left\{{\nabla }_{jj}({\tilde{\theta }}^{[k-1]})+N\cdot E_{jj}({\tilde{\theta }}^{[k-1]})\right\}}^{-1}\left\{U_j({\tilde{\theta }}^{[k-1]})-N\cdot E_{jj}({\tilde{\theta }}^{[k-1]}){\tilde{\theta }}_j^{[k-1]}\right\},$$

$$\mathrm{Loss}({\lambda }_{\alpha })=\sum _{j=1}^p{\left[\sum _{i=1}^N\left\{\frac{I_{B,i}}{{\pi }_B\{X_i^{\mathrm{T}}\hat{\alpha }({\lambda }_{\alpha })\}}-\frac{I_{A,i}}{{\pi }_{A,i}}\right\}{X}_{i,j}\right]}^2,$$

$$\mathrm{Loss}({\lambda }_{\beta })=\sum _{i=1}^N{I}_{B,i}{\left[Y_i-m\{X_i^{\mathrm{T}}\hat{\beta }({\lambda }_{\beta })\}\right]}^2,$$

$${\mathcal{M}}_{\alpha }=\{1\le j\le p:{\alpha }_j^{\ast }\ne 0\},{\mathcal{M}}_{\beta }=\{1\le j\le p:{\beta }_j^{\ast }\ne 0\},{\mathcal{M}}_{\theta }={\mathcal{M}}_{\alpha }\cup \{p+{\mathcal{M}}_{\beta }\}.$$

$$0<c_1\le {\lambda }_{\min }\left(\frac{1}{N}\sum _{i=1}^N{X}_i^{\mathrm{T}}{X}_i\right)\le {\lambda }_{\max }\left(\frac{1}{N}\sum _{i=1}^N{X}_i^{\mathrm{T}}{X}_i\right)\le c_2<\infty ,$$

$$P(|U_j^{\mathrm{p}}(\tilde{\theta })|=0,j\in {\mathcal{M}}_{\theta })$$

$$P\left(|U_j^{\mathrm{p}}(\tilde{\theta })|\le \frac{{\lambda }_{\theta }}{\log n},j\in {\mathcal{M}}_{\theta }^c\right)$$

$$P\left({\tilde{\theta }}_{{\mathcal{M}}_{\theta }^c}=0\right)$$

$${\tilde{\theta }}_{{\mathcal{M}}_{\theta }}-{\theta }_{{\mathcal{M}}_{\theta }}^{\ast }=O_P(\sqrt{s_{\theta }/n}),$$

$$\{(\hat{\alpha }-{\alpha }^{\ast }{)}^{\mathrm{T}},(\hat{\beta }-{\beta }^{\ast }{)}^{\mathrm{T}}\}=O_p(\sqrt{s_{\theta }/n}).$$

$$n^{1/2}\left\{{\hat{\mu }}_{\text{p-dr}}(\hat{\alpha },\hat{\beta })-\mu \right\}$$

$$=$$

$$n^{1/2}E[\left\{{\hat{\mu }}_{\text{p-dr}}({\alpha }^{\ast },{\beta }^{\ast })-\mu \right\}]=\frac{n^{1/2}}{N}\sum _{i=1}^{N}E\left\{\frac{I_{B,i}}{{\pi }_B(X_i^{\mathrm{T}}{\alpha }^{\ast })}-1\mid X_i\right\}E\left\{Y_i-m(X_i^{\mathrm{T}}{\beta }^{\ast })\mid X_i\right\}.$$

Full text

Doubly Robust Inference when Combining Probability and Non-probability Samples with High-dimensional Data

Shu Yang1 , Jae Kwang Kim2, and Rui Song1

1 Department of Statistics, North Carolina State University,

2Department of Statistics, Iowa State University syang24@ncsu.edu, Department of Statistics, North Carolina 27695, U.S.A.

Summary. We consider integrating a non-probability sample with a probability sample which provides high-dimensional representative covariate information of the target population. We propose a two-step approach for variable selection and finite population inference. In the first step, we use penalized estimating equations with folded-concave penalties to select important variables and show the selection consistency for general samples. In the second step, we focus on a doubly robust estimator of the finite population mean and re-estimate the nuisance model parameters by minimizing the asymptotic squared bias of the doubly robust estimator. This estimating strategy mitigates the possible first-step selection error and renders the doubly robust estimator root-$n$ consistent if either the sampling probability or the outcome model is correctly specified.

Keywords: Data integration; Double robustness; Generalizability; Penalized estimating equation; Variable selection

1 Introduction

Probability sampling is regarded as the gold-standard in survey statistics for finite population inference. Fundamentally, probability samples are selected under known sampling designs and therefore are representative of the target population. However, many practical challenges arise in collecting and analyzing probability sample data such as cost, time duration, and increasing non-response rates (Keiding and Louis, 2016). As the advancement of technology, non-probability samples become increasingly available for research purposes, such as remote sensing data, web-based volunteer samples, etc. Although non-probability samples do not contain information on the sampling mechanism, they provide rich information about the target population and can be potentially helpful for finite population inference. These complementary features of probability samples and non-probability samples raise the question of whether it is possible to develop data integration methods that leverage the advantages of both data sources.

Existing methods for data integration can be categorized into three types. The first type is the so-called propensity score adjustment (Rosenbaum and Rubin, 1983). In this approach, the probability of a unit being selected into the non-probability sample, which is referred to as the propensity or sampling score, is modeled and estimated for all units in the non-probability sample. The subsequent adjustments, such as propensity score weighting or stratification, can then be used to adjust for selection biases; see, e.g., Lee and Valliant (2009), Valliant and Dever (2011), Elliott and Valliant (2017) and Chen, Li and Wu (2018). Stuart et al. (2011; 2015) and Buchanan et al. (2018) use propensity score weighting to generalize results from randomized trials to a target population. O’Muircheartaigh and Hedges (2014) propose propensity score stratification for analyzing a non-randomized social experiment. One notable disadvantage of the propensity score methods is that they rely on an explicit propensity score model and are biased and highly variable if the model is misspecified (Kang and Schafer, 2007). The second type uses calibration weighting (Deville and Särndal, 1992, Kott, 2006, Chen, Valliant and Elliott, 2018, Chen et al., 2019). This technique calibrates auxiliary information in the non-probability sample with that in the probability sample, so that after calibration the weighted distribution of the non-probability sample is similar to that of the target population (DiSogra et al., 2011). The third type is mass imputation, which imputes the missing values for all units in the probability sample. In the usual imputation for missing data analysis, the respondents in the sample constitute a training dataset for developing an imputation model. In the mass imputation, an independent non-probability sample is used as a training dataset, and imputation is applied to all units in the probability sample; see, e.g., Breidt et al. (1996), Rivers (2007), Kim and Rao (2012), Chipperfield et al. (2012), Bethlehem (2016), and Yang and Kim (2018).

Let $X\in\mathbb{R}^{p}$ be a vector of auxiliary variables (including an intercept) that are available from two data sources, and let $Y$ be a general-type study variable of interest. We consider combining a probability sample with $X$, referred to as Sample A, and a non-probability sample with $(X,Y)$, referred to as Sample B, to estimate $\mu$ the population mean of $Y$. Because the sampling mechanism of a non-probability sample is unknown, the target population quantity is not identifiable in general. Researchers rely on an identification strategy that requires a non-informative sampling assumption imposed on the non-probability sample. To ensure this assumption holds, researchers should control for all covariates that are predictors of both sampling and the outcome variable. In practice, subject matter experts will recommend a rich set of potential useful variables but will not identify the exact variables to adjust for. In the presence of many auxiliary variables, variable selection becomes important, because existing methods may become unstable or even infeasible, and irrelevant auxiliary variables can introduce a large variability in estimation. There is a large literature on variable selection methods for prediction, but little work on variable selection for data integration that can successfully recognize the strengths and the limitations of each data source and utilize all information captured for finite population inference. Gao and Carroll (2017) propose a pseudo-likelihood approach to combining multiple non-survey data with high dimensionality; this approach requires all likelihoods be correctly specified and therefore is sensitive to model misspecification. Chen, Valliant and Elliott (2018) propose a model-assisted calibration approach using LASSO; this approach relies on a correctly specified outcome model. Up to our knowledge, robust inference has not been addressed in the context of data integration with high-dimensional data.

We propose a doubly robust variable selection and estimation strategy that harnesses the representativeness of the probability sample and the outcome information in the non-probability sample. The double robustness entails that the final estimator is consistent for the true value if either the probability of selection into the non-probability sample, referred to as the sampling score, or the outcome model is correctly specified, not necessarily both (a double robustness condition); see, e.g., Bang and Robins (2005), Tsiatis (2006), Cao et al. (2009), and Han and Wang (2013). To handle potentially high-dimensional covariates, our strategy separates the variable selection step and the estimation step for finite population mean to achieve two different goals.

In the first step, we select a set of variables that are important predictors of either the sampling score or the outcome model by penalized estimating equations. Following most of the empirical literature, we assume the sampling score follows a logistic regression model with the unknown parameter $\alpha\in\mathbb{R}^{p}$ and the outcome follows a generalized linear model (accommodating different types of the outcome) with the unknown parameter $\beta\in\mathbb{R}^{p}$. Importantly, we separate the estimating equations for $\alpha$ and $\beta$ in order to achieve stability in variable selection under the double robustness condition. Specifically, we construct the estimating equation for $\alpha$ by calibrating the weighted average of $X$ from Sample B, weighted by the inverse of the sampling score, to the design weighted average of $X$ from Sample A (i.e., a design estimate of population mean of $X$). We construct the estimating equation for $\beta$ by minimizing the standard least squared error loss under the outcome model. To establish the selection properties, we consider the “large $n$, diverging $p$” framework. To the best of our knowledge, the asymptotic properties of penalized estimating estimation based on survey data have not been studied in the literature. Our major technical challenge is that under the finite population framework, the sampling indicator of Sample A may not be independent even under simple random sampling. To overcome this challenge, we construct martingale random variables with a weak dependence that allows applying Bernstein inequality. This construction is innovative and crucial in establishing our new selection consistency result.

In the second step, we consider a doubly robust estimator of $\mu$, $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha},\widehat{\beta})$, and re-estimate $({\alpha},{\beta})$ based on the joint set of covariates selected from the first step. We propose using different estimating equations for $(\alpha,\beta)$, derived by minimizing the asymptotic squared bias of $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha},\widehat{\beta})$. This estimation strategy is not new; see, e.g., Kim and Haziza (2014) for missing data analyses in low-dimensional data; however, we demonstrate its new role in high-dimensional data to mitigate the possible selection error in the first step. In essence, our strategy for estimating $(\alpha,\beta)$ renders the first order term in the Taylor expansion of $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha},\widehat{\beta})$ with respect to $(\alpha,\beta)$ to be exactly zero, and the remaining terms are negligible under regularity conditions. This estimating strategy makes the doubly robust estimator root-$n$ consistent if either the sampling probability or the outcome model is correctly specified. This also enables us to construct a simple and consistent variance estimator allowing for doubly robust inferences. Importantly, the proposed estimator allows model misspecification of either the sampling score or the outcome model. In the existing high-dimensional causal inference literature, the doubly robust estimators have been shown to be robust to selection errors using penalization (Farrell, 2015) or approximation errors using machine learning (Chernozhukov et al., 2018). However, this double robustness feature requires both nuisance models to be correctly specified. We relax this requirement allowing one of the nuisance models to be misspecified. We clarify that even though the set of variables for estimation may include the variables that are solely related to the sampling score but not the outcome and therefore may harm efficiency of estimating $\mu$ (De Luna et al., 2011, Patrick et al., 2011), it is important to include these variables for $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha},\widehat{\beta})$ to achieve consistency in the case when the outcome model is misspecified and the sampling score model is correctly specified; see Section 6.

The paper proceeds as follows. Section 2 provides the basic setup of the paper. Section 3 presents the proposed two-step procedure for variable selection and doubly robust estimation of the finite population mean. Section 4 describes the computation algorithm for solving penalized estimating equations. Section 5 presents the theoretical properties for variable selection and doubly robust estimation. Section 6 reports simulation results that illustrate the finite-sample performance of the proposed method. In Section 7, we present an application to analyze a non-probability sample collected by the Pew Research Centre. We relegate all proofs to the supplementary material.

2 Basic Setup

2.1 Notation: Two Samples

Let $\mathcal{U}=\{1,\ldots,N\}$ be the index set of $N$ units for the finite population, with $N$ being the known population size. The finite population consists of $\mathcal{F}_{N}=\{(X_{i},Y_{i}):i\in\mathcal{U}\}$. Let the parameter of interest be the finite population mean $\mu=N^{-1}\sum_{i=1}^{N}Y_{i}$. We consider two data sources: one from a probability sample, referred to as Sample A, and the other one from a non-probability sample, referred to as Sample B. Table 1 illustrates the observed data structure. Sample A consists of observations $\mathcal{O}_{A}=\{(d_{A,i}=\pi_{A,i}^{-1},X_{i}):i\in\mathcal{A}\}$ with sample size $n_{A},$ where $\pi_{A,i}=P(i\in\mathcal{A})$ is known throughout Sample A, and Sample B consists of observations $\mathcal{O}_{B}=\{(X_{i},Y_{i}):i\in\mathcal{B}\}$ with sample size $n_{B}$. We define $I_{A,i}$ and $I_{B,i}$ to be the indicators of selection to Sample A and Sample B, respectively. Although the non-probability sample contains rich information on $(X,Y)$, the sampling mechanism is unknown, and therefore we cannot compute the first-order inclusion probability for Horvitz–Thompson estimation. The naive estimators without adjusting for the sampling process are subject to selection biases (Meng, 2018). On the other hand, although the probability sample with sampling weights represents the finite population, it does not observe the study variable of interest.

2.2 An Identification Assumption

Before presenting the proposed methodology for integrating the two data sources, we first discuss the identification assumption. Let $f(Y\mid X)$ be the conditional distribution of $Y$ given $X$ in the superpopulation model $\zeta$ that generates the finite population. We make the following primary assumption.

Assumption 1

(i) The sampling indicator $I_{B}$ of Sample B and the response variable $Y$ is independent given $X$; i.e. $P(I_{B}=1\mid X,Y)=P(I_{B}=1\mid X)$, referred to as the sampling score $\pi_{B}(X)$, and (ii) $\pi_{B}(X)>N^{\gamma-1}\delta_{B}>0$ for all $X$, where $\gamma\in(2/3,1]$.

Assumption 1 (i) implies that $E(Y\mid X)=E(Y\mid X,I_{B}=1)$, denoted by $m(X)$, can be estimated based solely on Sample B. Assumption 1 (ii) specifies a lower bound of $\pi_{B}(X)$ for the technicality in Section 5. A standard condition in the literature imposes a strict positivity in the sense that $\pi_{B}(X)>\delta_{B}>0$; however, it implies that $n_{B}^{-1}=O(N^{-1})$, which may be restrictive in survey sampling. Here, we relax this condition and allow $n_{B}^{-1}=O(N^{-\gamma})$, where $\gamma$ can be strictly less than $1$.

Assumption 1 is a key assumption for identification. Under Assumption 1, $E(\mu)$ is identifiable based on Sample A by $E\{I_{A}m(X)\}$ or Sample B by $E\{I_{B}Y/\pi_{B}(X)\}$. However, this assumption is not verifiable from the observed data. To ensure this assumption holds, researchers often consider many potentially predictors for the sampling indicator $I_{B}$ or the outcome $Y$, resulting in a rich set of variables in $X$.

2.3 Existing Estimators

In practice, the sampling score function $\pi_{B}(X)$ and the outcome mean function $m(X)$ are unknown and need to be estimated from the data. Let $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)$ and $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ be the posited models for $\pi_{B}(X)$ and $m(X)$, respectively, where $\alpha$ and $\beta$ are unknown parameters. Researchers have proposed various estimators for $\mu$ requiring different model assumptions and estimation strategies. We provide examples below and discuss their properties and limitations.

Example 1** (Inverse probability of sampling score weighting)**

Given an estimator $\widehat{\alpha}$, the inverse probability of sampling score weighting estimator is

[TABLE]

There are different approaches to obtain $\widehat{\alpha}$. Following Valliant and Dever (2011), one can obtain $\widehat{\alpha}$ by fitting the sampling score model based on the blended data $\mathcal{O}_{A}\cup\mathcal{O}_{B}=\{(d_{A,i},X_{i},I_{i}=0):i\in\mathcal{A}\}\cup\{(X_{i},I_{i}=1):i\in\mathcal{B}\}$, weighted by the design weights from Sample A. The resulting estimator $\widehat{\alpha}$ is valid if the size of Sample B is relatively small (Valliant and Dever, 2011). Elliott and Valliant (2017) propose an alternative strategy based on the Bayes rule: $\pi_{B}(X)\propto P(I_{A}=1\mid X)O_{B}(X)$, where $O_{B}(X)=P(I_{B}=1\mid X)/P(I_{B}=0\mid X)$ is the odds of selection into Sample B among the blended sample. This approach does not require the size of Sample B to be small; however, if $X$ does not correspond to the design variables for Sample A, it requires positing an additional model for $P(I_{A}=1\mid X)$. More importantly variable selection based on this approach is not straightforward in the setting with a high-dimensional $X$. To obtain $\hat{\alpha}$, we use the following estimating equation for $\alpha$:

[TABLE]

for some $h(X_{i};\alpha)$ such that (2) has a unique solution. Kott (2019) advocated using $h(X;\alpha)=X$ and Chen, Li and Wu (2018) advocated using $h(X;\alpha)=\pi(X;\alpha)\cdot X$. The justification for $\widehat{\mu}_{\mathrm{IPW}}$ relies on the correct specification of $\pi_{B}(X)$ and the consistency of $\widehat{\alpha}$. If $\pi_{B}(X_{i}^{\mathrm{\scriptscriptstyle T}}\alpha)$ is misspecified or $\widehat{\alpha}$ is inconsistent, $\widehat{\mu}_{\mathrm{IPW}}$ is biased.

Example 2** (Outcome regression based on Sample A)**

The outcome regression estimator is

[TABLE]

where $\widehat{\beta}$ is obtained by fitting the outcome model based solely on $\mathcal{O}_{B}=\{(X_{i},Y_{i}):i\in\mathcal{B}\}$ under Assumption 1.

The justification for $\widehat{\mu}_{\mathrm{reg}}$ relies on the correct specification of $m(X^{T}\beta)$ and the consistency of $\widehat{\beta}$. If $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ is misspecified or $\widehat{\beta}$ is inconsistent, $\widehat{\mu}_{\mathrm{reg}}$ can be biased.

Example 3** (Calibration weighting)**

The calibration weighting estimator is

[TABLE]

where $\{\omega_{i}:i\in\mathcal{S}_{B}\}$ satisfies constraint (i) $\sum_{i\in\mathcal{S}_{B}}\omega_{i}X_{i}=\sum_{i\in\mathcal{S}_{A}}d_{A,i}X_{i},$ or constraint (ii) $\sum_{i\in\mathcal{S}_{B}}\omega_{i}m(X_{i};\widehat{\beta})=\sum_{i\in\mathcal{S}_{A}}d_{A,i}m(X_{i};\widehat{\beta})$ in a model-assisted approach (McConville et al., 2017, Chen, Valliant and Elliott, 2018, Chen et al., 2019).

The justification for $\widehat{\mu}_{\mathrm{cal}}$ subject to constraint (i) relies on the linearity of the outcome model, i.e., $m(X)=X^{\mathrm{\scriptscriptstyle T}}\beta^{*}$ for some $\beta^{*}$, or the linearity of the inverse probability of sampling weight, i.e., $\pi_{B}(X)^{-1}=X^{\mathrm{\scriptscriptstyle T}}\alpha^{*}$ for some $\alpha^{*}$ (Fuller, 2009; Theorem 5.1). The linearity conditions are unlikely to hold for non-continuous variables. In these cases, $\widehat{\mu}_{\mathrm{cal}}$ is biased. The justification for $\widehat{\mu}_{\mathrm{cal}}$ subject to constraint (ii) relies on $m(X;\beta)$ being correctly specified in the data integration problem.

Example 4** (Doubly robust estimator)**

The doubly robust estimator is

[TABLE]

The estimator $\widehat{\mu}_{\mathrm{dr}}$ is doubly robust with fixed-dimensional $X$ (Chen, Li and Wu, 2018), in the sense that it achieves the consistency if either $\pi_{B}(X_{i}^{\mathrm{\scriptscriptstyle T}}\alpha)$ or $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ is correctly specified, but not necessarily both. The double robustness is attractive; therefore, we shall investigate the potential of $\widehat{\mu}_{\mathrm{dr}}$ in high-dimensional setup.

3 Methodology in High-dimensional Data

A major challenge arises in the presence of a large number of covariates, not all of them are necessary for making inference of the population mean of the outcome. This necessitates variable selection. For simplicity of exposition, we introduce the following notation. For any vector $\alpha\in\mathbb{R}^{p}$, denote the number of nonzero elements in $\alpha$ as $||\alpha||_{0}=\sum_{j=1}^{p}I(\alpha_{j}\neq 0)$, the $L_{1}$-norm as $||\alpha||_{1}=\sum_{j=1}^{p}|\alpha_{j}|$, the $L_{2}$-norm as $||\alpha||_{2}=\sqrt{\sum_{j=1}^{p}\alpha_{j}^{2}}$, and the $L_{\infty}$-norm as $||\alpha||_{\infty}=\max_{j=1}^{p}|\alpha_{j}|$. For any $\mathcal{J}\subseteq\{1,\ldots,p\}$, let $\alpha_{\mathcal{J}}$ be the sub-vector of $\alpha$ formed by elements of $\alpha$ whose indexes are in $\mathcal{J}$. Let $\mathcal{J}^{c}$ be the complement of $\mathcal{J}$. For any $\mathcal{J}_{1},\mathcal{J}_{2}\subseteq\{1,\ldots,p\}$ and matrix $\Sigma\in\mathbb{R}^{p\times p}$, let $\Sigma_{\mathcal{J}_{1},\mathcal{J}_{2}}$ be the sub-matrix of $\Sigma$ formed by rows in $\mathcal{J}_{1}$ and columns in $\mathcal{J}_{2}$. Following the literature on variable selection, we can first standardize the covariates so that approximately they have variances equal to one to stabilize the variable selection procedure. We make the following modeling assumptions.

Assumption 2** (Sampling score model)**

The sampling mechanism of Sample B, $\pi_{B}(X)$, follows a logistic regression model $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)$; i.e., $\text{logit}\left\{\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)\right\}=X^{\mathrm{\scriptscriptstyle T}}\alpha$ for $\alpha\in\mathbb{R}^{p}$.

Assumption 3** (Outcome model)**

The outcome mean function $m(X)$ follows a generalized linear regression model; i.e., $m(X)=m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ for $\beta\in\mathbb{R}^{p}$, where $m(\cdot)$ is a link function by an abuse the notation.

Define $\alpha^{*}$ to be the $p$-dimensional parameter that minimizes the Kullback-Leibler divergence

[TABLE]

and $\beta^{*}=\arg\min_{\beta}E\left[\{Y-m(X^{\mathrm{\scriptscriptstyle T}}\beta)\}^{2}\right].$

In Assumption 2, we adopt the logistic regression model for the sampling score following most of the empirical literature; but our framework can be extended to the case with other parameter models such as the probit model. The models $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)$ and $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ are working models and they may be misspecified. If the sampling score model is correctly specified, $\pi_{B}(X)=\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha^{*})$. If the outcome model is correctly specified, $m(X)=m(X^{\mathrm{\scriptscriptstyle T}}\beta^{*})$.

The proposed procedure consists of two steps: the first step selects important variables in the sampling score model and the outcome model, and the second step focuses on doubly robust estimation of the population mean.

In the first step, we propose solving penalized estimating equations for variable selection. Using (2) with $h(X;\alpha)=X$, we define the estimating function for $\alpha$ as

[TABLE]

To select important variables in $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$, under Assumption 1, we have $E(Y\mid X)=E(Y\mid X,I_{B}=1)$. Therefore, we define the estimating function for $\beta$ as

[TABLE]

Let $U(\theta)=(U_{1}(\alpha)^{\mathrm{\scriptscriptstyle T}},U_{2}(\beta)^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}$ be the joint estimating function for $\theta=(\alpha^{\mathrm{\scriptscriptstyle T}},\beta^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}$. When $p$ is large, following Johnson et al. (2008), we consider solving the penalized estimating function

[TABLE]

for $(\alpha,\beta)$, where $q_{\lambda_{\alpha}}(\alpha)=\{q_{\lambda_{\alpha}}(|\alpha_{0}|),\ldots,q_{\lambda_{\alpha}}(|\alpha_{p}|)\}^{\mathrm{\scriptscriptstyle T}}$ and $q_{\lambda_{\beta}}(\beta)=\{q_{\lambda_{\beta}}(|\beta_{0}|),\ldots,q_{\lambda_{\beta}}(|\beta_{p}|)\}^{\mathrm{\scriptscriptstyle T}}$ are some continuous functions, $q_{\lambda_{\alpha}}(|\alpha|)\mathrm{sign}(\alpha)$ is the element-wise product of $q_{\lambda_{\alpha}}(\alpha)$ and $\mathrm{sign}(\alpha),$ and $q_{\lambda_{\beta}}(|\beta|)\mathrm{sign}(\beta)$ is the element-wise product of $q_{\lambda_{\beta}}(\beta)$ and $\mathrm{sign}(\beta)$. We let $q_{\lambda}(x)=\mathrm{d}p_{\lambda}(x)/\mathrm{d}x$, where $p_{\lambda}(x)$ is some penalization function. Although the same discussion applies to different non-concave penalty functions, we specify $p_{\lambda}(x)$ to be a folded-concave smoothly clipped absolute deviation (SCAD) penalty function (Fan and Lv, 2011). Accordingly, we have

[TABLE]

for $a>0$, where $(\cdot)_{+}$ is the truncated linear function; i.e., if $x\geq 0$, $(x)_{+}=x$, and if $x<0$, $(x)_{+}=0$. We use $a=3.7$. Fan and Li (2001) demonstrate that with $a=3.7$, the SCAD selector has a good performance based on simulation. Thereafter, this choice has become standard in the literature and become a default choice in many softwares such as “ncvreg” in R. We select the variables if the corresponding estimates of (6) are nonzero in either the sampling score or the outcome model, indexed by $\mathcal{C}$.

Remark 1

To help understand the penalized estimating equation, we discuss two scenarios. If $|\alpha_{j}|$ is large, then $q_{\lambda_{\alpha}}(|\alpha_{j}|)$ is zero, and therefore $U_{1,j}(\alpha)$ is not penalized. Whereas, if $|\alpha_{j}|$ is small but nonzero, then $q_{\lambda_{\alpha}}(|\alpha_{j}|)$ is large, and therefore $U_{1,j}(\alpha)$ is penalized with a penalty term. The penalty term then forces $\widehat{\alpha}_{j}$ to be zero and excludes the $j$th element in $X$ from the final selected set of variables. The same discussion applies to $U_{2}(\beta)$ and $q_{\lambda_{\beta}}(|\beta|)$.

In the second step, we consider the estimator of the population mean $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha},\widehat{\beta})$ in (5) with $(\widehat{\alpha},\widehat{\beta})$ re-estimated based on $X_{\mathcal{C}}$. As we will show in Section 5, $\mathcal{C}$ contains the true important variables in either the sampling score model or the outcome model with probability approaching one (the oracle property). Therefore, if either the sampling score model or the outcome model is correctly specified, the asymptotic bias of $\widehat{\mu}_{\mathrm{dr}}(\alpha^{*},\beta^{*})$ is zero; however, if both models are misspecified, the asymptotic bias of $\widehat{\mu}_{\mathrm{dr}}(\alpha^{*},\beta^{*})$ is

[TABLE]

In order to minimize $\{\text{a.bias}(\alpha,\beta)\}^{2}$, we consider the estimating function

[TABLE]

and the corresponding empirical estimating function

[TABLE]

for estimating $(\alpha,\beta)$, constrained on $\{(\alpha^{\mathrm{\scriptscriptstyle T}},\beta^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}\in\mathbb{R}^{2p}:\alpha_{\mathcal{C}^{c}}=0,\beta_{\mathcal{C}^{c}}=0\}$. Equation (9) is doubly robust in the sense that $J(\alpha^{*},\beta^{*})$ is unbiased if either $\pi(X^{\mathrm{\scriptscriptstyle T}}\alpha)$ or $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ is correctly specified, not necessarily both (Kim and Haziza, 2014).

Remark 2

*The two steps use different estimating functions (6) and (9), respectively, for selection and estimation with the following advantages. First, (6) separates the selection for $\alpha$ and $\beta$ in $U_{1}(\alpha)$ and $U_{2}(\beta)$, so it stabilizes the selection procedure if either the sampling score model or the outcome model is misspecified. Second, using (9) for estimation leads to an attractive feature for inference about $\mu$. We clarify that although the joint estimating function (9) is motivated by minimizing the asymptotic bias $\text{a.bias}(\alpha^{*},\beta^{*})$ if both nuisance models are misspecified, we do not expect that the proposed estimator for $\mu$ is unbiased in this case. Instead, we show the advantage of (9) in the case if either the sampling probability or the outcome model is correctly specified in high-dimensional data. It is well-known that post-selection inference is notoriously difficult even when both models are correctly specified because the estimation step is based on a random set of variables being selected. We show that our estimating strategy based on (9) mitigates the possible first-step selection error and renders $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha},\widehat{\beta})$ root-$n$ consistent if either the sampling probability or the outcome model is correctly specified in high-dimensional data. Heuristically this is achieved because the first Taylor expansion term is set to be zero due to (8). We relegate the details to Section 5.

In summary, our two-step procedure for variable selection and estimation is as follows. spacing

• Step$\ 1.$

To facilitate joint selection of variables for the sampling score and outcome, solve the penalized joint estimating equations $U^{\mathrm{p}}(\alpha,\beta)=0$ in (6), denoted by $(\widetilde{\alpha},\widetilde{\beta})$. Let $\widehat{\mathcal{M}}_{\alpha}=\{j:\widetilde{\alpha}_{j}\neq 0\}$ and $\widehat{\mathcal{M}}_{\beta}=\{j:\widetilde{\beta}_{j}\neq 0\}$.

• Step$\ 2.$

Let the set of variables for estimation be $\mathcal{C}=\widehat{\mathcal{M}}_{\alpha}\cup\widehat{\mathcal{M}}_{\beta}$. Obtain the proposed estimator as

[TABLE]

where $\widehat{\alpha}$ and $\widehat{\beta}$ are obtained by solving the joint estimating equations (9) for $\alpha$ and $\beta$ with $\alpha_{\mathcal{C}^{c}}=0$ and $\beta_{\mathcal{C}^{c}}=0$.

Remark 3

Variable selection circumvents the instability or infeasibility of direct estimation of $(\alpha,\beta)$ with high-dimensional $X$. Moreover, in Step 2 for estimation, we consider a union of covariates $X_{\mathcal{C}}$, where $\mathcal{C}=\widehat{\mathcal{M}}_{\alpha}\cup\widehat{\mathcal{M}}_{\beta}$. It is worth comparing this choice with two other common choices in the literature. First, one considers separate sets of variables for the two models; i.e., the sampling score is fitted based on $\widehat{\mathcal{M}}_{\alpha}$, and the outcome model is fitted based on $\widehat{\mathcal{M}}_{\beta}$. However, we note that in the joint estimating equation (9), $J_{1}(\alpha,\beta)$ and $J_{2}(\alpha,\beta)$ should have the same dimension, otherwise, it is possible that (9) does not guarantee that there exists a solution. This is obvious if one considers a linear outcome model. Moreover, Brookhart et al. (2006) and Shortreed and Ertefaie (2017) show that including variables that are related to the outcome in the propensity score model will increase the precision of the estimated average treatment effect without increasing bias. This implies that an efficient variable selection and estimation method should take into account both sampling-covariate and outcome-covariate relationships. As a result, $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha},\widehat{\beta})$ may have a better performance than the oracle estimator which uses the true important variables in the sampling score and the outcome model. This is particularly true when one of the models is misspecified. Our simulation study in Section 6 demonstrates that $\widehat{\mu}_{\mathrm{dr}}$ with variable selection has a similar performance as the orcale estimator for the continuous outcome and outperforms the oracle estimator for the binary outcome. Second, many authors have suggested that including predictors that are solely related to the sampling score but not the outcome may harm efficiency (De Luna et al., 2011, Patrick et al., 2011). However, this strategy is effective when both the sampling score and outcome models are correctly specified. When the sampling score model is correctly specified but the outcome model is misspecified, restricting the variables to be the outcome predictors may render the sampling score “misspecified” by using the wrong set of variables. The simulation study suggests that $\widehat{\mu}_{\mathrm{p}-dr}$ restricted to the set of variables in $\widehat{\mathcal{M}}_{\beta}$ is not doubly robust.

4 Computation

In this section, we discuss the computation for solving the penalized estimating function (6). Following Johnson et al. (2008), we use an iterative algorithm that combines the Newton–Raphson algorithm for solving estimating equation and the minorization-maximization algorithm for non-convex penalty of Hunter and Li (2005).

First, by the minorization-maximization algorithm, the penalized estimator $\widetilde{\theta}=(\widetilde{\alpha},\widetilde{\beta})$ solving (6) satisfies

[TABLE]

for $\epsilon$ is a predefined small number. In our implementation, we choose $\epsilon$ to be $10^{-6}$.

Second, we solve (11) by the Newton-Raphson algorithm. It may be challenging to implement the Newton-Raphson algorithm directly, because it involves inverting a large matrix. For ease and stability in those cases, we can use a coordinate decent algorithm (Friedman et al., 2007) by cycling through and updating each of the coordinates.

Following most of the empirical literature, we assume that $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)$ follows a logistic regression model. Define $m^{(k)}(t)=\mathrm{d}^{k}m(t)/\mathrm{d}^{k}t$ for $k\geq 1$. We denote

[TABLE]

and

[TABLE]

Let $\theta$ start at an initial value $\widetilde{\theta}^{[0]}$. With the other coordinates fixed, the $k$th Newton-Raphson update for $\theta_{j}$ ($j=1,\ldots,2p$), the $j$th element of $\theta$, is

[TABLE]

where $\nabla_{jj}(\theta)$ and $E_{jj}(\theta)$ are the $j$th diagonal elements in $\nabla(\theta)$ and $E(\theta)$, respectively. The procedure cycles through all the $2p$ elements of $\theta$ and is repeated until convergence.

We use $K$-fold cross-validation to select the tuning parameter $(\lambda_{\alpha},\lambda_{\beta})$. To be specific, we partition both samples into approximately $K$ equal sized subsets and pair subsets of Sample A and subsets of Sample B randomly. Of the $K$ pairs, we retain one single pair as the validation data and the remaining $K-1$ pairs as the training data. We fit the models based on the training data and estimate the loss function based on the validation data. We repeat the process $K$ times, with each of the $K$ pairs used exactly once as the validation data. Finally, we aggregate the $K$ estimated loss function. We select the tuning parameter as the one that minimizes the aggregated loss function over a pre-specified grid.

Because the weighting estimator uses the sampling score $\pi_{B}(X)$ to calibrate the distribution of $X_{\mathcal{C}}$ between Sample B and the target population, we use the following loss function for selecting $\lambda_{\alpha}$:

[TABLE]

where $\widetilde{\alpha}(\lambda_{\alpha})$ is the penalized estimator $\widetilde{\alpha}$ with the tuning parameter $\lambda_{\alpha}$. We use the prediction error loss function for selecting $\lambda_{\beta}$:

[TABLE]

where $\widetilde{\beta}(\lambda_{\beta})$ is the penalized estimator $\widetilde{\beta}$ with the tuning parameter $\lambda_{\beta}$.

5 Asymptotic Results for Variable Selection and Estimation

We establish the asymptotic properties for the proposed double variable selection and doubly robust estimation. We can establish theoretical results for general sampling mechanisms for Sample A requiring specific regularity conditions. In this section, for technical convenience, we assume that Sample A is collected by simple random sampling or Poisson sampling with the following regularity conditions.

Assumption 4

For all $1\leq i\leq N$, $\pi_{A,i}\geq N^{\gamma-1}\delta_{A}>0$, where $\gamma\in(2/3,1]$.

Similar to Assumption 1 (ii), we relax the strict positivity on $\pi_{A,i}$ and render $n_{A}=O(N^{\gamma})$ for $\gamma$ possibly strictly less than $1$. Let $n=\min(n_{A},n_{B})$, which is $O(N^{\gamma})$ under Assumptions 1 and 4.

Let the support of model parameters be

[TABLE]

Define $s_{\alpha}=||\alpha^{*}||_{0}$, $s_{\beta}=||\beta^{*}||_{0}$ , $s_{\theta}=s_{\alpha}+s_{\beta}$, and $\lambda_{\theta}=\min(\lambda_{\alpha},\lambda_{\beta})$.

Assumption 5

The following regularity conditions hold.

(A1)

The parameter $\theta$ belongs to a compact subset in $\mathbb{R}^{2p}$, and $\theta^{*}$ lies in the interior of the compact subset.

(A2)

$\{X_{i}:i\in\mathcal{U}\}$* are fixed and uniformly bounded.*

(A3)

There exist constants $c_{1}$ and $c_{2}$ such that

[TABLE]

where $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ are the minimum and the maximum eigenvalue of a matrix, respectively.

(A4)

*Let $\epsilon_{i}(\beta)=Y_{i}-m(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)$ be the *$i$th residual. There exists a constant $c_{3}$ such that $E\{|\epsilon_{i}(\beta^{*})|^{2+\delta}\}\leq c_{3}$ for all $1\leq i\leq N$ and some $\delta>0$. There exist constants $c_{4}$ and $c_{5}$ such that $E[\exp\{c_{4}|\epsilon_{i}(\beta^{*})|\}\mid X_{i}]\leq c_{5}$ for all $1\leq i\leq N$.

(A5)

$m^{(1)}(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)$, $m^{(2)}(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)$, and $m^{(3)}(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)$ are uniformly bounded away from $\infty$ on $\mathcal{N}_{\theta,\tau}=\{\theta\in\mathbb{R}^{2p}:||\theta_{\mathcal{M}_{\theta}}-\theta_{\mathcal{M}_{\theta}}^{*}||\leq\tau\sqrt{s_{\theta}/n},\theta_{\mathcal{M}_{\theta}^{c}}=0\}$ for some $\tau>0$.

(A6)

$\min_{j\in\mathcal{M}_{\alpha}}|\alpha_{j}^{*}|/\lambda_{\alpha}\rightarrow\infty$* and $\min_{k\in\mathcal{M}_{\beta}}|\beta_{k}^{*}|/\lambda_{\beta}\rightarrow\infty,$ as $n\rightarrow\infty$.*

(A7)

$s_{\theta}=o(n^{1/3})$, $\lambda_{\alpha},\lambda_{\beta}\rightarrow 0$, $(\log n)^{2}=o(n\lambda_{\theta}^{2})$, $\log(p)=o\left\{n\lambda_{\theta}^{2}/\left(\log n\right){}^{2}\right\}$, $ps_{\theta}^{4}(\log n)^{6}=o(n^{3}\lambda_{\theta}^{2})$, $ps_{\theta}^{4}(\log n)^{8}=o(n^{4}\lambda_{\theta}^{4})$, as $n\rightarrow\infty$.

These assumptions are typical in the penalization literature. (A2) specifies a fixed design which is well suited under the finite population inference framework. (A4) holds for Gaussian distribution, sub-Gaussian distribution, and so on. (A5) holds for common models. (A7) specifies the restrictions on the dimension of covariates $p$ and the dimension of the true nonzero coefficients $s_{\theta}$. To gain insight, when the true model size $s_{\theta}$ is fixed, (A7) holds for $p=O(n)$, i.e., $p$ can be the same size as $n$.

We establish the asymptotic properties of the penalized estimating equation procedure.

Theorem 1

Under Assumptions 1–5, there exists an approximate penalized solution $\widetilde{\theta}$, which satisfies the selection consistency properties:

[TABLE]

and

[TABLE]

as $n\rightarrow\infty.$

Results (15) and (16) imply that $U(\widetilde{\theta})=o_{P}\left(\lambda_{\theta}/\log n\right)$. Results (17) and (18) imply that with probability approaching to one, the penalized estimating equation procedure would not over-select irrelevant variables and estimate the true nonzero coefficients at the $\sqrt{s_{\theta}/n}$ convergence rate, which is the so-called oracle property of variable selection.

Remark 4

It is worth discussing the relationship of Theorem 1 to existing variable selection methods in the survey literature. Based on a single probability sample source, McConville et al. (2017) propose a model-assisted survey regression estimator of finite-population totals using the LASSO to improve the efficiency. Chen, Valliant and Elliott (2018) and Chen et al. (2019) propose model-assisted calibration estimators using the LASSO based on nonprobability samples integrating with auxiliary known totals or probability samples, respectively. However, their methods require the working outcome model to include sufficient population information and therefore are not doubly robust. To the best of our knowledge, our paper is the first to propose doubly robust inference of finite population means after variable selection.

We now establish the asymptotic properties of $\widehat{\mu}_{\text{p-dr}}(\widehat{\alpha},\widehat{\beta})$. Define a sequence of events $\mathcal{D}_{n}=\{\mathcal{M}_{\theta}\subset\mathcal{C}\}$, where we emphasize that $\mathcal{D}_{n}$ depends on $n$ although we suppress the dependence of $\mathcal{M}_{\theta}$ and $\mathcal{C}$ on $n$. Following the same argument for (18), given the event $\mathcal{D}_{n}$, we have $\{(\widehat{\alpha}-\alpha^{*})^{\mathrm{\scriptscriptstyle T}},(\widehat{\beta}-\beta^{*})^{\mathrm{\scriptscriptstyle T}}\}=O_{p}(\sqrt{s_{\theta}/n}).$ Combining with $P(\mathcal{D}_{n})\rightarrow 1$, we have

[TABLE]

By Taylor expansion,

[TABLE]

where $\widehat{\mu}_{\text{p-dr}}(\alpha,\beta)$ is defined in (10). Equation (28) follows because we solve (9) for $(\alpha,\beta)$. Equation (29) follows because of (19) and Assumption 5 (A7). As a result, the way for estimating $(\alpha^{*},\beta^{*})$ leads to the asymptotic equivalence between $\widehat{\mu}_{\text{p-dr}}(\widehat{\alpha},\widehat{\beta})$ and $\widehat{\mu}_{\text{p-dr}}(\alpha^{*},\beta^{*})$.

Moreover, we show that $\widehat{\mu}_{\text{p-dr}}(\alpha^{*},\beta^{*})$ is asymptotically unbiased for $\mu$ under the double robustness condition. We note that

[TABLE]

If $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)$ is correctly specified, then $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha^{*})=\pi_{B}(X)$ and therefore (5) is zero; if $m(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)$ is correctly specified, then $m(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta^{*})=m(X_{i})$ and therefore (5) is zero.

Following the variance decomposition strategy of Shao and Steel (1999), the asymptotic variance of the linearized term is

[TABLE]

where the conditional distribution in $E(\cdot\mid I_{B},X,Y)$ and $V(\cdot\mid I_{B},X,Y)$ is the sampling distribution for Sample A. The first term $V_{1}$ is the sampling variance of the Horvitz–Thompson estimator. Thus,

[TABLE]

For the second term $V_{2}$, note that

[TABLE]

Thus,

[TABLE]

Theorem 2 below summarizes the asymptotic properties of $\widehat{\mu}_{\text{p-dr}}$.

Theorem 2

Under Assumptions 1–5, if either $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)$ or $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ is correctly specified,

[TABLE]

as $n\rightarrow\infty$, where $V=\lim_{n\rightarrow\infty}(V_{1}+V_{2}),$ $V_{1}$ and $V_{2}$ are defined in (30) and (31), respectively.

To estimate $V_{1}$, we can use the design-based variance estimator applied to $m(X_{i}^{\mathrm{\scriptscriptstyle T}}\widehat{\beta})$ as

[TABLE]

To estimate $V_{2},$ we further express $V_{2}$ as

[TABLE]

Let $\sigma^{2}(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta^{*})=E\left[\left\{Y_{i}-m(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta^{*})\right\}^{2}\right]$, and let $\widehat{\sigma}^{2}(X_{i})$ be a consistent estimator of $\sigma^{2}(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta^{*})$. We can then estimate $V_{2}$ by

[TABLE]

By the law of large numbers, $\widehat{V}_{2}$ is consistent for $V_{2}$ regardless whether one of $\pi_{B,i}(X_{i}^{\mathrm{\scriptscriptstyle T}}\alpha)$ or $\pi_{B,i}(X_{i}^{\mathrm{\scriptscriptstyle T}}\beta)$ is misspecificed, and therefore it is doubly robust.

Theorem 3** (Double robustness of $\widehat{V}$)**

Under Assumptions 1–5, if either $\pi_{B}(X^{\mathrm{\scriptscriptstyle T}}\alpha)$ or $m(X^{\mathrm{\scriptscriptstyle T}}\beta)$ is correctly specified, $\widehat{V}=\widehat{V}_{1}+\widehat{V}_{2}$ is consistent for $V$.

6 Simulation Study

6.1 Setup

In this section, we evaluate the finite-sample performance of the proposed procedure. We first generate a finite population $\mathcal{F}_{N}=\{(X_{i},Y_{i}):i=1,\ldots N\}$ with $N=10,000$, where $Y_{i}$ is a continuous or binary outcome variable, and $X_{i}=(1,X_{1,i},\ldots,X_{p-1,i})^{\mathrm{\scriptscriptstyle T}}$ is a $p$-dimensional vector of covariates with the first component being $1$ and other components independently generated from standard normal with mean [math] and variance $1$. We set $p=50$. From the finite population, we select a non-probability sample $\mathcal{B}$ of size $n_{B}\approx 2,000$, according to the inclusion indicator $I_{B,i}\sim$Ber$(\pi_{B,i})$. We select a probability sample $\mathcal{A}$ of the average size $n_{A}=500$ under Poisson sampling with $\pi_{A,i}\propto(0.25+|X_{1i}|+0.03|Y_{i}|)$. The parameter of interest is the population mean $\mu=N^{-1}\sum_{i=1}^{N}Y_{i}$.

For the non-probability sampling probability, we consider both linear and nonlinear sampling score models

PSM  I:

$\text{logit}(\pi_{B,i})=\alpha_{0}^{\mathrm{\scriptscriptstyle T}}X_{i}$, where $\alpha_{0}=(-2,1,1,1,1,0,0,0,\ldots,0)^{\mathrm{\scriptscriptstyle T}}$,

PSM  II:

$\text{logit}(\pi_{B,i})=3.5+\alpha_{0}^{\mathrm{\scriptscriptstyle T}}\log(X_{i}^{2})-\sin(X_{3,i}+X_{4,i})-X_{5,i}-X_{6,i}$, where $\alpha_{0}=(0,0,0,3,3,3,3,0,\ldots,0)^{\mathrm{\scriptscriptstyle T}}$.

For generating a continuous outcome variable $Y_{i}$, we consider both linear and nonlinear outcome models with $\beta_{0}=(1,0,0,1,1,1,1,0\ldots,0)^{\mathrm{\scriptscriptstyle T}}$:

OM  I:

$Y_{i}=\beta_{0}^{\mathrm{\scriptscriptstyle T}}X_{i}+\epsilon_{i},$ $\epsilon_{i}\sim\mathcal{N}(0,1)$,

OM  II:

$Y_{i}=1+\exp\left\{3\sin\left(\beta_{0}^{\mathrm{\scriptscriptstyle T}}X_{i}\right)\right\}+X_{5,i}+X_{6,i}+\epsilon_{i},$ $\epsilon_{i}\sim\mathcal{N}(0,1)$.

For generating a binary outcome variable $Y_{i}$, we consider both linear and nonlinear outcome models with $\beta_{0}=(1,0,0,3,3,3,3,0\ldots,0)^{\mathrm{\scriptscriptstyle T}}$,

OM  I:

$Y\sim$Ber$\{\pi_{Y}(X)\}$ with logit$\{\pi_{Y}(X)\}=\beta_{0}^{\mathrm{\scriptscriptstyle T}}X$,

OM  II:

$Y\sim$Ber$\{\pi_{Y}(X)\}$ with logit$\{\pi_{Y}(X)\}=2-\log\left\{\left(\beta_{0}^{\mathrm{\scriptscriptstyle T}}X\right)^{2}\right\}+2X_{5,i}+2X_{6,i}$.

We consider the following estimators:

Naive,$\ \hat{\mu}_{\mathrm{naive}}$,

the naive estimator using the simple average of $Y_{i}$ from Sample B, which provides the degree of the selection bias;

Oracle,

$\hat{\mu}_{\mathrm{ora}}$, the doubly robust estimator $\widehat{\mu}_{\mathrm{dr}}(\widehat{\alpha}_{\mathrm{ora}},\widehat{\beta}_{\mathrm{ora}}),$ where $\widehat{\alpha}_{\mathrm{ora}}$ and $\widehat{\beta}_{\mathrm{ora}}$ are based on the joint estimation restricting to the known important covariates for comparison purpose;

p-ipw,$\ \hat{\mu}_{\text{p-ipw}}$,

the penalized inverse probability of sampling weighting estimator $\hat{\mu}_{\mathrm{IPW}}=$ $N^{-1}\sum_{i\in\mathcal{B}}\widehat{\pi}_{B,i}^{-1}Y_{i},$ where $\widehat{\pi}_{B,i}=P(I_{B,i}=1\mid X_{i}^{\mathrm{\scriptscriptstyle T}}\widehat{\alpha})$, and $\widehat{\alpha}$ is obtained by a weighted penalized regression of $I_{B,i}$ on $X_{i}$ using the combined sample of A and B, weighted by the design weights;

p-reg,$\ \hat{\mu}_{\text{p-reg}}$,

the penalized regression estimator $\hat{\mu}_{\text{p-reg}}=N^{-1}\sum_{i\in\mathcal{A}}d_{A,i}m(X;\widehat{\beta})$, where $\widehat{\beta}$ is obtained by a penalized regression of $Y_{i}$ on $X_{i}$ based on Sample B;

p-dr0,$\ \hat{\mu}_{\text{p-dr}}$,

the penalized double estimating equation estimator based on the set of outcome predictors $\widehat{\mathcal{M}}_{\beta}$;

p-dr,$\ \hat{\mu}_{\text{p-dr}}$,

the proposed penalized double estimating equation estimator based on the union of sampling and outcome predictors $\widehat{\mathcal{M}}_{\alpha}\cup\widehat{\mathcal{M}}_{\beta}$;.

We also note that $\widehat{\mu}_{\mathrm{dr}}$ without variable selection is severely biased and unstable and therefore is excluded for comparison.

6.2 Simulation Results

All simulation results are based on $500$ Monte Carlo runs. Table 2 reports the selection performance of the proposed penalized estimating equation approach in terms of the proportion of the proposed procedure under-selecting (Under), over-selecting (Over), the average false negatives (FN: the average number of selected covariates that have the true zero coefficients), and the average false positives (FP: the average number of selected covariates that have the true zero coefficients). The proposed procedure selects all covariates with nonzero coefficients in both outcome model and the sampling score model under the true model specification. Moreover, the number of false positives is small under the true model specification.

Figure 1 displays the estimation simulation results for the continuous outcome. The naive estimator $\hat{\mu}_{\mathrm{naive}}$ shows large biases across scenarios. The oracle estimator $\hat{\mu}_{\mathrm{ora}}$ is doubly robust, in the sense that if either the outcome or the sampling score is correctly specified, it is unbiased. The penalized inverse probability of sampling weighting estimator $\hat{\mu}_{\text{p-ipw}}$ shows larges biases except for Scenario (ii). The weighted estimator $\widehat{\alpha}$ is based on the blended sample combining Sample A and Sample B, where the units in Sample A are weighted by the known sampling weights and the units in Sample B are weighted by $1$. This approach is justifiable only if the sampling rate of Sample B is relatively small compared the the population size. The penalized regression estimator $\hat{\mu}_{\text{p-reg}}$ is only singly robust. When the outcome model is misspecified as in Scenarios (ii) and (iv), it shows large biases. The proposed penalized double estimating equation estimator $\hat{\mu}_{\text{p-dr}}$ based on $\widehat{\mathcal{M}}_{\alpha}\cup\widehat{\mathcal{M}}_{\beta}$ is doubly robust, and its performance is comparable to the oracle estimator that requires knowing the true important variables. Moreover, $\hat{\mu}_{\text{p-dr}}$ is slightly more efficient than $\hat{\mu}_{\mathrm{ora}}$. This efficiency gain is due to using the union of covariates selected for the sampling score model and the outcome model. This phenomenon is consistent with the findings in Brookhart et al. (2006) and Shortreed and Ertefaie (2017). The proposed penalized double estimating equation estimator $\widehat{\mu}_{\mathrm{p-dr0}}$ based on $\widehat{\mathcal{M}}_{\beta}$ is slightly more efficient than $\widehat{\mu}_{\mathrm{p-dr}}$ based on $\widehat{\mathcal{M}}_{\alpha}\cup\widehat{\mathcal{M}}_{\beta}$ in Scenario (i) when both the outcome and sampling score models are correctly specified; however, $\widehat{\mu}_{\mathrm{p-dr0}}$ has a large bias in Scenario (ii) when the outcome model is misspecified and therefore is not doubly robust anymore; see Remark 3.

Figure 2 displays the estimation results for the binary outcome. The same discussion above applies here. Moreover, when the outcome model is incorrectly specified, the oracle estimator has a large variability. In this case, the proposed estimator outperforms the oracle estimator, because the variable selection step helps to stabilize the estimation performance.

Table 3 reports the simulation results for the coverage properties for the continuous outcome and binary outcome. Under the double robustness condition (i.e., if either the outcome model or the sampling score model is correctly specified), the coverage rates are close to the nominal coverage; while if both models are misspecified, the coverage rates are off the nominal coverage.

7 An Application

We analyze two datasets from the 2005 Pew Research Centre (PRC, http://www.pewresearch.org/) and the 2005 Behavioral Risk Factor Surveillance System (BRFSS). The goal of the PRC study was to evaluate the relationship between individuals and community (Chen, Li and Wu, 2018, Kim et al., 2018). The 2005 PRC dataset is from a non-probability sample provided by eight different vendors, which consists of $n_{B}=9,301$ subjects. We focus on two study variables, a continuous $Y_{1}$ (days had at least one drink last month) and a binary $Y_{2}$ (an indicator of voted local elections). The 2005 BRFSS sample is a probability sample, which consists of $n_{A}=441,456$ subjects with survey weights. This dataset does not have measurements on the study variables of interest; however, it contains a rich set of common covariates with the PRC dataset listed in Figure 3. To illustrate the heterogeneity in the study populations, Figure 3 contrasts the covariate means from the PRC data and the design-weighted covariate means (i.e., the estimated population covariate means) from the BRFSS dataset. The covariate distributions from the PRC sample and the BRFSS sample are considerably different, e.g., age, education (high school or less), financial status (no money to see doctors, own house), retirement rate, and health (smoking). Therefore, the naive analyses of the study variables based on the PRC dataset are subject to selection biases.

We compute the naive and proposed estimators. To apply the proposed method, we assume the sampling score to be a logistic regression model, the continuous outcome to be a linear regression model, and the binary outcome model to be a logistic regression model adjusting for all available covariates. Using cross validation, the double selection procedure identifies 18 important covariates (all available covariates except for the northeast region) in the sampling score and the binary outcome model, and it identifies 15 important covariates (all available covariates except for black, indicator of smoking everyday, the northeast region and the south region).

Table 4 presents the point estimate and the standard error. For estimating the standard error, because the second-order inclusion probabilities are unknown, following the survey literature, we approximate the variance estimator in (32) by assuming the survey design is single-stage Poisson sampling. We find significant differences in the results between our proposed estimator and the corresponding naive estimator. As demonstrated by the simulation in Section 6, the naive estimator may be biased due to selection biases, and the proposed estimator utilizes a probability sample to correct for such biases. From the results, on average, the target population had at least one drink for $4.84$ days over the last month, and $71.8\%$ of the target population voted in local elections.

Acknowledgment

Dr. Yang is partially supported by the National Science Foundation grant DMS 1811245, National Cancer Institute grant P01 CA142538, and Oak Ridge Associated Universities. Dr. Kim is partially supported by the National Science Foundation MMS-1733572. An R package ”IntegrativeFPM” that implements the proposed method is available at https://github.com/shuyang1987/IntegrativeFPM.

SUPPLEMENTARY MATERIAL

Supplementary material provides technical details and proofs.