Comments on “Statistical inference with non-probability survey samples”
Section 2. Additional approaches to combining data from probability and non-probability surveys

• Previous
• Next

Dr. Wu’s paper follows the general prescription of 1) using model estimation and subsequent calibration to probability-sample-estimated covariate distributions, 2) developing propensity score estimates based on discrepancies between the probability- and non-probability sample data, and 3) doubly-robust methods that combine 1) and 2) in a manner such that only one of the two underlying models needs to be correct.

2.1  Propensity score estimators

Rivers (2007) appears to have been the first to suggest estimating propensity score using logistic regression with membership in the non-probability sample as the outcome and taking the reciprocal of the resulting propensity scores to use as inclusion weights. This approach was formalized further in Valliant and Dever (2011). Separately, using simple results from Bayes’ theorem and discriminant analysis first described in Elliott and Davis (2005), Elliott, Resler, Flannagan and Rupp (2010) and Elliott (2013) developed a somewhat different estimator of the form

$${\widehat{\pi }}_i^A\left(x_i\mathrm{,}\alpha \right)\mathrm{<mo>=</mo>}\widehat{P}\left(i\in S_A\right)\propto P\left(i\in S_B\right)\frac{\widehat{P}\left(i\in S_A|i\in S_A\text{or}i\in S_B\mathrm{,}{x}_i\mathrm{,}\alpha \right)}{\widehat{P}\left(i\in S_B|i\in S_A\text{or}i\in S_B\mathrm{,}{x}_i\mathrm{,}\alpha \right)}\mathrm{.}(2.1)$$

$\widehat{P}\left(i\in S_A|i\in S_A\text{or}i\in S_B\mathrm{,}{x}_i\mathrm{,}\alpha \right)$ can be obtained using logistic regression, or using one of the suite of machine learning-type approaches such as support vector machines (Soentpiet, 1999), targeted maximum likelihood estimation (Van Der Laan and Rubin, 2006), or Bayesian Additive Regression Trees (BART) (Chipman, George and McCulloch, 2010), and $\widehat{P}\left(i\in S_A|i\in S_B\text{or}i\in S_B\mathrm{,}{x}_i\mathrm{,}\alpha \right)$ obtained as $1-\widehat{P}\left(i\in S_A|i\in S_A\text{or}i\in S_B\mathrm{,}{x}_i\mathrm{,}\alpha \right).$ In principle $P\left(i\in S_B\right)\mathrm{<mo>=</mo>}1/d_i^B$ is known since sampling probabilities are known for all elements of the population, including those in the non-probability sample, but in practice analysts with access only to public use data may have to estimate this as well. (In addition, $d_i^B$ may include calibration and non-response adjustments that are not known for the non-probability sample elements.) This last point is critical as use of the probability sample to develop propensity scores using only the discrepancies between the non-probability sample and the probability sample will be biased unless the probability sample used an equal probability (epsem) design, as noted by Wu.

In contrast, Chen, Li and Wu (2020) shows that using a pseudo-likelihood approach to estimating ${\widehat{\pi }}_i^A\left(x_i\mathrm{,}\alpha \right)$ directly from the population likelihood for the indicators $I\left(i\in S_A\right)$ as a function of $x_i$ yields an estimator that does not require $P\left(i\in S_B\right)$ for elements in the non-probability sample under the restriction that ${\pi }_i^A\left(x_i\mathrm{,}\alpha \right)$ follows a generalized linear model with a canonical link, i.e., logistic regression.

(None of these approaches actually has the correct intercept to obtain a true propensity score; however, as noted in Wu, weighted estimation usually uses Hájek-type estimators [using weights to estimate a population total for denominators; Hájek, 1971] so that propensity scores estimated up to a normalizing constant are sufficient.)

2.2  Doubly-robust estimators

If inference is focused on a particular variable $Y$ available only in the non-probability sample, we can return to the model-assisted estimators that date back to Cassel, Särndal and Wretman (1976), which posit a model for the expectation $E\left(y_i|x_i\right)\mathrm{<mo>=</mo>}{m}_i.$ Combining this with propensity score estimates of the probability of being in the non-probability sample (which we are treating as an “unknown probability sample” $–$ more about this under Assumptions below) yields estimators of the form

$$\frac{1}{{\widehat{N}}^A}{\displaystyle \sum _{i\in S_A}}\frac{y_i-{\widehat{m}}_i}{{\widehat{\pi }}_i^A}+\frac{1}{{\widehat{N}}^B}{\displaystyle \sum _{i\in S_B}}{d}_i^B{\widehat{m}}_i(2.2)$$

corresponding to ${\widehat{\mu }}_{\text{DR}2}$ of (4.11) in Wu. The intuition is that any bias due to the model misspecification in estimation of $m_i$ in $\frac{1}{{\widehat{N}}^B}{\displaystyle {\sum }_{i\in S_B}}{d}_i^B{\widehat{m}}_i$ will be equal to and opposite in sign of $\frac{1}{{\widehat{N}}^A}{\displaystyle {\sum }_{i\in S_A}}\frac{y_i-{\widehat{m}}_i}{{\widehat{\pi }}_i^A}$ if the model for ${\pi }_i^A$ is correctly specified. Conversely, if the model for ${\pi }_i^A$ is misspecified but $m_i$ is correctly specified, $y_i-{\widehat{m}}_i$ will be iid with mean zero and consequently $\frac{1}{{\widehat{N}}^A}{\displaystyle {\sum }_{i\in S_A}}\frac{y_i-{\widehat{m}}_i}{{\widehat{\pi }}_i^A}$ will also have mean 0, yielding an unbiased estimator. Chen, Valliant and Elliott (2019) used LASSO for prediction in combination with generalized regression estimators (McConville, Breidt, Lee and Moisen, 2017) when $X$ is of high dimension. As Wu notes, Wu and Sitter (2001) show the equivalence between GREG applied to predicted values and DR estimators of the form in (2.2), which indicates that the Chen et al. (2019) approach was equivalent to (2.2) with LASSO estimation for $m_i$ and an assumption of simple random sampling for the non-probability sample.

A disadvantage of using (2.1) as opposed to Chen et al. (2020) as the estimator of ${\pi }_i^A,$ and thus of $d_i^A,$ is the requirement that the probability sample weights $d_i^B$ be known or at least estimated for the non-probability sample. An advantage of using (2.1), is that non-linear models and machine learning methods can be used in estimation. Rafei, Flannagan and Elliott (2020) uses BART to estimate both $m_i$ and ${\pi }_i^A,$ reducing the impact of potential model misspecification. Simulations showed considerable improvement in bias and variance reduction over the method of Chen et al. (2020) when the linear models is misspecified. Variance estimation can proceed by adapting Rubin’s multiple imputation rules: from $M$ independent draws from BART, the mean of the variances computed treating the draw of $d_i^A$ as known using standard complex sample design estimators and added to $\frac{M+1}{M}$ times the variance of the point estimates computed across the draws of $d_i^A$ yield an approximately unbiased variance estimator.

An alternative approach to doubly-robust estimation uses the fact that the propensity score is the coarsest possible “balancing score” that contains all of the information about the association between the sampling indicator and the outcome of interest. This has led to the development of mean estimators that use smooth functions of weights to produce consistent estimators that can be more efficient when weights are highly variable or only weakly related to the outcome (Elliott and Little, 2000; Zheng and Little, 2005). Zhou, Elliott and Little (2019) extended this idea into the causal inference setting in non-randomized settings, in which probability of assignment to a treatment or exposure (propensity score) is estimated as a function of covariates $P_Z\left(x_i\mathrm{,}\alpha \right)$ using logistic regression, and then non-observed potential outcomes $Y^z$ under treatment arm $z_i^{{}^{\prime }}\ne z_i$ for observed treatment $z_i$ are imputed from

where $P^{*}$ is the logit transformation of $P,$ $s({\widehat{P}}_Z^{*}|{\theta }_Z)$  denotes a penalized spline with fixed knots (Eilers and Marx, 1996) of propensity, and $g_Z({\widehat{P}}^{*}\mathrm{,}{x}_i|{\beta }_Z)$ is a general function of covariates including the propensity scores. The resulting estimator is doubly robust in the sense that if either $P_Z\left(x_i\mathrm{,}\alpha \right)$ or $E\mathrm{<mo\; stretchy="false">(</mo>}{Y}^z\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{g}_Z\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{P}}^{*}\text{}\mathrm{,}{x}_i|{\beta }_Z\mathrm{<mo\; stretchy="false">)</mo>}$ is correctly specified, $Y^{\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{<mo\; stretchy="false">)</mo>}}$ will be approximately unbiased; see Zhang and Little (2009). This can be implemented in the non-probability setting by replacing ${\widehat{P}}_Z\mathrm{<mo\; stretchy="false">(</mo>}{x}_i\mathrm{,}\alpha \mathrm{<mo\; stretchy="false">)</mo>}$ in the mean model for (2.3) with ${\widehat{\pi }}_i^A$ estimated using (2.1) to obtain a draw of $Y_i^{\mathrm{<mo\; stretchy="false">(</mo>}b\mathrm{<mo\; stretchy="false">)</mo>}}.$ (Note this requires obtaining ${\widehat{\pi }}_i^A$ for the probability sample elements requiring prediction.) Inference can proceed by obtaining $b\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,B$ draws from the posterior distribution of the estimated population quantity of interest, e.g., for the population mean

$$Y^{\mathrm{<mo\; stretchy="false">(</mo>}b\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\frac{{\displaystyle {\sum }_{i\in S_R}{N}_i^{\mathrm{<mo\; stretchy="false">(</mo>}b\mathrm{<mo\; stretchy="false">)</mo>}}{Y}_i^{\mathrm{<mo\; stretchy="false">(</mo>}b\mathrm{<mo\; stretchy="false">)</mo>}}}+{\displaystyle {\sum }_{i\in S_A}\left(y_i-Y_i^{\mathrm{<mo\; stretchy="false">(</mo>}b\mathrm{<mo\; stretchy="false">)</mo>}}\right)}}{N}$$

where now $N_i^{\mathrm{<mo\; stretchy="false">(</mo>}b\mathrm{<mo\; stretchy="false">)</mo>}}$ is a estimate of the population represented by the weight $d_i^R$ obtained from a finite population Bayesian bootstrap (Little and Zheng, 2007); more complete FBPP extensions to complex sample designs that include clustering and stratification are available in Dong, Elliott and Raghunathan (2014).

As in the estimation of (2.1), the non-parametric (spline) component of (2.3) can be replaced with other machine-learning estimators; see Chapter 4 of Rafei (2021) for implementation using Gaussian processes. Also, extensions to non-normal models are direct, although not necessarily computational easy.

2.3  Poststratified estimators

Wu also describes the use of poststratified estimators in the context of quota sampling, which is not only a very old form of non-probability sampling but indeed the standard before Neyman made the case for stratified random sampling (Neyman, 1934). Wu’s Section 5 suggests a robust alternative to the propensity score estimates obtained by ordering observations in the probability sample by ${\widehat{\pi }}_i,$ stratifying into $K$ strata based on this ordering, and computing the predicted proportion of the population belonging to the $k^{\text{th}}$ stratum as proportion of the sample weights $W_k$ in this stratum using the probability sample, with

$${\widehat{\mu }}_{\text{PST}}\mathrm{<mo>=</mo>}{\displaystyle \sum _k}{\widehat{W}}_k{\overline{y}}_k(2.4)$$

where ${\overline{y}}_k$ is the mean within the $k^{\text{th}}$ stratum in the non-probability sample. Wu notes the tradeoff between choosing $K$ to be large enough to retain homogeneity within units but small enough to obtain stable estimates of ${\overline{y}}_k,$ suggesting 30 as the old “rule of thumb” for “large [enough] sample sizes”. I would add that a more formal approach discussed in Little (1986) suggests a method to generate strata (there in the context of non-response adjustment) that minimizes mean square error by maximizing the between-stratum-to-within-stratum variance. It would seem such an approach would be appropriate to consider in the non-probability post-stratified estimator as well.

A more direct approach to obtain estimates using a post-stratified type estimator is multilevel regression and poststratification (Wang, Rothschild, Goel and Gelman, 2015; Downes and Carlin, 2020). Here only data from the non-probability sample is used in the outcome model:

$$E\left(Y_{k\mathrm{<mo\; stretchy="false">[</mo>}i\mathrm{<mo\; stretchy="false">]</mo>}}\right)\mathrm{<mo>=</mo>}{\beta }_0+x_k^T\beta +{\displaystyle \sum _j}{a}_{l\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}}^j(2.5)$$

where $k\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,K$ indexes the poststratum developed from $j\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,J$ variables, $a_{l\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}}^j~N\mathrm{<mo\; stretchy="false">(</mo>0,}{\sigma }_j^2\mathrm{<mo\; stretchy="false">)</mo>}$ for $l\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,L_j$ and $l\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}$ maps the postratum cell $k$ to the appropriate category $l$ of variable $j.$ The poststratifed estimator is still given by (2.4) with ${\widehat{W}}_k$ now replaced with known population totals $W_k;$ posterior inference is obtained though posterior draws of ${\beta }_0,$ $\beta ,$ and $a_{l\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}}^j$ to obtain a draw of

Though not technically doubly-robust, it has been shown to work well in some applications where $J$ is large enough to capture all of the important discrepancies between the probability and non-probability sample, and the non-probability sample is sufficiently large to allow reasonably accurate estimation of $a_{l\mathrm{<mo\; stretchy="false">[</mo>}k\mathrm{<mo\; stretchy="false">]</mo>}}^j.$ In the absence of known joint distributions of a high dimensional $X,$ this approach has the weakness of relying on estimated distributions, which are unstable. A possible alternative might be replace the simple ${\overline{y}}_k$ with (2.5) in Wu’s poststratified estimator (2.4), using the fact that the sampling weights $d_i^R$ summarize the information about $X$ in the probability sample similar to that of the propensity score for non-probability sample.

• Previous
• Article main page
• Next