Statistical inference with non-probability survey samples
Section 6. Variance estimation

• Previous
• Next

Variance estimation under the two sample $S_A$ and $S_B$ setup involves at least two different sources of variation. The probability sampling design for the reference sample $S_B$ remains one of the sources regardless of the approaches used for non-probability survey samples. Estimation of the variance component due to the use of $S_B$ requires either suitable variance approximation formulas or replication weights as part of the dataset from the reference probability sample. Our discussion in this section assumes that a design-based variance estimator for the survey weighted point estimator based on $S_B$ is available.

6.1  Variance estimation for mass imputation estimators

Variance estimation for the model-based prediction estimator ${\widehat{\mu }}_y$ involves first deriving the asymptotic variance formula for $\text{Var}\left({\widehat{\mu }}_y-{\mu }_y\right)$ under the assumed outcome regression model or the imputation model $\xi$ and the probability sampling design $p,$ and then using plug-in estimators for various unknown population quantities.

The mass imputation estimator ${\widehat{\mu }}_{y\text{MI}}={\widehat{N}}_B^{-1}{\displaystyle {\sum }_{i\in S_B}{d}_i^B{y}_i^{*}}$ given in (3.5) is a special type of model-based prediction estimator, where the model $\xi$ refers to the one used for imputation and is not necessarily the same as the outcome regression model. The imputation method plays a key role in deriving the asymptotic variance formula, and the variance estimator needs to be constructed accordingly. Noting that ${\widehat{\mu }}_{y\text{MI}}$ is a Hájek type estimator due to the use of the estimated population size ${\widehat{N}}_B,$ derivations of the asymptotic variance formula start with putting the true value $N$ in first and then dealing with ${\widehat{\mu }}_{y\text{MI}}$ as a ratio estimator. Kim et al. (2021) considered variance estimation for ${\widehat{\mu }}_y=N^{-1}{\displaystyle {\sum }_{i\in S_B}{d}_i^B{y}_i^{*}},$ where $y_i^{*}=m\left(x_i\mathrm{,}\widehat{\beta }\right)$ is the imputed value for $y_i$ based on the semiparametric model (3.1). The asymptotic variance formula is developed in two steps. First, a linearized version of ${\widehat{\mu }}_y$ is obtained by using a Taylor series expansion at ${\beta }^{*},$ where ${\beta }^{*}$ is the probability limit of $\widehat{\beta }$ such that $\widehat{\beta }={\beta }^{*}+O_p\left(n_A^{-1/2}\right).$ Second, two variance components are derived for $\text{Var}\left({\widehat{\mu }}_y-{\mu }_y\right)$ based on the linearized version using the semiparametric model (3.1) and the sampling design for $S_B.$ The process is tedious, which is the case for most model-based variance estimation methods. A bootstrap variance estimator turns out to be more attractive for practical applications. See Kim et al. (2021) for further details.

6.2  Variance estimation for IPW estimators

The commonly used IPW estimator ${\widehat{\mu }}_{\text{IPW}2}$ given in (4.8) is valid under the assumed model $q$ for the propensity scores. An explicit asymptotic variance formula for ${\widehat{\mu }}_{\text{IPW}2}$ can be derived under the joint $qp$ -framework when the propensity scores are estimated using the pseudo maximum likelihood method or an estimating equation based method as discussed in Section 4.1. The theoretical tool is the sandwich-type variance formula for point estimators defined as the solution to a combined system of estimating equations for both ${\mu }_y$ and ${\alpha }_0.$

Consider the parametric form ${\pi }_i^A=\pi \left(x_i\mathrm{,}\alpha \right)$ for the propensity scores, where the model parameters $\alpha$ are estimated through the estimating equations (4.4) with user-specified functions $h(x,\alpha ).$ The first major step in deriving the asymptotic variance formula for ${\widehat{\mu }}_{\text{IPW}2}$ is to write down the system of joint estimating equations for both ${\mu }_y$ and ${\alpha }_0.$ Let $\eta ={\left(\mu \mathrm{,}{\alpha }^{\prime }\right)}^{\prime }$ be the vector of the combined parameters. The estimator $\widehat{\eta }={\left({\widehat{\mu }}_{\text{IPW}2}\mathrm{,}{\alpha }^{\prime }\right)}^{\prime }$ is the solution to the system of joint estimating equations ${\Phi }_n\left(\eta \right)=0,$ where

$${\Phi }_n\left(\eta \right)\mathrm{<mo>=</mo>}\left(\begin{array}{c}{N}^{-1}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^N{R}_i\left(y_i-\mu \right)/{\pi }_i^A}\\ N^{-1}{\displaystyle {\sum }_{i=1}^N{R}_{i}h\left(x_i\mathrm{,}\alpha \right)}-N^{-1}{\displaystyle {\sum }_{i\in S_B}{d}_i^B{\pi }_i^{A}h\left(x_i\mathrm{,}\alpha \right)}\end{array}\right).(6.1)$$

The factor $N^{-1}$ is redundant but useful in facilitating asymptotic orders. The estimating functions defined by (6.1) are unbiased under the joint $qp$ -framework, i.e., $E_{qp}\left\{\Phi ({\eta }_0)\right\}=0,$ where ${\eta }_0={\left({\mu }_y\mathrm{,}{\alpha }_0^{\prime }\right)}^{\prime }.$ There are two major consequences from the unbiasedness of the estimating equations system. First, consistency of the estimator $\widehat{\eta }$ can be argued using the theory of general estimating functions similar to those presented in Section 3.2 of Tsiatis (2006). Second, the asymptotic variance-covariance matrix of $\widehat{\eta },$ denoted as $\text{AV}\left(\widehat{\eta }\right),$ has the standard sandwich form and is given by

$$\text{AV}\left(\widehat{\eta }\right)={\left[E\left\{{\varphi }_n({\eta }_0)\right\}\right]}^{-1}\text{Var}\left\{{\Phi }_n({\eta }_0)\right\}{\left[E{\left\{{\varphi }_n({\eta }_0)\right\}}^{\prime }\right]}^{-1},$$

where ${\varphi }_n(\eta )=\partial {\Phi }_n(\eta )/\partial \eta ,$ which depends on the forms of ${\pi }_i^A=\pi (x_i\mathrm{,}\alpha )$ and $h(x_i\mathrm{,}\alpha ).$ The term $\text{Var}\left\{{\Phi }_n({\eta }_0)\right\}$ consists of two components, one due to the propensity score model $q$ and the other from the probability sampling design for $S_B.$ More specifically, we have $\text{Var}\left\{{\Phi }_n({\eta }_0)\right\}\mathrm{<mo>=</mo>}{V}_q(A_1)+V_p(A_2),$ where $V_q(\cdot )$ denotes the variance under the propensity score model $q$ and $V_p(\cdot )$ represents the design-based variance under the probability sampling design $p,$ and

$$A_1=\frac{1}{N}{\displaystyle \sum _{i=1}^N{R}_i\left(\begin{array}{c}(y_i-\mu )/{\pi }_i^A\\ h(x_i\mathrm{,}\alpha )\end{array}\right)}\mathrm{,}{A}_2\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^B\left(\begin{array}{c}0\\ {\pi }_i^{A}h(x_i\mathrm{,}\alpha )\end{array}\right)}\mathrm{.}$$

The analytic expression for $V_q(A_1)$ follows immediately from $V_q(R_i)={\pi }_i^A(1-{\pi }_i^A)$ and the independence among $R_1\mathrm{,}\dots \mathrm{,}{R}_N.$ The design-based variance component $V_p(A_2)$ requires additional information on the survey design for $S_B$ or a suitable variance approximation formula with the given design.

The asymptotic variance formula for the IPW estimator ${\widehat{\mu }}_{\text{IPW}2}$ is the first diagonal element of the matrix $\text{AV}\text{(}\widehat{\eta }\text{)}\text{.}$ The final variance estimator for ${\widehat{\mu }}_{\text{IPW}2}$ can then be obtained by replacing various population quantities with sample-based moment estimators. Chen et al. (2020) presented the variance estimator with explicit expressions when ${\pi }_i^A=\pi (x_i\mathrm{,}\alpha )$ are modelled by the logistic regression and the $\widehat{\alpha }$ is obtained by the pseudo maximum likelihood method.

6.3  Variance estimation for doubly robust estimators

It turns out that variance estimation for the doubly robust estimator is a challenging problem. While double robustness is a desirable property for point estimation, it creates a dilemma for variance estimation. The estimator ${\widehat{\mu }}_{\text{DR}2}$ given in (4.11) is consistent if either the propensity score model $q$ or the outcome regression model $\xi$ is correctly specified. There is no need to know which model is correctly specified, which is the most crucial part behind double robustness. This ambiguous feature, however, becomes a problem for variance estimation. The asymptotic variance formula under the model $q$ is usually different from the one under the model $\xi ,$ and consequently, it is difficult to construct a consistent variance estimator with unknown scenarios on model specifications.

There have been several strategies proposed in the literature on variance estimation for the doubly robust estimators. A naive approach is to use the variance estimator derived under the assumed propensity score model $q$ and take the risk that such a variance estimator might have non-negligible biases under the outcome regression model. One good news is that, under the propensity score model, the estimation of the parameters $\beta$ for the outcome regression model has no impact asymptotically on the variance of doubly robust estimators. This can be seen by using ${\widehat{\mu }}_{\text{DR}1}$ of (4.10) as an example. Let ${\widehat{m}}_i=m\left(x_i\mathrm{,}\widehat{\beta }\right),$ where $\widehat{\beta }$ is obtained based on the working model (3.1) which is not necessarily correct. Let ${\beta }^{*}$ be the probability limit of $\widehat{\beta }$ such that $\widehat{\beta }={\beta }^{*}+O_p\left(n_A^{-1/2}\right)$ regardless of the true outcome regression model (White, 1982). Let $m_i^{*}=m\left(x_i\mathrm{,}{\beta }^{*}\right)$ and $a(x,\beta )=\partial m(x,\beta )/\partial \beta .$ It can be seen that

$$\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^B{\widehat{m}}_i}-\frac{1}{N}{\displaystyle \sum _{i\in S_A}\frac{{\widehat{m}}_i}{{\widehat{\pi }}_i^A}}=\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^B{m}_i^{*}}-\frac{1}{N}{\displaystyle \sum _{i\in S_A}\frac{m_i^{*}}{{\widehat{\pi }}_i^A}}+{\left\{B\left({\beta }^{*}\right)\right\}}^{\prime }\left(\widehat{\beta }-{\beta }^{*}\right)+o_p\left(n_A^{-1/2}\right),$$

where

$$B\left({\beta }^{*}\right)=\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^{B}a\left(x_i\mathrm{,}{\beta }^{*}\right)}-\frac{1}{N}{\displaystyle \sum _{i\in S_A}\frac{a\left(x_i\mathrm{,}{\beta }^{*}\right)}{{\widehat{\pi }}_i^A}}.(6.2)$$

Since the two terms on the right hand side of (6.2) are both consistent estimators of $N^{-1}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^{N}a\left(x_i\mathrm{,}{\beta }^{*}\right)},$ we conclude that $B\left({\beta }^{*}\right)\mathrm{<mo>=</mo>}{o}_p(1)$ and

$$\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^B{\widehat{m}}_i}-\frac{1}{N}{\displaystyle \sum _{i\in S_A}\frac{{\widehat{m}}_i}{{\widehat{\pi }}_i^A}}\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^B{m}_i^{*}}-\frac{1}{N}{\displaystyle \sum _{i\in S_A}\frac{m_i^{*}}{{\widehat{\pi }}_i^A}}+o_p\left(n_A^{-1/2}\right)\mathrm{.}$$

It follows that

$${\widehat{\mu }}_{\text{DR}1}=\frac{1}{N}{\displaystyle \sum _{i\in S_A}\frac{y_i-m_i^{*}}{{\widehat{\pi }}_i^A}}+\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^B{m}_i^{*}}+o_p\left(n_A^{-1/2}\right)\mathrm{.}$$

The same arguments apply to ${\widehat{\mu }}_{\text{DR}2}.$ We can treat $\widehat{\beta }$ as if it is fixed in deriving the asymptotic variance for ${\widehat{\mu }}_{\text{DR}1}$ and ${\widehat{\mu }}_{\text{DR}2}$ under the assumed propensity score model. The techniques described in Section 6.2 can be directly used where the first estimating function in (6.1) is replaced by the one for defining ${\widehat{\mu }}_{\text{DR}1}$ or ${\widehat{\mu }}_{\text{DR}2}.$ See Theorem 2 of Chen et al. (2020) for further details. The variance estimator derived under the propensity score model, however, is generally biased under the outcome regression model.

Chen et al. (2020) also described a technique using the original idea presented in Kim and Haziza (2014) for the construction of the so-called doubly robust variance estimator. The technique is a delicate one with some theoretical attractiveness but has various issues for practical applications. We use ${\widehat{\mu }}_{\text{DR}1}$ as an example to illustrate the steps for the construction of the doubly robust variance estimator. Let

$$\widehat{\mu }\left(\alpha ,\beta \right)=\frac{1}{N}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^N{R}_i\frac{y_i-m\left(x_i\mathrm{,}\beta \right)}{\pi \left(x_i\mathrm{,}\alpha \right)}}+\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^{B}m\left(x_i\mathrm{,}\beta \right)}.$$

It follows that ${\widehat{\mu }}_{\text{DR}1}\mathrm{<mo>=</mo>}\widehat{\mu }(\widehat{\alpha },\widehat{\beta })$ if $\widehat{\alpha }$ and $\widehat{\beta }$ are from the original estimation methods. The first step is to modify the estimation of $\alpha$ and $\beta$ such that $\widehat{\alpha }$ and $\widehat{\beta }$ are obtained as solutions to

$$\frac{\partial \widehat{\mu }\left(\alpha ,\beta \right)}{\partial \alpha }\mathrm{<mo>=</mo>}0\text{and}\frac{\partial \widehat{\mu }\left(\alpha ,\beta \right)}{\partial \beta }=0.(6.3)$$

Under the logistic regression model $q$ where $\text{logit}\left\{\pi (x_i\mathrm{,}\alpha )\right\}=x_i^{\prime }\alpha$ and the linear regression model $\xi$ where $m\left(x_i\mathrm{,}\beta \right)=x_i^{\prime }\beta ,$ the equation system (6.3) becomes

$$\frac{1}{N}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^N{R}_i\left\{\frac{1}{\pi \left(x_i\mathrm{,}\alpha \right)}-1\right\}\left(y_i-x_i^{\prime }\beta \right)x_i}=0,(6.4)$$

$$\frac{1}{N}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^N\frac{R_i{x}_i}{\pi \left(x_i\mathrm{,}\alpha \right)}}-\frac{1}{N}{\displaystyle \sum _{i\in S_B}{d}_i^B{x}_i}=0.(6.5)$$

The estimating equations in (6.5) are unbiased under the joint $qp$ -framework. They are identical to (4.5) discussed in Section 4.1.2. The estimating equations in (6.4) are also unbiased under the outcome regression model, but they are different from the quasi score equations given in (3.2). The estimators $\widehat{\alpha }$ and $\widehat{\beta }$ obtained as solutions to (6.4) and (6.5) are less stable than those from standard methods. In addition, the equations system (6.4) and (6.5) will not have a solution if $\alpha$ and $\beta$ are not of the same dimension, since the number of equations in (6.4) is decided by the dimension of $\alpha$ and the number of equations in (6.5) is the same as the dimension of $\beta .$ The final estimator ${\widehat{\mu }}_{\text{DR}}=\widehat{\mu }\left(\widehat{\alpha },\widehat{\beta }\right)$ also suffers from efficiency losses when $\alpha$ and $\beta$ are estimated by solving (6.4) and (6.5).

The reason behind the use of the equations system (6.3) is purely technical. It can be shown through a first order Taylor series expansion that the estimators $\widehat{\alpha }$ and $\widehat{\beta }$ obtained from (6.3) have no impact asymptotically on the variance of ${\widehat{\mu }}_{\text{DR}}=\widehat{\mu }\left(\widehat{\alpha },\widehat{\beta }\right).$ This technical maneuver enables that simple explicit expressions for the variance $V_{qp}\left({\widehat{\mu }}_{\text{DR}}\right)$ under the $qp$ framework and for the prediction variance $V_{\xi p}\left({\widehat{\mu }}_{\text{DR}}-{\mu }_y\right)$ under the $\xi p$ framework can easily be obtained. Construction of the doubly robust variance estimator for ${\widehat{\mu }}_{\text{DR}}$ starts with the plug-in estimator for $V_{qp}\left({\widehat{\mu }}_{\text{DR}}\right)$ under the propensity scores model $q.$ A bias-correction term is then added to obtain a valid estimator for $V_{\xi p}\left({\widehat{\mu }}_{\text{DR}}-{\mu }_y\right)$ under the outcome regression model $\xi .$ The happy ending of the story is that the bias-correction term has the analytic form $N^{-2}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^N\left(R_i/{\pi }_i^A-1\right){\sigma }_i^2}$ where ${\sigma }_i^2=E_{\xi }\left(y_i|x_i\right),$ which is negligible under the propensity score model $q.$ The bias-corrected variance estimator is valid under either the propensity score model or the outcome regression model.

A doubly robust variance estimator for the commonly used ${\widehat{\mu }}_{\text{DR}2}$ is not available in the literature. A practical solution is to use bootstrap methods. Chen et al. (2022) demonstrated that standard with-replacement bootstrap procedures applied separately to $S_A$ and $S_B$ provide doubly robust confidence intervals using the pseudo empirical likelihood approach to non-probability survey samples when the reference sample is selected by single stage unequal probability sampling designs. Complications will arise when the probability sample $S_B$ uses stratified multi-stage sampling methods, a known challenge for variance estimation with complex surveys. Construction of doubly robust variance estimators for the doubly robust estimator ${\widehat{\mu }}_{\text{DR}2}$ under general settings deserves efforts in future research.

• Previous
• Article main page
• Next