3. Variance estimation for the one-step calibration estimator

Phillip S. Kott and Dan Liao

Previous | Next

In this section, we let

$$t_y={\displaystyle {\sum }_R{w}_k{y}_k}={\displaystyle {\sum }_R{d}_k\alpha \left(g^T{x}_k\right)y_k}$$

be the calibration-weighted estimator for $T_y,$  where $w_k=d_k\alpha \left(g^T{x}_k\right)$  when $k\in R$  is the calibration weight, and $w_k$  is conveniently defined to be 0 when $k\notin R.$  The weight-adjustment function $\alpha (\cdot )$  is defined implicitly by equation (2.4), and $g$  is again chosen so that the calibration equation (2.5) holds for either $\theta =0$  or $1.$

We propose the following estimator for the variance $t_y:$

$$v\left(t_y\right)={\displaystyle \sum _{k,j\in S}\left(1-\frac{{\pi }_k{\pi }_j}{{\pi }_{kj}}\right)\left[d_k\left(\theta z_k^{T}b+{\alpha }_k{e}_k\right)\right]\left[d_j\left(\theta z_j^{T}b+{\alpha }_j{e}_j\right)\right]}+{\displaystyle \sum _{k\in R}{d}_k\left({\alpha }_k^2-{\alpha }_k\right)e_k^2},(3.1)$$

where ${\pi }_{kj}$  is the joint selection probability of $k$  and $j$  under the original sampling design, ${\pi }_{kk}={\pi }_k=1/d_k,{\pi }_k=\alpha \left(g^T{x}_k\right)$  when $k\in R$  and 0 otherwise,

$$b={\left[{\displaystyle {\sum }_R{d}_k{\alpha }^{\prime }\left(g^T{x}_k\right)x_k{z}_k^T}\right]}^{-1}{\displaystyle {\sum }_R{d}_k{\alpha }^{\prime }\left(g^T{x}_k\right)x_k{y}_k},(3.2)$$

and $e_k=y_k-z_k^{T}b.$  We will show that $v\left(t_y\right)$  in equation (3.1) can be nearly unbiased in some sense if either a response model (Section 3.1) or prediction model holds (Section 3.2).

The variance estimator in equation (5.2) of Kott (2006) is identical to $v\left(t_y\right)$  in equation (3.1) when $\theta =0.$  The variance estimator in Kim and Haziza (2014) is also similar. Their prediction model is more general than the linear prediction model considered here.

This variance estimator $v\left(t_y\right)$  presupposes that the original sampling design is such that each element can only be drawn once. In Section 3.1, we see that when the probabilities of response are independent (Poisson), then under mild assumptions, $v\left(t_y\right)$  is a nearly unbiased estimator of the mean squared error of $t_y$  under the quasi-sampling design whether or not the prediction model, $\text{E}\left(y_k|x_k,z_k\right)=z_k^T\beta ,$  holds.

In Section 3.2, $v\left(t_y\right)$  is shown to be a nearly unbiased estimator for the combined prediction-model and original-sampling-design variance of $t_y$  as an estimator for $T_y$  whether or not the response model in equation (2.4) holds. Thus, $v\left(t_y\right)$  can be called a “simultaneous variance estimator”.

3.1 Variance estimation under the response model

For ease of exposition we will assume that the response model in equation (2.4) with a finite $u$  holds. Sufficient conditions for $v\left(t_y\right)$  to be a nearly unbiased estimator for the mean squared error of $t_y$  (by which the bias converges to 0 as the sample size grows arbitrary large) are

$${\pi }_{kj}\ge B_0>0(3.3)$$

$${\displaystyle \sum _{j=1}^N\left|\frac{{\pi }_{kj}}{{\pi }_k{\pi }_j}-1\right|\le B_1<\infty \text{ for every }k,}(3.4)$$

$$\frac{{\displaystyle {\sum }_{j=1}^N{\psi }_j^r}}{N}\le B_2<\infty \text{ where }{\psi }_j\text{is }{y}_j\text{ or any component of }{x}_j\text{ or }{z}_j,\text{while }r=1\text{ or }2,(3.5)$$

and $N^{-1}{\sum }_R{d}_k{\alpha }^{\prime }\left(g^T{x}_k\right)z_k{x}_k^T$  is of full rank and is bounded in probability as the sample size grows arbitrarily large.

From these, ${\alpha }^{\prime }\left(\varphi \right)=\left(1-\alpha \left(\varphi \right)/u\right)\exp \left(\varphi \right)/\left[\left(1+\exp \left(\varphi \right)/u\right)\right]$  being bounded when $u$  is finite, and the Cauchy-Schwarz inequality $\left({\left(\sum a_k{b}_k\right)}^2\le \sum a_k^2\sum b_k^2\right),$  it is not hard to see not only that $g$  is a consistent estimator for $\gamma ,$  but also that $b$  in equation (3.2) (which can be rendered $b={\left[N^{-1}{\sum }_R{d}_k{\alpha }^{\prime }\left(g^T{x}_k\right)x_k{z}_k^T\right]}^{-1}{N}^{-1}{\sum }_R{d}_k{\alpha }^{\prime }\left(g^T{x}_k\right)x_k{y}_k)$ has a probability limit, call it $b^{*},$  whether or not the prediction model holds. Moreover, both $b-b^{*}$  and $g-\gamma$  are $O_p\left(1/\sqrt{n}\right).$

Observe that

$$\begin{array}{lll}\left(t_y-T_y\right)/N\hfill & =\hfill & \theta \left({\displaystyle {\sum }_S{d}_k{z}_k^T{b}^{*}}-{\displaystyle {\sum }_U{z}_k^T{b}^{*}}\right)/N\hfill \\ \hfill & +\hfill & \left[{\displaystyle {\sum }_R{d}_k\alpha \left(g^T{x}_k\right)e_k^{*}}-{\displaystyle {\sum }_R{d}_k\alpha \left({\gamma }^T{x}_k\right)e_k^{*}}\right]/N\hfill \\ \hfill & +\hfill & \left[{\displaystyle {\sum }_R{d}_k\alpha \left({\gamma }^T{x}_k\right)e_k^{*}}-{\displaystyle {\sum }_U{e}_k^{*}}\right]/N,\hfill \end{array}$$

where $e_k^{*}=y_k-z_k^T{b}^{*}.$  The insertion of the ${\alpha }^{\prime }(\cdot )$  into the “regression coefficient” $b$  allows us to ignore the contribution to quasi-design mean squared error of the second term in this sum, $Q={\sum }_R{d}_k\left[\alpha \left(g^T{x}_k\right)-\alpha \left({\gamma }^T{x}_k\right)\right]e_k^{*}/N.$  That is because ${\sum }_R{d}_k{\alpha }^{\prime }\left({\gamma }^T{x}_k\right)x_k{e}_k=0$  is true by definition, which implies ${\sum }_R{d}_k{\alpha }^{\prime }\left({\gamma }^T{x}_k\right)x_k{e}_k^{*}$  is ${\text{O}}_p\left(1/\sqrt{n}\right)$  under our assumptions. Moreover, since $\alpha \left(g^T{x}_k\right)-\alpha \left({\gamma }^T{x}_k\right)={\alpha }^{\prime }\left(c_k\right){\left(g-\gamma \right)}^T{x}_k$  is also ${\text{O}}_p\left(1/\sqrt{n}\right),$   $Q={\left(g-\gamma \right)}^T{\sum }_R{d}_k{\alpha }^{\prime }\left(c_k\right)x_k{e}_k^{*}$  is ${\text{O}}_p\left(1/n\right),$  which is asymptotically ignorable relative to the two ${\text{O}}_p\left(1/\sqrt{n}\right)$  components of $\left(t_y-T_y\right)/N.$

With the contribution of $Q$  eliminated from consideration, an idealized, but not calculable, nearly unbiased estimator for the quasi-design mean squared error of $t_y$  is

$$v_{I1}\left(t_y\right)={\displaystyle \sum _{k,j\in S}\left(1-\frac{{\pi }_k{\pi }_j}{{\pi }_{kj}}\right)}\left[d_k\left(\theta z_k^T{b}^{*}+e_k^{*}\right)\right]\left[d_j\left(\theta z_j^T{b}^{*}+e_j^{*}\right)\right]+{{\displaystyle \sum _{k\in R}\left(\frac{d_k{e}_k^{*}}{p_k}\right)}}^2\left(1-p_k\right),(3.6)$$

where the first term on the right estimates the mean squared error before nonresponse (if any) and the second the added variance due to nonresponse.

An alternative nearly unbiased idealized mean squared error estimator, closer to being calculable, is

$$v_{I2}\left(t_y\right)={\displaystyle \sum _{k,j\in S}\left(1-\frac{{\pi }_k{\pi }_j}{{\pi }_{kj}}\right)}\left[d_k\left(\theta z_k^T{b}^{*}+\frac{R_k}{p_k}{e}_k^{*}\right)\right]\left[d_j\left(\theta z_j^T{b}^{*}+\frac{R_j}{p_j}{e}_j^{*}\right)\right]+{{\displaystyle \sum _{k\in R}{d}_k\left(\frac{e_k^{*}}{p_k}\right)}}^2\left(1-p_k\right),(3.7)$$

where again $R_k=1$  when $k\in R,0$  otherwise. Since the $\left(R_k/p_k\right)e_k^{*}$  are independent under the response model with mean $e_k^{*}$  and variance ${\left(e_k^{*}/p_k\right)}^2{p}_k\left(1-p_k\right),\text{E}\left[\left(R_k/p_k\right)e_k^{*}\left(R_j/p_j\right)e_j^{*}\right]=e_k^{*}{e}_j^{*}$  when $k\ne j.$  By contrast, the following holds when $k=j:$

$$\begin{array}{ll}\left(1-{\pi }_k\right)E\left[{\left(d_k\frac{R_k}{p_k}{e}_k^{*}\right)}^2\right]\hfill & =\left(1-{\pi }_k\right)\left[{\left(d_k{e}_k^{*}\right)}^2+{\left(\frac{d_k{e}_k^{*}}{p_k}\right)}^2{p}_k\left(1-p_k\right)\right]\hfill \\ \hfill & =\left(1-{\pi }_k\right){\left(d_k{e}_k^{*}\right)}^2+{\left(\frac{d_k{e}_k^{*}}{p_k}\right)}^2{p}_k\left(1-p_k\right)-d_k{\left(\frac{e_k^{*}}{p_k}\right)}^2{p}_k\left(1-p_k\right).\hfill \end{array}$$

The first summation on the right-hand side of equation (3.7) has terms where $k\ne j$  and terms where $k=j,$  the latter of which causes the second summation in (3.7) to differ from the second summation on the right-hand side of equation (3.6). Note that the expectation under the response model of ${\sum }_R{d}_k{\left(e_k^{*}/p_k\right)}^2\left(1-p_k\right)$  in the second summation on the right-hand side of (3.7) is ${\sum }_S{d}_k{\left(e_k^{*}/p_k\right)}^2{p}_k\left(1-p_k\right).$

Finally, $v_{I2}\left(t_y\right)$  can be replaced by the asymptotically identical, but computable, $v\left(t_y\right)$  in equation (3.1) since ${\sum }_{j\in S}\left(1-{\pi }_k{\pi }_j/{\pi }_{kj}\right)$  is bounded for all $k$  under assumptions (3.3) and (3.4), allowing $e_k$  and ${\alpha }_k$  to be substituted for the unknown $e_k^{*}$  and $1/p_k,$  respectively (because $e_k^{*}-e_k$  and ${\alpha }_k-1/p_k$  are ${\text{O}}_p\left(1/\sqrt{n}\right)$  for all $k).$

3.2 Variance estimation under the prediction model

Matters are a bit simpler when we assume a prediction model holds but not necessarily the response model in equation (2.4). Suppose $\text{E}\left(y_k|x_k,z_k\right)=z_k^T\beta ,$  whether or not $k$  is sampled or responds when sampled, and the ${\epsilon }_k=y_k-z_k^T\beta$  are uncorrelated random variables with variances equal to ${\sigma }_k^2=z_k^T\eta ,$  where $\eta$  need not be specified other than having finite components.

The mean squared error of $t_y$  as an estimator for $T_y$  under that prediction model is the sum of the prediction variance of $t_y$  as an estimator for $T_y,{\sum }_R\left(w_k^2-w_k\right){\sigma }_k^2$  (see, for example, Kott 2009, page 69), and the squared bias, ${\left({\sum }_S{x}_k^T\beta -{\sum }_U{x}_k^T\beta \right)}^2,$  the latter being zero when $\theta =0.$  The combined variance of $t_y$  as an estimator for $T_y$  under the prediction model and original sample design is

$$V_C=\theta {\text{Var}}_D\left({\displaystyle {\sum }_S{x}_k^T\beta }\right)+{\text{E}}_D\left[{\displaystyle {\sum }_S\left(w_k^2-w_k\right){\sigma }_k^2}\right],$$

where the subscript $D$  denotes that the operation (variance or expectation) is with respect to the original sampling design. Recall $w_k=0$  for $k\ne R.$

To see that $v\left(t_y\right)$  in equation (3.1) provides a nearly unbiased estimator for $V_C,$  observe first that

$$e_k=y_k-z_k^{T}b={\epsilon }_k-z_k^T{\left[N^{-1}{\displaystyle {\sum }_R{d}_j{\alpha }^{\prime }\left(g^T{x}_j\right)x_j{z}_j^T}\right]}^{-1}{N}^{-1}{\displaystyle {\sum }_R{d}_j{\alpha }^{\prime }\left(g^T{x}_j\right)x_j{\epsilon }_j.}$$

Let ${\delta }_{kj}=1$  when $k=j$  and $0$ otherwise. Because the ${\epsilon }_k$  are uncorrelated, and $\text{E}\left({\epsilon }_k^2\right)={\sigma }_k=z_k^T\eta ,$  it is now not hard to show that $\text{E}\left(e_k{e}_j\right)={\delta }_{kj}{\sigma }_k^2+\text{O}\left(1/n\right)$  for almost every $k,j$  pair under the prediction model when $N^{-1}{\sum }_R{d}_k{\alpha }^{\prime }\left(g^T{x}_k\right)z_k{x}_k^T$  converges to an invertible matrix, and assumptions (3.3), (3.4), and

$$\frac{{\displaystyle \sum _{j=1}^N{\psi }_j^r}}{N}\le B_2<\infty \text{ where }{\psi }_j\text{ is any component of }{x}_j\text{ or }{z}_j,\text{ and }r\text{ = 1, 2, 3, or 4,}(3.8)$$

hold. Observe that the change from the assumptions in (3.5) to (3.8) makes the relative bias of $v\left(t_y\right)$  as an estimator for $V_C\left(\text{or }{\sum }_R\left(w_k^2-w_k\right){\sigma }_k^2\text{ when }\theta =0\right)\text{O}\left(1/n\right)$  rather than $\text{O}\left(1/\sqrt{n}\right).$

Previous | Next