Robust Bayesian Inference for Big Data: Combining Sensor-based Records with Traditional Survey Data

Denote by $U$ a finite population of size $N<\infty$. We consider the values of a scalar outcome variable, $y_{i}$, $i=1,2,...,N$ and $x_{i}=[1,x_{i1},x_{i2},...,x_{ip}]$, the values of a $p$-dimensional vector of relevant auxiliary variables, $X$. Let $S_{B}$ be a non-probability sample of size $n_{B}$ selected from $U$. The goal is to estimate an unknown finite population quantity, e.g. $Q(Y)$. Here, the quantity of interest is considered to be the finite population mean that is a function of the outcome variable, i.e. $Q(Y)=\overline{y}_{U}=\sum_{i=1}^{N}y_{i}/N$. Suppose $\delta^{B}_{i}=I(i\in S_{B})$ $(i=1,...,N)$ is the inclusion indicator variable of the “big data” survey $S_{B}$ of size $n_{B}$ in $U$. Further, we initially assume that given $X$, elements in $S_{B}$ are independent draws from $U$, but later, we relax this assumption by considering $S_{B}$ to have a single-stage clustering design as is the case in the real-data application of this article.

Suppose $S_{R}$ is a parallel reference survey of size $n_{R}$, for which the same set of covariates, $X$, has been measured. We also define $d_{i}=[1,d_{i1},d_{i2},...,d_{iq}]$, the values of a $q$-dimensional vector of design variables for the reference survey. We assume $y_{i}$ is unobserved for $i\in S_{R}$; otherwise inference could be directly drawn based on $S_{R}$. Also, let $\delta^{R}_{i}=I(i\in S_{R})$ denote the inclusion indicator variable associated with $S_{R}$ for $i\in U$. To avoid unnecessary complexity, we assume that units of $S_{R}$ are selected independently. Being a full probability sample implies that the selection mechanism in $S_{R}$ is ignorable given its design features, i.e. $f(\delta^{R}_{i}|y_{i},d_{i})=f(\delta^{R}_{i}|d_{i})$ for $i\in U$, where $d_{i}$ denotes a $q$-dimensional vector of associated design variables. Thus, one can define the selection probabilities and sampling weights in $S_{R}$ as $\pi^{R}_{i}=p(\delta^{R}_{i}=1|d_{i})$ and $w^{R}_{i}=1/\pi^{R}_{i}$, respectively, for $i\in U$, which we assume are known.

While $X$ and $D$ may overlap or correlate, we define $x^{*}_{i}=[x_{i},d_{i}]^{T}$, the $(p+q)$-dimensional vector of all auxiliary variables associated with $S_{B}$ and $S_{R}$. To be able to make unbiased inference for $S_{B}$, we consider the following assumptions for $S_{B}$:

1. C1.

   Positivity—$S_{B}$ actually does have a probabilistic sampling mechanism, albeit unknown. That means $p(\delta^{B}_{i}=1|x_{i})>0$ for all possible values of $x_{i}$ in $U$.

2. C2.

   Ignorability—the selection mechanism of $S_{B}$ is fully governed by $X$, which implies $Y\rotatebox[origin={c}]{90.0}{$\models$}\delta^{B}|X$. Then, for $i\in U$, the unknown “pseudo-inclusion” probability associated with $S_{B}$ is defined as $\pi^{A}_{i}=p(\delta^{B}_{i}=1|x_{i})$.

3. C3.

   Independence—conditional on $X^{*}$, $S_{R}$ and $S_{B}$ are selected independently, i.e. $\delta^{B}\rotatebox[origin={c}]{90.0}{$\models$}\delta^{R}|X^{*}$.

Note that the first two assumptions are collectively called “strong ignorability” by Rosenbaum and Rubin (1983). Considering C1-C3, the joint density of $y_{i}$, $\delta^{B}_{i}$ and $\delta^{R}_{i}$ can be factorized as below:

$$f(y_{i},\delta^{B}_{i},\delta^{R}_{i}|x^{*}_{i};\theta,\beta,\gamma)=f(y_{i}|x_{i}^{*};\theta)f(\delta^{B}_{i}|x_{i};\beta)f(\delta^{R}_{i}|d_{i};\gamma),\hskip 11.38109pt\forall i\in U$$

(2.1)

where $\Psi=(\theta,\beta,\gamma)$ are some distributional parameters. While $\theta$ and $\beta$ are unknown, $\gamma$ may be known as $S_{R}$ is a probability sample. A QR approach involves modeling $f(\delta^{B}_{i}|x^{*}_{i};\beta)$, whereas a PM approach deals with modeling $f(y_{i}|x^{*}_{i};\theta)$.

Now suppose $S_{B}$ and $S_{R}$ have trivial overlap, i.e. $p(\delta^{B}_{i}+\delta^{R}_{i}=2)\approx 0$. This assumption is reasonable when the sampling fraction in both samples is small. Note that under the ignorable assumption, the propensity model for $S_{B}$ depends on $X$ observed for the entire population. Thus, given the combined sample, $S=S_{B}\cup S_{R}$, with $n=n_{B}+n_{R}$ being the sample size, it is reasonable to expect that the pseudo-inclusion probabilities, $\pi^{B}_{i}$’s, are a function of both $x_{i}$ and $d_{i}$ for $i\in S$. Let $z_{i}=I(i\in S_{B}|\delta_{i}=1)$ be the indicator of subject $i$ belonging to the non-probability sample in the combined sample where $\delta_{i}=\delta^{B}_{i}+\delta^{R}_{i}$. Note that since $S_{B}\cap S_{R}=\emptyset$, $\delta_{i}$ can take values of either $0$ or $1$ as below:

$$\delta_{i}=\begin{cases}0,&\hskip 8.53581pt\text{if}\hskip 8.53581pt\delta^{R}_{i}=0\hskip 8.53581ptand\hskip 8.53581pt\delta^{B}_{i}=0\\ 1,&\hskip 8.53581pt\text{if}\hskip 8.53581pt\delta^{R}_{i}=1\hskip 17.07164ptor\hskip 8.53581pt\delta^{B}_{i}=1\end{cases}$$

Figure 1 illustrates the data structure in both the finite population and the combined sample.

In QR, the non-probability sample is treated as if the self-selection mechanism of population units mimics a stochastic process, but with unknown selection probabilities. Then, attempts are made to estimate these missing quantities in $S_{B}$ based on external information. Conditional on $x^{*}_{i}$, suppose $\pi^{B}_{i}$ follows a logistic regression model in the finite population. We have

$$\pi^{B}(x_{i};\beta)=p(\delta^{B}_{i}=1|x_{i};\beta)=\frac{\exp\{\beta^{T}_{0}x_{i}\}}{1+\exp\{\beta^{T}_{0}x_{i}\}},\hskip 14.22636pt\forall i\in U$$

(2.2)

where $\beta$ is a vector of model parameters in $U$. Using a modified pseudo-maximum likelihood approach (PMLE), Chen et al. (2020) demonstrate that, given $S$, a consistent estimate of $\beta$ can be obtained by solving the following estimating equation with respect to $\beta$:

$$U(\beta)=\sum_{i=1}^{n_{B}}x_{i}-\sum_{i=1}^{n_{R}}\pi^{B}(x_{i};\beta)x_{i}/\pi^{R}_{i}=0$$

(2.3)

The estimates of the $\pi^{B}_{i}$’s, which we also call propensity scores (PS), are obtained by plugging the solution of Eq. 2.3, i.e. $\widehat{\beta}$, into Eq. 2.2. It is important to note that the proposed PS estimator by Chen et al. (2020) depends implicitly on $d_{i}$ in addition to $x_{i}$, because we know that $\pi^{R}_{i}$ is a determinstic function of $d_{i}$ for $i\in U$. Under certain regularity conditions, the authors show that the inverse PS weighted (IPSW) mean from $S_{B}$ yields a consistent and asymptotically unbiased estimate for the population mean.

Obviously, the possible solutions of Eq. 2.3 are not a typical output of logistic regression procedures in the existing standard software. With one additional assumption, which is mutual exclusiveness of the two samples, i.e. $S_{B}\cap S_{R}=\emptyset$, we show that estimating $\pi^{B}_{i}$’s can be reduced to modeling $Z_{i}$ for $i\in S$ instead of modeling $\delta^{B}_{i}$ for $i\in U$. Intuitively, one can view the selection process of the $i$-th population unit in $S_{B}$ as being initially selected in the joint sample ($\delta_{i}=1$) and then being selected in $S_{B}$ given the combined sample ($Z_{i}=1$). By conditioning on $x^{*}_{i}$, the selection probabilities in $S_{B}$ are factorized as

	$\displaystyle p(\delta^{B}_{i}=1|x^{*}_{i})$	$\displaystyle=p(\delta^{B}_{i}=1,\delta_{i}=1|x^{*}_{i})$		(2.4)
		$\displaystyle=p(\delta^{B}_{i}=1|\delta_{i}=1,x^{*}_{i})p(\delta_{i}=1|x^{*}_{i})$
		$\displaystyle=p(Z_{i}=1|x^{*}_{i})p(\delta_{i}=1|x^{*}_{i})\hskip 14.22636pti\in S$

Note that the last expression in Eq. 2.4 results from the definition of $Z_{i}$ given $S$. The same factorization can be derived for the selection probabilities in $S_{R}$. Thus, we have

$\displaystyle p(\delta^{R}_{i}=1|x^{*}_{i})=p(Z_{i}=0|x^{*}_{i})p(\delta_{i}=1|x^{*}_{i})$

(2.5)

By dividing the two sides of the equations 2.4 and 2.5, one can get rid of $p(\delta_{i}=1|x^{*}_{i})$ and obtain the pseudo-selection probabilities in $S_{B}$ as below:

$$p(\delta^{B}_{i}=1|x^{*}_{i})=p(\delta^{R}_{i}=1|x^{*}_{i})\frac{p(Z_{i}=1|x^{*}_{i})}{p(Z_{i}=0|x^{*}_{i})}$$

(2.6)

It is clear that $p(\delta^{R}_{i}=1|x^{*}_{i})=\pi^{R}_{i}$ as $x^{*}_{i}$ contains $d_{i}$ and the sampling design of $S_{R}$ is known given $d_{i}$. As will be seen in Section 2.4, conditioning on $d_{i}$ is vital for the DR estimator, as Chen’s method is limited to situations where the dimension of the auxiliary variables must be the same in QR and PM.

Note that Eq. 2.6 is identical to the pseudo-weighting formula Elliott and Valliant (2017) derive for a non-probability sample. Unlike the PMLE approach, modeling $Z_{i}$ in $S$ can be performed using the standard binary logistic regression or any alternative classification methods, such as supervised machine learning algorithms. Under a logistic regression model, we have

$$p(Z_{i}=1|x^{*}_{i})=\frac{\exp\{\beta^{T}_{1}x_{i}^{*}\}}{1+\exp\{\beta^{T}_{1}x_{i}^{*}\}}$$

(2.7)

where $\beta_{1}$ denotes the vector of model parameters being estimated via maximum likelihood estimation (MLE). Hence, in situations where $\pi^{R}_{i}$ is known or can be calculated for $i\in S_{B}$, the estimate of $\pi^{B}_{i}$ for $i\in S_{B}$ is given by

$$\widehat{\pi}^{B}_{i}=\pi^{R}_{i}\exp\{\widehat{\beta}^{T}_{1}x_{i}^{*}\}=\pi^{R}_{i}\frac{p_{i}(\widehat{\beta}_{1})}{1-p_{i}(\widehat{\beta}_{1})}$$

(2.8)

where $\widehat{\beta}_{1}$ denotes the MLE estimate of the logistic regression model parameters, and $p_{i}(\widehat{\beta}_{1})$ is a shorthand of $p(Z_{i}=1|x^{*}_{i};\widehat{\beta}_{1})$. Intuitively, one can envision that the first factor in 2.8 treats $S_{B}$ as if it is selected under the design of $S_{R}$, and the second factor attempts to balance the distribution of $x$ in $S_{B}$ with respect to that in $S_{R}$.

Having $\pi^{B}_{i}$ estimated based on 2.8 for all $i\in S_{B}$, one can construct the Hájek-type pseudo-weighted estimator for the finite population mean as below:

$$\widehat{\overline{y}}_{PAPW}=\frac{1}{\widehat{N}_{B}}\sum_{i=1}^{n_{B}}\frac{y_{i}}{\widehat{\pi}^{B}_{i}}$$

(2.9)

where $\widehat{N}_{B}=\sum_{i=1}^{n_{B}}1/\pi^{B}_{i}$. Hereafter, we refer to the estimator in 2.9 as propensity-adjusted probability weighting (PAPW). Under mild regularity conditions, the ignorable assumption in $S_{B}$ given $x$, the logistic regression model and the additional assumption of $S_{B}\cap S_{R}=\emptyset$, Appendix 8.1 shows that this estimator is consistent and asymptotically unbiased for $\overline{y}_{U}$. Further, when $\pi_{i}^{R}$ is known, the sandwich-type variance estimator for $\widehat{\overline{y}}_{PAPW}$ is given by

	$\displaystyle\widehat{Var}\left(\widehat{\overline{y}}_{PAPW}\right)=\frac{1}{N^{2}}\sum_{i=1}^{n_{B}}\big{\{}1-\widehat{\pi}^{B}_{i}\big{\}}\left(\frac{y_{i}-\widehat{\overline{y}}_{PAPW}}{\widehat{\pi}^{B}_{i}}\right)^{2}$	$\displaystyle-2\frac{\widehat{b}^{T}}{N^{2}}\sum_{i=1}^{n_{B}}\big{\{}1-p_{i}(\widehat{\beta}_{1})\big{\}}\left(\frac{y_{i}-\widehat{\overline{y}}_{PAPW}}{\widehat{\pi}^{B}_{i}}\right)x^{*}_{i}$		(2.10)
		$\displaystyle+\widehat{b}^{T}\left[\frac{1}{N^{2}}\sum_{i=1}^{n}p_{i}(\widehat{\beta}_{1})x^{*}_{i}x_{i}^{*T}\right]\widehat{b}$

where

$$\widehat{b}^{T}=\bigg{\{}\frac{1}{N}\sum_{i=1}^{n_{B}}\left(\frac{y_{i}-\widehat{\overline{y}}_{PAPW}}{\widehat{\pi}^{B}_{i}}\right)x^{*T}_{i}\bigg{\}}\bigg{\{}\frac{1}{N}\sum_{i=1}^{n}p_{i}(\widehat{\beta}_{1})x^{*}_{i}x_{i}^{*T}\bigg{\}}^{-1}$$

(2.11)

where $\widehat{\pi}^{B}_{i}$ is the estimated pseudo-selection probability based on 2.9 for $i\in S_{B}$. See Appendix 8.1 for the derivation.

In situations where $\pi^{R}_{i}$ is unknown for $i\in S_{B}$, Elliott and Valliant (2017) suggest predicting this quantity for units of the non-probability sample. Note that, in this situation, it is no longer required to condition on $d_{i}$ in addition to $x_{i}$. Treating $\pi^{R}_{i}$ as a random variable for $i\in S_{B}$ conditional on $x_{i}$, one can obtain this quantity by regressing the $\pi^{R}_{i}$’s on the $x_{i}$’s in the reference survey. We have

	$\displaystyle p(\delta^{R}_{i}=1|x_{i})$	$\displaystyle=\int_{0}^{1}p(\delta^{R}_{i}=1|\pi^{R}_{i},x_{i})p(\pi^{R}_{i}|x_{i})d\pi^{R}_{i}$		(2.12)
		$\displaystyle=\int_{0}^{1}\pi^{R}_{i}p(\pi^{R}_{i}|x_{i})d\pi^{R}_{i}$
		$\displaystyle=E(\pi^{R}_{i}|x_{i})\hskip 14.22636pti\in S_{R}.$

However, since the outcome is continuous bounded taking values within $(0,1)$, fitting a $Beta$ regression model is recommended (Ferrari and Cribari-Neto, 2004). Note that, $\pi^{R}_{i}$ is fixed given $d_{i}$ as $S_{R}$ is a probability sample, but conditional on $x_{i}$, $\pi^{R}_{i}$ can be regarded as a random variable.

Rafei et al. (2020) call this approach propensity-adjusted probability prediction (PAPP). This two-step derivation of pseudo-inclusion probabilities is especially useful, as it separates sampling weights in $S_{R}$ from the propensity model computationally. When the true model is unknown, this feature enables us to fit a broader and more flexible range of models, such as algorithmic tree-based methods. It is worth noting that modeling $E(\pi^{R}_{i}|x_{i})$ does not impose an additional ignorable assumption in $S_{R}$ given $x$, because in the extreme case if $\delta^{R}_{i}\rotatebox[origin={c}]{90.0}{$\models$}x_{i}$, that means weighted and unweighted distributions of $x$ are identical in $S_{R}$, and therefore the $\pi^{R}_{i}$’s can be safely ignored in propensity modeling.

An alternative approach to deal with selectivity in Big Data is modeling $f(y|x^{*})$ (Smith, 1983). In a fully model-based fashion, this essentially involves imputing $y$ for the population non-sampled units, $U-S_{B}$. When $x^{*}$ is unobserved for non-sampled units, it is recommended that a synthetic population is generated by undoing the selection mechanism of $S_{R}$ through a non-parametric Bayesian bootstrap method using the design variables in $S_{R}$ (Dong et al., 2014; Zangeneh and Little, 2015). In the non-probability sample context, Elliott and Valliant (2017) propose an extension of the General Regression Estimator (GREG) when only summary information about $x^{*}$, such as totals, is known regarding $U$. In situations where an external probability sample is available with $x^{*}$ measured, an alternative is to limit the outcome prediction to the units in $S_{R}$, and then, use design-based approaches to estimate the population quantity (Rivers, 2007; Kim et al., 2018).

However, to the best of our knowledge, none of the prior literature distinguish the role of $D$ from $X$ in the conditional mean structure of the outcome, while the likelihood factorization in Eq. 2.1 indicates that predicting $y$ requires conditioning not only on $x$ but also on $d$. Suppose $U$ is a realization of a repeated random sampling process from a super-population under the following model:

$\displaystyle E(y_{i}|x^{*}_{i};\theta)=m(x^{*}_{i};\theta)\hskip 14.22636pt\forall i\in U$

(2.13)

where $m(x^{*}_{i};\theta)$ can be either a parametric model with $m$ being a continuous differentiable function or an unspecified non-parametric form. Under the ignorable condition where

$\displaystyle f(y_{i}|x^{*}_{i},z_{i}=1;\theta)=f(y_{i}|x^{*}_{i},z_{i}=0;\theta)=f(y_{i}|x_{i},d_{i};\theta)$

(2.14)

an MLE estimate of $\theta$ can be obtained by regressing $Y$ on $X^{*}$ given $S_{B}$. The predictions for units in $S_{R}$ are then given by

$\displaystyle\widehat{y}_{i}=E(y_{i}|x^{*}_{i},z_{i}=0;\widehat{\theta})=m(x^{*}_{i};\widehat{\theta})\hskip 14.22636pt\forall i\in S_{R}$

(2.15)

where $m(x^{*}_{i};\widehat{\theta})=\widehat{\theta}^{T}x^{*}_{i}$. Once $y$ is imputed for all units in the reference survey, the population mean can be estimated through the Hájek formula as below:

$\displaystyle\widehat{\overline{y}}_{PM}=\frac{1}{\widehat{N}_{R}}\sum_{i=1}^{n_{R}}\frac{\widehat{y}_{i}}{\pi^{R}_{i}}$

(2.16)

where $\widehat{y}_{i}=m(x^{*}_{i};\widehat{\theta})$ for $i\in S_{R}$, $\widehat{N}_{R}=\sum_{i=1}^{n_{R}}w^{R}_{i}$ and $\pi^{R}_{i}$ is the selection probability for subject $i\in S$.

The asymptotic properties of the estimator in 2.16, including consistency and unbiasedness, have been investigated by Kim et al. (2018). Note that in situations where $\pi^{R}_{i}$ is available for $i\in S_{B}$, one can use $w^{R}_{i}$ instead of the high-dimensional $d_{i}$ as a predictor in $m(.)$. This method is known as linear-in-the-weight prediction (LWP) (Scharfstein et al., 1999; Bang and Robins, 2005; Zhang and Little, 2011). However, since outcome imputation relies fully on extrapolation, even minor misspecification of the underlying model can be seriously detrimental to bias correction.

To reduce the sensitivity to model misspecification, Chen et al. (2020) reconcile the two aforementioned approaches, i.e. QR and PM, in a way that estimates remain consistent even if one of the two models is incorrectly specified. Their method involves an extension of the augmented inverse propensity weighting (AIPW) proposed by Robins et al. (1994). When $N$ is known, the expanded AIPW estimator takes the following form:

$$\overline{y}_{DR}=\frac{1}{N}\sum_{i=1}^{n_{B}}\frac{\{y_{i}-m(x^{*}_{i};\theta)\}}{\pi^{B}(x^{*}_{i};\beta)}+\frac{1}{N}\sum_{j=1}^{n_{R}}\frac{m(x^{*}_{i};\theta)}{\pi^{R}_{j}}$$

(2.17)

where given $x^{*}$, $\beta$ and $\theta$ are estimated using the modified PMLE and MLE methods mentioned in sections 2.2 and 2.3, respectively. The theoretical proof of the asymptotic unbiasedness of $\overline{y}_{DR}$ under the correct modeling of $\pi^{B}(x^{*}_{i};\beta)$ or $m(x^{*}_{i};\theta)$ is reviewed in Appendix 8.1.

To avoid using $\pi^{R}$ in modeling $\delta^{B}_{i}$ because of the PMLE restrictions we discussed in Section 2.2, in this study, we suggest estimating $\pi^{B}_{i}$ for $i\in S_{B}$ in Eq. 2.17 based on the PAPW/PAPP method depending on whether $\pi^{R}_{i}$ is available for $i\in S_{B}$ or not. As a result, in situations where $\pi^{R}_{i}$ is known for $i\in S_{B}$, our proposed DR estimator is given by

$$\widehat{\overline{y}}_{DR}=\frac{1}{N}\sum_{i=1}^{n_{B}}\frac{1}{\pi^{R}_{i}}\left[\frac{1-p_{i}(\beta_{1})}{p_{i}(\beta_{1})}\right]\{y_{i}-m(x^{*}_{i};\theta)\}+\frac{1}{N}\sum_{j=1}^{n_{R}}\frac{m(x^{*}_{i};\theta)}{\pi^{R}_{j}}$$

(2.18)

wher $\pi^{B}(x^{*}_{i};\beta)$ is substituted using Eq. 2.8. We demonstrate that this form of the AIPW estimator is identical to that defined by Kim and Haziza (2014) in the non-response adjustment context under probability surveys. Assuming that $y_{i}$ is fully observed for $i\in S_{R}$, let us define the following HT-estimator for the population mean:

$$\widehat{\overline{y}}_{U}=\frac{1}{N}\sum_{i=1}^{n_{R}}\frac{y_{i}}{\pi^{R}_{i}}$$

(2.19)

Now, one can easily conclude that

	$\displaystyle\widehat{\overline{y}}_{DR}$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{n}\frac{1}{\pi^{R}_{i}}\left[Z_{i}\left(\frac{1-p_{i}(\beta_{1})}{p_{i}(\beta_{1})}\right)\{y_{i}-m(x^{*}_{i};\theta)\}+(1-Z_{i})m(x^{*}_{i};\theta)\right]$		(2.20)
		$\displaystyle=\widehat{\overline{y}}_{U}+\frac{1}{N}\sum_{i=1}^{n}\frac{1}{\pi^{R}_{i}}\left[\frac{Z_{i}}{p_{i}(\beta_{1})}-1\right]\big{\{}y_{i}-m(x^{*}_{i};\theta)\big{\}}$

where $p_{i}(\beta_{1})=p(Z_{i}=1|x_{i}^{*};\beta_{1})$. The formula in 2.20 is similar to what is derived by Kim and Haziza (2014). Therefore, the rest of the theoretical proof of asymptotic unbiasedness, i.e. $\widehat{\overline{y}}_{DR}-\overline{\widehat{y}}_{U}=O_{p}(n_{B}^{-1/2})$, in Kim and Haziza (2014) should hold for our modified AIPW estimator in 2.18 as well.

To preserve the DR property for both the point and variance estimator of $\overline{y}_{DR}$, as suggested by Kim and Haziza (2014), one can solve the following estimating equations simultaneously given $S$ to obtain the estimate of $(\beta_{1},\theta)$. The aim is to cancel the first-order derivative terms in the Taylor-series expansion of $\widehat{\overline{y}}_{DR}-\widehat{\overline{y}}_{U}$ under QR and PM. These estimating equations are given by

	$\displaystyle\frac{\partial}{\partial\beta_{1}}\left[\widehat{\overline{y}}_{DR}-\widehat{\overline{y}}_{U}\right]$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{n}\frac{Z_{i}}{\pi^{R}_{i}}\left[\frac{1}{p_{i}(\beta_{1})}-1\right]\{y_{i}-m(x^{*}_{i};\theta)\}x^{*}_{i}=0$		(2.21)
	$\displaystyle\frac{\partial}{\partial\theta}\left[\widehat{\overline{y}}_{DR}-\widehat{\overline{y}}_{U}\right]$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{n}\frac{1}{\pi^{R}_{i}}\left[\frac{Z_{i}}{p_{i}(\beta_{1})}-1\right]\dot{m}(x^{*}_{i};\theta)=0$

where $\dot{m}$ is the derivative of $m$ with respect to $\beta_{1}$. Under a linear regression model, $\dot{m}(x^{*}_{i})=x^{*}_{i}$. Therefore, given the same regularity conditions, ignorability in $S_{B}$, the logistic regression model as well as the additional imposed assumption of $S_{B}\cap S_{R}=\emptyset$, one can show that the proposed DR estimator is consistent and approximately unbiased given that either the QR or PM model holds.

It is important to note that the system of equations in 2.21 may not have unique solutions unless the dimension of covariates in QR and PM is identical. Therefore, the AIPW estimator by Chen et al. (2020) may not be applicable here, as our likelihood factorization suggests that conditioning on $d_{i}$ is necessary at least for the PM. Furthermore, when $\pi^{R}_{i}$ is known for $i\in S_{B}$, one can replace the $q$-dimensional $d_{i}$ with the $1$-dimensional $w^{R}_{i}$ in modeling both QR and PM. Bang and Robins (2005) shows that estimators based on a linear-in-weight prediction model remains consistent.

A fully Bayesian approach specifies a model for the joint distribution of selection indicator, $\delta^{B}_{i}$, and the outcome variable, $y_{i}$, for $i\in U$ (McCandless et al., 2009; An, 2010). This requires multiply generating synthetic populations and fitting the QR and PM models on each of them repeatedly (Little and Zheng, 2007; Zangeneh and Little, 2015), which can be computationally expensive under a Big Data setting. While joint modeling may result in good frequentist properties (Little, 2004), feedback occurs between the two models (Zigler et al., 2013). This can be controversial in the sense that PS estimates should not be informed by the outcome model (Rubin, 2007). Here, we are interested in modeling the PS and the outcome separately through the two-step framework proposed by Kaplan and Chen (2012). The first step involves fitting Bayesian models to multiply impute the PS and the outcome by randomly subsampling the posterior predictive draws, and Rubin’s combining rules are utilized as the second step to obtain the final point and interval estimates. This method not only is computationally efficient as it suffices to fit the models once and on the combined sample but also cuts the undesirable feedback between the models as they are fitted separately. Bayesian modeling can be performed either parametrically or non-parametrically.

To obtain an unbiased variance estimate for the the proposed DR estimator, one needs to account for three sources of uncertainty: (i) the uncertainty due to estimated pseudo-weights in $S_{B}$, (ii) the uncertainty due to the predicted outcome in both $S_{B}$ and $S_{R}$, and (iii) the uncertainty due to sampling itself. In the following, we consider two scenarios: