2 A general procedure for constructing fully efficient replication weights

Jae Kwang Kim and Changbao Wu

Previous | Next

In principle, we can construct replication weights for any measurable sampling design, using the method outlined in Fay (1984) and Fay and Dippo (1989), such that the resulting replication variance estimators are algebraically equivalent to the standard linearization variance estimators.

Let $U=\{1,2,\dots ,N\}$ be the set of $N$ units in the finite population and $S=\{1,2,\dots ,n\}$ be the set of $n$ units in the sample, selected according to a probability sampling design. Let $w_i=1/{\pi }_i$ be the basic design weight, where ${\pi }_i=P(i\in S)$ is the first order inclusion probability for unit $i.$

Let $y_i$ be the value of the study variable $y$ for unit $i$ and $t_y={\displaystyle {\sum }_{i=1}^N{y}_i}$ be the population total of interest. The Horvitz-Thompson estimator of $t_y$ is given by

${\widehat{t}}_y={\displaystyle \sum _{i\in S}{w}_i{y}_i}.$(2.1)

The estimator ${\widehat{t}}_y$ given in (2.1) is also called the expansion estimator, with the basic design weight $w_i$ denoting the number of units in the population represented by unit $i$ in the sample. The standard variance estimator of ${\widehat{t}}_y$ can be written as

$v={\displaystyle \sum _{i\in S}{\displaystyle \sum _{j\in S}{\Omega }_{ij}{y}_i{y}_j}},$(2.2)

where ${\Omega }_{ij}=({\pi }_{ij}-{\pi }_i{\pi }_j)/({\pi }_{ij}{\pi }_i{\pi }_j)$ and ${\pi }_{ij}=P(i,j\in S)$ is the second order joint inclusion probability for $(ij).$ It is assumed that ${\pi }_{ij}>0$ for all $(ij).$ Note that ${\pi }_{ii}={\pi }_i.$ The standard variance estimator $v$ is often viewed as fully efficient since it is the Horvitz-Thompson estimator of the design-based variance $V({\widehat{t}}_y).$

Let $\Delta =({\Omega }_{ij})$ be an $n\times n$ matrix. We can re-write (2.2) as $v=y^{\prime }\Delta y,$ where $y=(y_1,y_2,\dots ,y_n{)}^{\prime }$ is the vector of sampled $y_i’s.$ The matrix $\Delta$ is nonnegative definite and can be decomposed as

$\Delta ={\displaystyle \sum _{k=1}^p{\lambda }_k{\delta }_k{{\delta }^{\prime }}_k}$(2.3)

for some ${\lambda }_k>0$ and some $n\text{–}$ dimensional vectors ${\delta }_k,k=1,2,\dots ,p.$ The most well-known decomposition (2.3) is given by the spectral decomposition where ${\delta }_k$ is the eigenvector associated with the eigenvalue ${\lambda }_k.$ In practice, very small eigenvalues are often ignored for computational reasons. For stratified sampling, the matrix $\Delta$ is block-diagonal so the computational burden may be alleviated. However, we do not restrict (2.3) to the spectral decomposition. Any decomposition satisfying (2.3) can be used.

Suppose that we want to express the fully efficient variance estimator $v$ given by (2.2) as a replication variance estimator in the form of

$v_R={\displaystyle \sum _{k=1}^L{c}_k{({\widehat{t}}_y^{(k)}-{\widehat{t}}_y)}^2,}$(2.4)

where ${\widehat{t}}_y^{(k)}={\displaystyle {\sum }_{i\in S}{w}_i^{(k)}{y}_i,w^{(k)}=(w_1^{(k)},\dots ,w_n^{(k)}{)}^{\prime }}$ is the $k^{\text{th}}$ set of replication weights, $c_k>0$ is the factor associated with the $k^{\text{th}}$ set of replication weights and $L$ is the total number of replications; see Kim, Navarro and Fuller (2006) for further discussion.

The form given by (2.4) does not include all replication variance estimators. For instance, Campbell (1980) provided a jackknife variance estimator where the pseudovalues are derived based on the von Mises approximation to the parameter of interest. Nevertheless, most replication variance estimators can be put in this form.

We have the following result on the construction of $w^{(k)}$ for $v_R$ based on the decomposition (2.3).

Theorem 1. The fully efficient variance estimator $v$ and the replication variance estimator $v_R$ are algebraically identical if we let $L=p$ and $w^{(k)}=w+{({\lambda }_k/c_k)}^{1/2}{\delta }_k,$ where $w=(w_1,\dots ,w_n{)}^{\prime }$ is the set of original basic design weights.

Proof. The proof follows directly from the fact that $v=y^{\prime }\Delta y={\displaystyle {\sum }_{k=1}^p{\lambda }_k}{({{\delta }^{\prime }}_{k}y)}^2$ and that ${\widehat{t}}_y^{(k)}-{\widehat{t}}_y=(w^{(k)}-w{)}^{\prime }y={({\lambda }_k/c_k)}^{1/2}{{\delta }^{\prime }}_{k}y.$

The choices of $c_k’\text{s}$ can be arbitrary and bear no impact on the validity and efficiency of the replication variance estimators. However, certain choices of $c_k$ will result in replication weights with negative values, which is undesirable as it may produce negative replicates for the parameters that are always positive. In practical situations one can always choose relatively large $c_k$ to avoid negative values for replication weights. In our simulation study (Case I) reported in Section 5, the problem of negative replication weights can be eliminated with the choice of $c_k=1.$

The replication variance estimator $v_R={\displaystyle {\sum }_{k=1}^L{c}_k}{({\widehat{t}}_y^{(k)}-{\widehat{t}}_y)}^2$ with $L=p$ and replication weights $w^{(k)}=w+{({\lambda }_k/c_k)}^{1/2}{\delta }_k$ is fully efficient for an arbitrary variable $y.$ It also provides fully efficient variance estimator for $\widehat{\theta }$ when $\theta$ is a smooth function of population means or totals. Practical implementation of the method depends crucially on two related issues: (i) the feasibility of the decomposition of the $n\times n$ matrix $\Delta$ specified in (2.3); and (ii) the number of sets of replication weights required to achieve the full efficiency determined by $p=\text{rank}(\Delta ).$

As for the first issue, modern advances in computational power and improved performances of available software packages make it possible to do the spectral decomposition with relatively large $n.$ For instance, on a 12-CPU unix machine with 96 gigabytes of memory, the R function eigen() can handle matrices of sizes at least as large as $n=$ 4,000. Note that the computational task involved here is for survey runners at the data preparation stage and is not for users of the data files. As for the second issue, the value of $p$ is related to the given sampling design. For simple random sampling without replacement, we have

$\Delta =N^2(1-n/N){(n(n-1))}^{-1}(I_n-1_n{1^{\prime }}_n/n),$

where $I_n$ is the $n\times n$ identity matrix and $1_n=(1,1,\dots ,1{)}^{\prime }$ is the $n\times 1$ vector of $1’\text{s}.$ It follows that $p=\text{rank}(\Delta )=\text{trace}(I_n-1_n{1^{\prime }}_n/n)=n-1.$  This is typically the case for single stage unequal probability sampling designs. For stratified simple random sampling, we have $p=n-H,$ where $H$ is the total number of strata.

It should be noted that $p\le n$ for any sampling design and the exact value of $p$ is not required for the proposed procedure to be implemented. However, since the values of $p$ and $n$ have the same order of magnitude, the proposed method requires a large number of replicates whenever $n$ is large. Under the current practices in sample surveys, the fully efficient replication weights described above become immediately implementable if $p\le$ 500 and the second order inclusion probabilities ${\pi }_{ij}$ are available. When $p$ is large, a two-stage procedure to be described in Section 3 can be used to produce a small number $L_0$ sets of replication weights for public-use data files.

In some cases, the spectral decomposition (2.3) can be avoided. For example, Deville (1999) argued that the variance estimator of ${\widehat{t}}_y$ under unequal probability sampling designs with fixed sample size can be approximated by

$$v\doteq c{\displaystyle \sum _{i\in S}(1-{\pi }_i)}{\left(\frac{y_i}{{\pi }_i}-{\tilde{t}}_y\right)}^2\text{ (2.5)}$$

where $c={\left(1-{\displaystyle {\sum }_{i\in S}{b}_i^2}\right)}^{-1},b_i=(1-{\pi }_i)/{\displaystyle {\sum }_{k\in S}(1-{\pi }_k)}$ and ${\tilde{t}}_y={\displaystyle {\sum }_{i\in S}{b}_i}(y_i/{\pi }_i).$ More generally, we consider the following form of matrix $\Delta$ in $v=y^{\prime }\Delta y,$ where

$\Delta ={\Delta }_0-{\Delta }_{0}X{(X^{\prime }{\Delta }_{0}X)}^{-1}{X}^{\prime }{\Delta }_0$(2.6)

where ${\Delta }_0=\text{diag}\{{\lambda }_1,\dots ,{\lambda }_n\},{\lambda }_i>0$  for all $i=1,2,\dots ,n,X^{\prime }=(x_1,\dots ,x_n)$ and $x_i$ is a vector of design and auxiliary variables. Many elementary single-stage sampling designs take the form (2.6) for variance estimation. In particular, Deville's formula in (2.5) can be expressed as $v\doteq y^{\prime }\Delta y$ with $\Delta$ given by (2.6), where ${\lambda }_i=c{\pi }_i^{-2}(1-{\pi }_i)$ in ${\Delta }_0$ and $x_i={\pi }_i.$ The conditional Poisson sampling design to be discussed in Section 5 also takes the form (2.6) where $x_i$ are the design variables in the design constraint ${\displaystyle {\sum }_{i\in S}{\pi }_i^{-1}}{x}_i={\displaystyle {\sum }_{i=1}^N{x}_i}.$

For the matrix given by (2.6), it can be shown that

$y^{\prime }\Delta y=(y-X\widehat{\beta }{)}^{\prime }{\Delta }_0(y-X\widehat{\beta }),$

where $\widehat{\beta }={(X^{\prime }{\Delta }_{0}X)}^{-1}{X}^{\prime }{\Delta }_{0}y.$ Thus, we have

$y^{\prime }\Delta y={\displaystyle \sum _{k=1}^n{\lambda }_k}{(y_k-{x^{\prime }}_k\widehat{\beta })}^2,$(2.7)

which is useful in deriving an expression for replication variance estimator in the form given by (2.4). The fully efficient variance estimator $v$ in (2.7) and the replication variance estimator $v_R$ in (2.4) are algebraically identical if we let $L=n$ and $w^{(k)}=w+{({\lambda }_k/c_k)}^{1/2}{\delta }_k,$ where $w=(w_1,\dots ,w_n{)}^{\prime }$ is the set of original basic design weights and ${\delta }_k=({\delta }_{1k},\dots ,{\delta }_{nk}{)}^{\prime }$ with

               ${\delta }_{ik}=\{\begin{array}{ll}-1+{x^{\prime }}_k{(X^{\prime }{\Delta }_{0}X)}^{-1}{x}_i{\lambda }_i\hfill & \text{if }i=k\hfill \\ {x^{\prime }}_k{(X^{\prime }{\Delta }_{0}X)}^{-1}{x}_i{\lambda }_i\hfill & \text{otherwise}\text{.}\hfill \end{array}$

The proof follows directly from the fact that ${{\delta }^{\prime }}_{k}y=-y_k+{x^{\prime }}_k\widehat{\beta },y^{\prime }\Delta y={\displaystyle {\sum }_{k=1}^n{\lambda }_k}{({{\delta }^{\prime }}_{k}y)}^2$ and that ${\widehat{t}}_y^{(k)}-{\widehat{t}}_y=(w^{(k)}-w{)}^{\prime }y={({\lambda }_k/c_k)}^{1/2}{{\delta }^{\prime }}_{k}y.$

Previous | Next