3 Sparse and efficient replication weights

Jae Kwang Kim and Changbao Wu

Previous | Next

Large-scale complex surveys usually have a relatively large sample size ranging from a few hundreds to many thousands. The fully efficient replication weights described in Section 2 or replication weights constructed by some existing methods such as the jackknife or the bootstrap methods would involve a very large number of sets of weights. Although valid replication weights provide enormous convenience to the users of survey data, who are not necessarily the survey runners, the burden of manipulating a data set with hundreds or even thousands of replicate weights can be enormous. As a result, how to achieve efficient replication variance estimation with a relatively small number of replicate weights is a question with both theoretical and practical value.

We propose two strategies to construct sparse and efficient replication weights. We start with a large number $L$ sets of replication weights. These initial weights may be produced using the general method described in Section 2 or by existing methods. Suppose they can be viewed as fully efficient. The first strategy is to select a small number $L_0$ sets of weights from the initial large number $L$ sets of weights using a probability sampling method. The small number $L_0$ satisfies the desired sparsity and the random selection procedure guarantees validity of the resulting variance estimators. The second strategy is to achieve efficiency through a novel weight-calibration procedure. The $L_0$ sets of calibrated replication weights provide fully efficient variance estimators for variables used in the calibration and also highly efficient variance estimators for variables related to calibration variables.

3.1  Achieve sparsity and efficiency through random sampling

Suppose that the fully efficient replication variance estimator is given by $v_R={\displaystyle {\sum }_{k=1}^L{c}_k{({\widehat{t}}_y^{(k)}-{\widehat{t}}_y)}^2,}$ with replication weights constructed by using Theorem 1. Observe that $v_R$ can be viewed as a finite population total. If we want to use $L_0(<L)$ sets of replication weights to provide valid inference for variance estimation, the following simple strategy can be used. First, select $L_0$ sets of weights from the original $L$ sets of weights by simple random sampling without replacement. For notational simplicity and without loss of generality, we denote the selected sets of weights by $w^{(j)},j=1,2,\dots ,L_0.$ Then, calculate the replication variance estimator of ${\widehat{t}}_y$ based on the $L_0$ sets of weights as

$$v_R^{(1)}=\frac{L}{L_0}{\displaystyle \sum _{j=1}^{L_0}{c}_j}{({\widehat{t}}_y^{(j)}-{\widehat{t}}_y)}^2.\text{ (3.1)}$$

The variance estimator $v_R^{(1)}$ is still unbiased for an arbitrary variable $y,$ since $E^{*}(v_R^{(1)})={\displaystyle {\sum }_{k=1}^L{c}_k}{({\widehat{t}}_y^{(k)}-{\widehat{t}}_y)}^2=v_R,$ where $E^{*}(\cdot )$ denotes the expectation under the random selection of $L_0$ sets of weights.

An alternative form of the replication variance estimator based on the $L_0$ sets of weights can be derived as follows. Noting that ${\widehat{t}}_y^{(k)}-{\widehat{t}}_y={({\lambda }_k/c_k)}^{1/2}{{\delta }^{\prime }}_{k}y,$ we can re-write the fully efficient variance estimator as

$v_R=m\left\{{\displaystyle \sum _{k=1}^L{\lambda }_k{({{\delta }^{\prime }}_{k}y)}^2}/{\displaystyle \sum _{k=1}^L{\lambda }_k}\right\},$

where $m={\displaystyle {\sum }_{k=1}^L{\lambda }_k.}$ The ${\delta }_k’\text{s}$ are orthogonal eigenvectors satisfying $\left|{\delta }_k\right|=1$ under spectral decomposition and ${{\delta }^{\prime }}_{k}y$ are projections of $y$ onto the $n\text{–}$ dimensional unit-sphere. It is very natural to use the following weighted version for the variance estimator of ${\widehat{t}}_y$ based on the $L_0$ randomly selected sets of weights:

$$v_R^{(2)}=m\left\{{\displaystyle \sum _{j=1}^{L_0}{\lambda }_j{({{\delta }^{\prime }}_{j}y)}^2}/{\displaystyle \sum _{j=1}^{L_0}{\lambda }_j}\right\}=\frac{m}{m_0}{\displaystyle \sum _{j=1}^{L_0}{c}_j{({\widehat{t}}_y^{(j)}-{\widehat{t}}_y)}^2,}\text{ (3.2)}$$

where $m_0={\displaystyle {\sum }_{j=1}^{L_0}{\lambda }_j.}$ Noting that $v_R^{(2)}$ is a ratio estimator of $v_R,$ it is usually more efficient than $v_R^{(1)}\text{}.$

A third version of the replication variance estimator can be constructed by first selecting $L_0$ sets of weights with unequal probabilities and then using a Horvitz-Thompson estimator of $v_R.$ Note that we can view $v_R={\displaystyle {\sum }_{k=1}^L{\lambda }_k}{({{\delta }^{\prime }}_{k}y)}^2$ as a population total and ${\lambda }_k$ as a size variable. We select $L_0$ sets of weights from the original $L$ sets of weights with inclusion probabilities ${\eta }_k$ proportional to ${\lambda }_k.$ The resulting variance estimator of ${\widehat{t}}_y$ is given by

$v_R^{(3)}={\displaystyle \sum _{j=1}^{L_0}\frac{c_j}{{\eta }_j}}{({\widehat{t}}_y^{(j)}-{\widehat{t}}_y)}^2,$(3.3)

where ${\eta }_j=L_0{\lambda }_j/{\displaystyle {\sum }_{k=1}^L{\lambda }_k}.$ It turns out that the eigenvalues ${\lambda }_k$ differ substantially in magnitude and using ${\lambda }_k$ as size measure leads to a large portion of the $L$ sets of weights being selected with probability one. In the simulation study described in Section 5, we also included a fourth version of the replication variance estimator, denoted as $v_R^{(4)},$ with ${\lambda }_k^{1/2}$ as the size measure and ${\eta }_j=L_0{\lambda }_j^{1/2}/{\displaystyle {\sum }_{k=1}^L{\lambda }_k^{1/2}}.$

Another possible version of the replication variance estimator is to simply select $L_0$ sets of weights corresponding to the $L_0$ largest values of ${\lambda }_k$ and then use $v_R^{(2)}.$ Simulation results, not reported here, showed that the resulting variance estimator is severely biased and shouldn't be used in practice.

3.2  Achieve sparsity and efficiency through weight-calibration

We now discuss a novel approach of achieving sparsity without losing the efficiency of the variance estimators for some key variables. Suppose $L_0$ is the desired replication size, which is much smaller than the sample size $n.$ For example, the Natural Resources Inventory Survey (sponsored by the US department of Agriculture) used $L_0=$ 29 while the PSU sample size can be as large as $n=$ 300,000. We present a weight-calibration technique that not only allows the use of a small $L_0,$ but also provides fully efficient variance estimators for key population parameters. Our proposed strategy for constructing the smaller number $L_0$ sets of replication weights consists of the following four steps:

Step 1. Choose a set of key variables for which full efficiency of the variance estimator is desired. Let $z_i=(z_{i1},\dots ,z_{im}{)}^{\prime }$ be the vector of key variables for the $i^{\text{th}}$ unit included in the survey data file, where $m\le L_0.$ Among them can be important auxiliary variables and study variables as well as design variables. Let ${\widehat{t}}_z={\displaystyle {\sum }_{i\in S}{w}_i{z}_i.}$ Let $v_1({\widehat{t}}_z)$ be an $m\times m$ estimated variance-covariance matrix for ${\widehat{t}}_z$ computed by the standard linearization method or by a replication method that is fully efficient.

Step 2. Construct an initial $L_0$ sets of replication weights that produce an asymptotically unbiased variance estimator. These initial replicates can be obtained by a bootstrap method with $L_0$ replicates, or by the delete-a-group jackknife method of Kott (2001), or by the sampling method described in Section 3.1. Let $w_0^{(k)}=(w_{10}^{(k)},\dots ,w_{n0}^{(k)}{)}^{\prime },k=1,\dots ,L_0$  be the initial sets of weights. Denote ${\widehat{t}}_{y0}^{(k)}={\displaystyle {\sum }_{i\in S}{w}_{i0}^{(k)}{y}_i}$ and let

$v_0={\displaystyle \sum _{k=1}^{L_0}{c}_{k0}{({\widehat{t}}_{y0}^{(k)}-{\widehat{t}}_y)}^2}$(3.4)

be the replication variance estimator based on the $L_0$ sets of weights.

We can apply the variance formula (3.4) to the vector of key variables $z$ to get $v_0({\widehat{t}}_z)={\displaystyle {\sum }_{k=1}^{L_0}{c}_{k0}({\widehat{t}}_{z0}^{(k)}-{\widehat{t}}_z)({\widehat{t}}_{z0}^{(k)}-{\widehat{t}}_z}{)}^{\prime },$ where ${\widehat{t}}_{z0}^{(k)}={\displaystyle {\sum }_{i\in S}{w}_{i0}^{(k)}{z}_i}.$ Note that $v_0({\widehat{t}}_z)$ is not as efficient as $v_1({\widehat{t}}_z)$ obtained in Step 1.

Step 3. Decompose the nonnegative definite variance-covariance matrix $v_1({\widehat{t}}_z)$ as

$v_1({\widehat{t}}_z)={\displaystyle \sum _{k=1}^m{\alpha }_k{q}_k{q^{\prime }}_k}$(3.5)

using the spectral decomposition or any other suitable methods. Let ${\alpha }_k=0$ for $k=m+1,\dots ,L_0$ and define

${\widehat{t}}_z^{(k)}={\widehat{t}}_z+{({\alpha }_k/c_{k0})}^{1/2}{q}_k,k=1,2,\dots ,L_0.$

It follows that the $L_0$ pseudo-replicates ${\widehat{t}}_z^{(k)}$ defined above satisfy

${\displaystyle \sum _{k=1}^{L_0}{c}_{k0}}({\widehat{t}}_z^{(k)}-{\widehat{t}}_z)({\widehat{t}}_z^{(k)}-{\widehat{t}}_z{)}^{\prime }=v_1({\widehat{t}}_z),$(3.6)

due to the decomposition to $v_1({\widehat{t}}_z)$ given in (3.5). It should be noted that (3.5) bears no relation to the decomposition to $\Delta$ described in Section 2 and the condition $m\le L_0$ is required to make (3.6) possible.

Step 4. Improve the efficiency of $v_0$ computed from (3.4) for an arbitrary $y$ variable through a weight-calibration procedure. For the $k^{\text{th}}$ set of initial weights $w_0^{(k)}=(w_{10}^{(k)},\dots ,w_{n0}^{(k)}{)}^{\prime },$ the calibrated weights $w_c^{(k)}=(w_{1c}^{(k)},\dots ,w_{nc}^{(k)}{)}^{\prime }$ minimize the chi-square distance measure

$\Phi (w_c^{(k)},w_0^{(k)})={\displaystyle \sum _{i\in S_k}{\tau }_i}{(w_{ic}^{(k)}-w_{i0}^{(k)})}^2/w_{i0}^{(k)}$(3.7)

subject to the constraint

${\displaystyle \sum _{i\in S}{w}_{ic}^{(k)}{z}_i={\widehat{t}}_z^{(k)},}$(3.8)

where $S_k=\left\{i|i\in S;w_{i0}^{(k)}>0\right\},$ the ${\tau }_i’\text{s}$ are known constants, and ${\widehat{t}}_z^{(k)}$ is the $k^{\text{th}}$ pseudo replicate of ${\widehat{t}}_z$ defined in Step 3.

The calibrated weights $w_c^{(k)}=(w_{1c}^{(k)},\dots ,w_{nc}^{(k)}{)}^{\prime },k=1,2,\dots ,L_0$  are used in (3.4) to compute the final replication variance estimator $v_C({\widehat{t}}_y)={\displaystyle {\sum }_{k=1}^{L_0}{c}_{k0}}{({\widehat{t}}_{yc}^{(k)}-{\widehat{t}}_y)}^2,$ where ${\widehat{t}}_{yc}^{(k)}={\displaystyle {\sum }_{i\in S}{w}_{ic}^{(k)}}{y}_i.$

The proposed weight-calibration procedure ensures that the replication estimator $v_C({\widehat{t}}_z)$ based on the $L_0$ sets of calibrated weights matches exactly the fully efficient estimator $v_1({\widehat{t}}_z),$ due to the calibration constraints (3.8) and the equation (3.6). Furthermore, the calibrated replication weights provide more efficient variance estimators for an arbitrary $y$ that is related to $z.$ To see this, we re-write ${\widehat{t}}_{y0}^{(k)}={\displaystyle {\sum }_{i\in S}{w}_{i0}^{(k)}}{y}_i$ as

${\widehat{t}}_{y0}^{(k)}={\widehat{t}}_{e0}^{(k)}+({\widehat{t}}_{z0}^{(k)}{)}^{\prime }\widehat{\beta },$(3.9)

where ${\widehat{t}}_{e0}^{(k)}={\displaystyle {\sum }_{i\in S}{w}_{i0}^{(k)}{\widehat{e}}_i,{\widehat{e}}_i}=y_i-{z^{\prime }}_i\widehat{\beta }$ and $\widehat{\beta }={\left\{{\displaystyle {\sum }_{i\in S}\text{}}{w}_i{z}_i{z^{\prime }}_i/{\tau }_i\right\}}^{-1}{\displaystyle {\sum }_{i\in S}{w}_i{z}_i{y}_i}/{\tau }_i.$ Let ${\widehat{t}}_e={\displaystyle {\sum }_{i\in S}{w}_i{\widehat{e}}_i.}$ The variance estimator of ${\widehat{t}}_y$ based on the initial $L_0$ sets of weights can be expressed as

$v_0({\widehat{t}}_y)={\displaystyle \sum _{k=1}^{L_0}{c}_{k0}}{({\widehat{t}}_{y0}^{(k)}-{\widehat{t}}_y)}^2$
$={\displaystyle \sum _{k=1}^{L_0}{c}_{k0}}{({\widehat{t}}_{e0}^{(k)}-{\widehat{t}}_e)}^2+{\displaystyle \sum _{k=1}^{L_0}{c}_{k0}}{(({\widehat{t}}_{z0}^{(k)}{)}^{\prime }\widehat{\beta }-({\widehat{t}}_z{)}^{\prime }\widehat{\beta })}^2+2{\displaystyle \sum _{k=1}^{L_0}{c}_{k0}}({\widehat{t}}_{e0}^{(k)}-{\widehat{t}}_e)(({\widehat{t}}_{z0}^{(k)}{)}^{\prime }\widehat{\beta }-({\widehat{t}}_z{)}^{\prime }\widehat{\beta })$
$=v_0({\widehat{t}}_e)+{\widehat{\beta }}^{\prime }{v}_0({\widehat{t}}_z)\widehat{\beta }+2{\widehat{\beta }}^{\prime }{\mathrm{cov}}_0({\widehat{t}}_e,{\widehat{t}}_z),$

where ${\mathrm{cov}}_0({\widehat{t}}_e,{\widehat{t}}_z)$ is the estimated covariance between ${\widehat{t}}_e^{}$ and ${\widehat{t}}_z^{}$ based on the initial $L_0$ sets of replication weights. In many designs, we can choose a suitable ${\tau }_i$ such that $\text{Cov}({\widehat{t}}_e,{\widehat{t}}_z)\doteq 0.$ This is the case, for instance, with the choice of ${\tau }_i={(w_i-1)}^{-1}$ or ${\tau }_i=w_i^{-1}$ under Poisson sampling. Fuller (1998) discussed the required conditions in the context of two-phase sampling. It follows that

$v_0({\widehat{t}}_y)\doteq v_0({\widehat{t}}_e)+\widehat{{\beta }^{\prime }}{v}_0({\widehat{t}}_z)\widehat{\beta }.$(3.10)

Using similar argument, it can be shown that the variance estimator based on the $L_0$ sets of calibrated weights satisfies

$v_C({\widehat{t}}_y)\doteq v_0({\widehat{t}}_e)+\widehat{{\beta }^{\prime }}{v}_1({\widehat{t}}_z)\widehat{\beta }.$(3.11)

The variance estimator $v_C({\widehat{t}}_y)$ given by (3.11) is generally more efficient than $v_0({\widehat{t}}_y)$ given by (3.10), due to the use of $v_1({\widehat{t}}_z)$ instead of $v_0({\widehat{t}}_z).$ The gain of efficiency depends on the relative magnitude of $\widehat{{\beta }^{\prime }}{v}_1({\widehat{t}}_z)\widehat{\beta }$ over $v_0({\widehat{t}}_e).$ If $y$ is highly correlated with $\widehat{y}=z^{\prime }\widehat{\beta },$ the variance of the residual term $v_0({\widehat{t}}_e)$ will be relatively small. In this case $v_C({\widehat{t}}_y)$ will be highly efficient. On the other hand, if there is no correlation between $y$ and $\widehat{y},$ then no improvement will be achieved by using the calibrated weights $w_c^{(k)};$ see also Theorem 3 in Section 4.

One of the drawbacks of the chi-square distance $\Phi (w_c^{(k)},w_0^{(k)})$ in Step 4 is that some of the resulting calibrated weights could take negative values. To avoid negative weights, we propose replacing the chi-square distance in (3.7) by the following minimum entropy distance

$$D(w_c^{(k)},w_0^{(k)})=-{\displaystyle \sum _{i\in S_k}{\tau }_i^{-1}}\left\{w_{i0}^{(k)}\log \left(\frac{w_{ic}^{(k)}}{w_{i0}^{(k)}}\right)-w_{ic}^{(k)}+w_{i0}^{(k)}\right\}\text{ (3.12)}$$

for two reasons. First, the calibrated weights $w_{ic}^{(k)}$ are guaranteed to be positive. Second, there exists a well-behaved computational algorithm for this constrained minimization problem. It can be shown that $w_c^{(k)}$ minimizing $D(w_c^{(k)},w_0^{(k)})$ subject to (3.8) are given by

$$w_{ic}^{(k)}=\frac{w_{i0}^{(k)}/{\tau }_i}{1+{\lambda }^{\prime }{z}_i},\text{ (3.13)}$$

where the Lagrange multiplier $\lambda$ is the solution to

$g(\lambda )={\displaystyle \sum _{i\in S}\frac{w_{i0}^{(k)}{z}_i/{\tau }_i}{1+{\lambda }^{\prime }{z}_i}}-{\widehat{t}}_z^{(k)}=0.$(3.14)

An efficient computational algorithm for finding the solution $\lambda$ to (3.14) can be found in Wu (2004) and a related R function can be obtained by a minor modification of the R function presented in Wu (2005).

Previous | Next