14  Multipurpose Surveys

If theoretical statisticians ignore the challenge of dealing with multipurpose surveys, then the gap between them and practical statisticians will become ever wider. The design and analysis of multivariate surveys should be one of the next major research areas. Smith (1976)

This chapter presents an introduction to the research motivation provided by one of the most influential figures in the world statistical scene, Professor Smith (Smith 1976), who discussed the foundations of sampling theory from its early years to the latest trends in prediction and estimation in finite populations. In his many articles, this author stated that in sampling, univariate problems, meaning estimation of an unknown parameter for a single characteristic of interest, are found in only a few application branches, limited to public opinion surveys, industrial acceptance sampling, and audit sampling. However, the vast majority of surveys conducted around the world are multipurpose, meaning they estimate several unknown parameters for several characteristics of interest. Professor Smith noted the limitation of the great sampling classics in not considering this type of study or including them in their pages, and he called on sampling theorists to carry out formal research on these topics, as shown by the motivation at the beginning of the chapter.

14.1 Introduction

Most applications in sample surveys involve multiple study variables. This brief section presents a framework for the joint estimation of the parameters of interest under some sampling designs. With respect to the sampling design, Holmberg (2002a) and Holmberg (2002b) developed the relevant theory for selecting probability samples in multipurpose surveys, and with respect to multiparameter estimation, Gutiérrez (2009) proposes a general estimation system based on classical results from linear model theory and linear algebra.

The purpose of a sampling study is to obtain information about a particular finite population by estimating population parameters such as means, totals, proportions, or ratios. However, most surveys involve not a single characteristic, but several characteristics of interest. Classical sampling books seem to omit the fact that a survey is rarely planned to estimate only one parameter, and the theory developed by sampling researchers focuses on finding sampling strategies that attempt to estimate one parameter. Many advantages have been developed around these topics, as seen in previous chapters; however, all of them are motivated under the assumption that the researcher is interested in estimating only one parameter. As Holmberg (2002a) states, “a typical survey in the economic sector involves several characteristics of interest and several target parameters… with multiple parameters of interest and multiple precision requirements, the statistician should choose a sampling design that takes these characteristics into account.”

A survey can be divided into two stages: the design stage and the estimation stage. Anders Holmberg’s work over the past decade has focused on the search for a sampling design that induces unequal inclusion probabilities and is optimal in the sense that it produces a significant increase in precision for each characteristic of interest. This chapter summarizes the sampling design proposals and presents a possible solution to the problem of multiparameter estimation through a matrix approach, providing the reader with a comprehensive approach to joint estimation in sampling. Although the results in this chapter are simple, they offer a powerful tool for planning sampling strategies in multipurpose surveys. First, the estimation approach is proposed for the case of multiple auxiliary information characteristics. Later, Holmberg’s research results are summarized with respect to the sampling design of a strategy involving several characteristics of interest. Finally, the chapter closes with a numerical example that reveals the matrix approach and its advantages in multipurpose surveys.

14.2 Estimation of Several Parameters

Suppose that the survey involves the study of \(Q\) characteristics of interest. Assume that the \(k\)-th element \((k\in U)\) is associated with a vector of \(Q\) characteristics of interest, \(\mathbf{y}_k=(y_{k1},\ldots,y_{kQ})\), whose values are unknown for the finite population. Thus, the following matrix will be called the matrix of interest.

\[ \mathbf{Y}_U=\left( \begin{array}{cccc} y_{11} & y_{12} & \ldots & y_{1Q} \\ \vdots & \vdots & \ddots & \vdots \\ y_{k1} & y_{k2} & \ldots & y_{kQ} \\ \vdots & \vdots & \ddots & \vdots \\ y_{N1} & y_{N2} & \ldots & y_{NQ} \\ \end{array} \right) =\left( \begin{array}{ccccc} \mathbf{y}^1 & \mathbf{y}^2 & \ldots & \mathbf{y}^Q \\ \end{array} \right) \]

Note that the entry \(y_{kq}\) refers to the value of the \(q\)-th characteristic of interest in the \(k\)-th element, with \(k\in U\) and \(q=1,\ldots,Q\). In a design-based inference context, \(\mathbf{y}^q\) is not considered a random vector, because its components are considered fixed but unknown parameters. Thus, the values of each characteristic of interest are not necessarily continuous, such as income, weight, or height, but may also be discrete, such as indicators of population subgroups like domains, strata, or post-strata. In this way, the matrix \(\mathbf{Y}_U\) can be viewed as a matrix of mixed values.

The objective is to estimate the \(Q\) components of the vector of totals defined by the following expression: \[ \mathbf{t}=(t_1,t_2,...,t_Q)'=\mathbf{Y}'_U\mathbf{1}_N, \]

where \(\mathbf{1}_N=(1,1,\ldots,1)_{N\times 1}'\) and \(t_q=\sum_{k\in U}y_{kq}\) is the population total of the \(q\)-th characteristic of interest. When the sample of size \(n\) is selected, \(y_{kq}\) is observed (\(k\in S\)), and the following matrix can be defined: \[ \mathbf{Y}_s=\left( \begin{array}{cccc} y_{11} & y_{12} & \ldots & y_{1Q} \\ \vdots & \vdots & \ddots & \vdots \\ y_{k1} & y_{k2} & \ldots & y_{kQ} \\ \vdots & \vdots & \ddots & \vdots \\ y_{n1} & y_{n2} & \ldots & y_{nQ} \\ \end{array} \right). \]

Note that when \(s=U\), \(\mathbf{Y}_U=\mathbf{Y}_s\). Thus, the matrix of inclusion probabilities is defined by the following expression: \[ \boldsymbol{\Pi}=diag(\pi_1, \pi_2, ..., \pi_n), \]

Along these lines, the Horvitz-Thompson estimator of the vector of totals \(\mathbf{t}\) is defined as \[ \widehat{\mathbf{t}}_{\pi}=(\widehat{t}_{1,\pi}, \widehat{t}_{2,\pi},...,\widehat{t}_{Q,\pi})'= \mathbf{Y}_s'\boldsymbol{\Pi}^{-1}\mathbf{1}_n, \]

with \(\mathbf{1}_N=(1,1,\ldots,1)_{n\times 1}'\), and \(\widehat{t}_{q,\pi}= \sum_{k \in s}y_{kq}/\pi_k\) is the Horvitz-Thompson estimator of \(t_q\). It is easy to prove that \(\widehat{\mathbf{t}}_{\pi}\) is an unbiased estimator of \(\mathbf{t}\), and its variance matrix is given by \[ \mathbf{V}(\widehat{\mathbf{t}}_{\pi})=E(\widehat{\mathbf{t}}_{\pi}-\mathbf{t})(\widehat{\mathbf{t}}_{\pi}-\mathbf{t})'. \]

Note that, if \(N\geq q\), then \(\mathbf{V}(\widehat{\mathbf{t}}_{\pi})\) will be a symmetric positive definite matrix whose \(qq'\) element is \[ \sum_{k\in U}\sum_{l\in U} \Delta_{kl}\frac{y_{kq}}{\pi_k}\frac{y_{lq'}}{\pi_l}, \]

with \(\Delta_{kl}=\pi_{kl}-\pi_k\pi_l\). If \(s\neq U\), it is impossible to calculate the value of the previous expression. However, if \(n\geq q\), the variance can be estimated using a symmetric positive definite matrix \(\widehat{\mathbf{V}}(\widehat{\mathbf{t}}_{\pi})\) whose \(qq'\) element is \[ \sum_{k\in S}\sum_{l\in s} \frac{\Delta_{kl}}{\pi_{kl}}\frac{y_{kq}}{\pi_k}\frac{y_{lq'}}{\pi_l}. \]

In some cases, the survey requirement is the estimation of the vector of population means given by \[ \bar{\mathbf{y}}=\frac{1}{N}\mathbf{t}. \]

Therefore, an unbiased estimator for \(\bar{\mathbf{y}}\) is \[ \bar{\mathbf{y}}_{\pi}=\frac{1}{N}\widehat{\mathbf{t}}_{\pi}, \]

whose variance matrix will be estimated unbiasedly by \(\frac{1}{N^2}\widehat{\mathbf{V}}(\widehat{\mathbf{t}}_{\pi})\). If the population size is unknown, it can be estimated unbiasedly using the principles of the Horvitz-Thompson estimator, such that \[ \widehat{N}_{\pi}=\mathbf{1}_n'\boldsymbol{\Pi}^{-1}\mathbf{1}_n. \]

Note that computational efficiency could be increased by incorporating this matrix approach, because several parameters of interest are estimated through a single algebraic operation.

14.3 Some Sampling Designs

This section introduces some examples of estimating several parameters of interest under the most common sampling designs in the theory.

TipResult

Under the Bernoulli sampling design, the vector of totals \(\mathbf{t}\) is estimated unbiasedly by \[ \widehat{\mathbf{t}}_{\pi}=\frac{1}{\pi}\mathbf{Y}_s'\mathbf{1}_n \]

and its variance matrix is estimated unbiasedly by \[ \widehat{\mathbf{V}}(\widehat{\mathbf{t}}_{\pi})=\frac{1}{\pi}\left(\frac{1}{\pi}-1\right) \mathbf{Y}_s\mathbf{Y}_s'. \]

TipResult

Although the simple random sampling without replacement design is not the most widely used in practice, it is used in the final sampling stages of complex designs. Under this sampling design, \(\mathbf{t}\) is estimated unbiasedly by \[ \widehat{\mathbf{t}}_{\pi}=\frac{N}{n}\mathbf{Y}_s'\mathbf{1}_n. \]

and its covariance matrix is estimated unbiasedly by \[ \widehat{\mathbf{V}}(\widehat{\mathbf{t}}_{\pi})=\frac{N^2}{n}\left(1-\frac{n}{N}\right) \mathbf{S}_y, \]

with \(\mathbf{S}_y\), the covariance matrix of the characteristics of interest calculated with the observations collected in the selected sample. On the other hand, \(\bar{\mathbf{y}}\) is estimated unbiasedly by \[ \bar{\mathbf{y}}_{\pi}=\frac{1}{N}\widehat{\mathbf{t}}_{\pi}=\frac{1}{n}\mathbf{Y}_s'\mathbf{1}_n. \]

and its covariance matrix is estimated unbiasedly by the following expression: \[ \widehat{\mathbf{V}}(\bar{\mathbf{y}}_{\pi})=\frac{1}{N^2}\widehat{\mathbf{V}}(\widehat{\mathbf{t}}_{\pi}). \]

14.3.1 Domain Estimation

If the survey requirements are related to estimating the absolute size of a domain or the total of one or several characteristics of interest in that domain, then the following methodological construction is proposed. Suppose that the population is partitioned into \(D\) domains such that \(U={U_1,\ldots,U_d,\ldots,U_D}\). Then, the domain indicator matrix is defined as \[ \mathbf{Z}=\left( \begin{array}{ccccc} z_{11} & \ldots & z_{1d} & \ldots & z_{1D} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ z_{k1} & \ldots & z_{kd} & \ldots & z_{kD} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ z_{n1} & \ldots & z_{nd} & \ldots & z_{nD} \\ \end{array} \right) \]

where the element \[ z_{kd}= \begin{cases} 1 & \text{if $k\in U_d$, and}\\ 0 & \text{otherwise} \end{cases} \]

The vector of absolute sizes of domain \(d\) is given by \[ \mathbf{N}_{d}=(N_{1},N_{2},...,N_{D})' \]

where \[ N_{d}=\sum_{k\in U}z_{kd}. \]

\(\mathbf{N}_{d}\) is estimated unbiasedly by the Horvitz-Thompson estimator as follows: \[ \widehat{\mathbf{N}}_{d}=(\widehat{N}_{1},\widehat{N}_{2},...,\widehat{N}_{D})'=\mathbf{Z}'\boldsymbol{\Pi}^{-1}\mathbf{1}_n, \]

its variance matrix is estimated unbiasedly by \(\widehat{\mathbf{V}}(\widehat{\mathbf{N}}_{d})\), which is defined analogously to (13.2.6).

In many cases, it is necessary to estimate the totals of characteristics of interest over all domains. Thus, the total of the \(q\)-th variable over all \(D\) domains of interest is given by \[ \mathbf{t}_{dq}=(t_{1q},t_{2q},...,t_{Dq})' \]

and one way to estimate it is given by the following expression: \[ \widehat{\mathbf{t}}_{dq\pi}=(\widehat{t}_{1q\pi},\widehat{t}_{2q\pi},...,\widehat{t}_{Dq\pi})' =(\mathbf{y}^q\mathbf{1}_D\odot\mathbf{Z})'\boldsymbol{\Pi}^{-1}\mathbf{1}_n \]

where \(\mathbf{y}^q\) denotes the \(q\)-th column of the matrix \(\mathbf{Y}_s\), \(\mathbf{1}_D=(1,\ldots,1)'_{D \times 1}\), and \(\odot\) denotes the Hadamard matrix product.

TipResult

Under simple random sampling without replacement, the Horvitz-Thompson estimator for the vector of absolute domain sizes and for the total of the \(q\)-th characteristic of interest across all \(D\) domains are given by \[ \widehat{\mathbf{N}}_{d}=(N/n)\mathbf{Z}'\mathbf{1}_n, \] \[ \widehat{\mathbf{t}}_{dq\pi}=(N/n)(\mathbf{y}^q\mathbf{1}_D\odot\mathbf{Z})\mathbf{1}_n. \] respectively.

14.3.2 Estimation in Stratified Designs

For stratified designs, the following framework is used. The finite population \(U\) is divided into \(H\) mutually exclusive groups or strata \(U_1,\ldots,U_h\ldots,U_H\). Note that before data collection, the membership of each element in each stratum is known. Thus, a random sample is selected in each of the \(H\) strata in the finite population. It is necessary to arrange the matrices to obtain estimates using the principles of the Horvitz-Thompson estimator. Therefore, the matrix \(\mathbf{Y}_s\) is partitioned into \(H\) blocks as follows: \[ \mathbf{Y}_s=\left( \begin{array}{cccc} \mathbf{Y}_{1} \\ \vdots \\ \mathbf{Y}_{h} \\ \vdots \\ \mathbf{Y}_{H} \\ \end{array} \right), \]

where \(\mathbf{Y}_h\) is a submatrix containing the values of each characteristic of interest for the elements belonging to the \(h\)-th stratum, with \(h=1,\ldots,H\). Note that \(\mathbf{Y}_s\in \mathfrak{R}^{Hn\times Q}\) and \(\mathbf{Y}_h\in \mathfrak{R}^{n_h\times Q}\). Defining \(\mathbf{n}=(n_1,\ldots,n_H)'\), then \(n=\mathbf{n}'\mathbf{1}_H=n_1+\cdots+n_H\).

As usual, the objective is to estimate the \(Q\) components of the vector of totals in the \(h\)-th stratum given by \[ \mathbf{t}_h=(t_{1h},t_{2h},...,t_{Qh})'=\mathbf{Y}_h'\mathbf{1}_{N_h}, \]

where \(N_h\) is the size of the \(h\)-th stratum. The population total can be written as \[ \mathbf{t}=(t_1,t_2,...,t_Q)'=\sum_{h=1}^{H}\mathbf{t}_h, \]

where \(\mathbf{t}_h\) is estimated unbiasedly by the following expression: \[ \widehat{\mathbf{t}}_{h\pi}=(\widehat{t}_{1h\pi},\widehat{t}_{2h\pi},...,\widehat{t}_{Qh\pi})'=\mathbf{Y}_h'\mathbf{1}_{n_h}, \]

with \(n_h\) the sample size in the \(h\)-th stratum. Of course, independence is assumed for the sampling design implemented in each stratum. Thus, the population total is given by \[ \widehat{\mathbf{t}}_{\pi}=(\widehat{t}_{1\pi},\widehat{t}_{2\pi},...,\widehat{t}_{Q\pi})'=\sum_{h=1}^{H}\widehat{\mathbf{t}}_h, \]

and its variance matrix can be written as \[ {\mathbf{V}}_{ST}(\widehat{\mathbf{t}}_{\pi})=\sum_{h=1}^{H}{\mathbf{V}}_h(\widehat{\mathbf{t}}_{\pi}) \]

which is estimated unbiasedly by \[ \widehat{{\mathbf{V}}}_{ST}(\widehat{\mathbf{t}}_{\pi})=\sum_{h=1}^{H}\widehat{{\mathbf{V}}}_h(\widehat{\mathbf{t}}_{\pi}). \]

TipResult

Under the stratified random sampling design, the Horvitz-Thompson estimator for the population total is \[ \widehat{\mathbf{t}}_{\pi}=\sum_{h=1}^H\frac{N_h}{n_h}\mathbf{Y}_h'\mathbf{1}_{n_h}, \]

and its covariance matrix is estimated unbiasedly by \[ \widehat{\mathbf{V}}_{STSI}(\widehat{\mathbf{t}}_{\pi})=\sum_{h=1}^H\frac{N_h^2}{n_h}\left(1-\frac{n_h}{N_h}\right)\mathbf{S}_{yh}, \]

with \(\mathbf{S}_{yh}\), the variance matrix of the characteristics of interest in the sample belonging to the \(h\)-th stratum.

14.4 Auxiliary Information

Assume that the \(k\)-th element \((k\in U)\) is associated with a vector of \(P\) auxiliary information characteristics contained in a vector \(\mathbf{x}_k\). The values of this vector \(\mathbf{x}_k=(x_{k1},\ldots,x_{kP})\) are assumed known for the finite population. Thus, we have the following matrix: \[ \mathbf{X}_U=\left( \begin{array}{cccc} x_{11} & x_{12} & \ldots & x_{1P} \\ \vdots & \vdots & \ddots & \vdots \\ x_{k1} & x_{k2} & \ldots & x_{kP} \\ \vdots & \vdots & \ddots & \vdots \\ x_{N1} & x_{N2} & \ldots & x_{NP} \\ \end{array} \right) =\left( \begin{array}{ccccc} \mathbf{x}^1 & \mathbf{x}^2 & \ldots & \mathbf{x}^P \\ \end{array} \right) \]

which will be called the auxiliary information matrix.

14.4.1 Some Relationships

It is possible to assume that there is an explicit linear relationship between each component of the characteristics of interest and the auxiliary information characteristics through a superpopulation model \(\xi_q\), \(q=1,\ldots,Q\), such that \[ \begin{aligned} \underset{(N\times 1)}{\mathbf{Y}^q}&=\underset{(N\times P)}{\mathbf{X}}\underset{(P\times 1)}{\boldsymbol{\beta}^q}+\underset{(N\times 1)}{\mathbf{\boldsymbol{\varepsilon}}^q}. \end{aligned} \]

The model \(\xi_q\) has the following properties: \[ \begin{aligned} E_{\xi_q}(\boldsymbol{\varepsilon}^q)&=\mathbf{0} \\ Var_{\xi_q}(\boldsymbol{\varepsilon}^q)&=\boldsymbol{\Sigma}_q. \end{aligned} \]

\(\Sigma_q\) establishes the variance structure of the vector \(\mathbf{\boldsymbol{\varepsilon}}^q\). Note that the previous relationships can be rewritten through a joint model \(\xi\) such that \[ \begin{aligned} \underset{(N\times Q)}{\mathbf{Y}}&=\underset{(N\times P)}{\mathbf{X}}\underset{(P\times Q)}{\boldsymbol{\beta}}+\underset{(N\times Q)}{\mathbf{\boldsymbol{\varepsilon}}}. \end{aligned} \]

This approach suggests that \(\mathbf{Y}\), \(\mathbf{X}\), and \(\boldsymbol{\varepsilon}\) are random matrices (Gupta and Nagar 1999) defined in the superpopulation model \(\xi\), for which \(\mathbf{Y}_U\) and \(\mathbf{X}_U\) are assumed to be mere realizations of those random matrices. More precisely, the model \(\xi\) has the following characteristics:

\[ \begin{aligned} E_{\xi}(\boldsymbol{\varepsilon})&=\underset{(N\times Q)}{\mathbf{0}} \\ Var_{\xi}(\vec \boldsymbol{\varepsilon})&=\underset{(NQ\times NQ)}{\boldsymbol{\Sigma}}=diag(\Sigma_1,\Sigma_2,\ldots,\Sigma_Q) \end{aligned} \]

Note that the subscript \(\xi\) refers to the expectation under the particular structure induced by that superpopulation model. In practical situations, it is common to assume \(\boldsymbol{\Sigma}_q=\sigma^2_q diag(c_{1q},\ldots,c_{Nq})\), where \(c_{kq}=f_q(x_{k1},\ldots,x_{kP})\) and \(f_q\) is a real-valued function.

The problem of estimating the parameter vector \(\boldsymbol{\beta}\) is considered briefly. Let \(D(\mathbf{X})\) be a translation-invariant dispersion measure such that \(D(\mathbf{X}+\mathbf{K})=D(\mathbf{X})\), with \(\mathbf{K}\) a matrix of constants. Then the estimate of \(\boldsymbol{\beta}\) will correspond to the vector that minimizes the previous dispersion measure. In particular, \(D(\cdot)\) could be given by the total multivariate variance defined as \[ \operatorname{trace}(\mathbf{Y}-\mathbf{X}\boldsymbol{\beta})'(\mathbf{Y}-\mathbf{X}\boldsymbol{\beta}). \tag{14.1}\]

With the previous choice and using the least squares method, (equation 14.1) is minimized by the following expression: \[ \mathbf{B}=(\mathbf{B}_1,\mathbf{B}_2,\ldots,\mathbf{B}_Q), \]

where \[ \mathbf{B}_q=(\mathbf{X}_U'\Sigma_q^{-1}\mathbf{X}_U)^{-1}(\mathbf{X}_U'\Sigma_q^{-1}\mathbf{Y}_U). \]

Note that, in order to calculate this estimate, all population values of both the matrix of characteristics of interest and the auxiliary information characteristics must be known.

14.4.2 Traditional Information

In real applications, only one sample is selected and it is not possible to calculate \(\mathbf{B}\). Therefore, this value must be estimated using the information available in the selected or realized random sample. It can be shown that the following expression corresponds to an asymptotically unbiased estimator for \(\mathbf{B}\): \[ \widehat{\mathbf{B}}=(\widehat{\mathbf{B}_1},\widehat{\mathbf{B}_2},\ldots,\widehat{\mathbf{B}_Q}), \]

where \[ \widehat{\mathbf{B}_q}=(\mathbf{X}_s'\mathbf{A}_q^{-1}\mathbf{X}_s)^{-1}(\mathbf{X}_s'\mathbf{A}_q^{-1}\mathbf{Y}_s), \]

\(q=1,\ldots,Q\), \(\mathbf{X_s}\) is similarly defined as in (13.2.3), and \[ \mathbf{A}_q=\boldsymbol{\Pi}^{1/2}\Sigma_q\boldsymbol{\Pi}^{1/2}. \]

Thus, the multiple general regression estimator for the vector of population totals is defined as \[ \widehat{\mathbf{t}}_{Mgreg}= \widehat{\mathbf{t}}_{\mathbf{y}\pi}+\widehat{\mathbf{B}}'(\mathbf{t}_{\mathbf{x}}-\widehat{\mathbf{t}}_{\mathbf{x}\pi}), \]

with \(\widehat{\mathbf{t}}_{\mathbf{y}\pi}\) and \(\widehat{\mathbf{t}}_{\mathbf{x}\pi}\) the Horvitz-Thompson estimators of \(\mathbf{t_y}\) and \(\mathbf{t}_{\mathbf{x}}\), respectively. Note that \(\widehat{\mathbf{B}}_q\) can also be written as \[ \begin{aligned} \widehat{\mathbf{B}}_q&=(\mathbf{X}_s'\mathbf{D}_{\lambda}\mathbf{X}_s)^{-1}\mathbf{X}_s\mathbf{D}_{\lambda}\mathbf{Y}_s\\ &=\left(\sum_{k \in s}\mathbf{x}_k\lambda_k^q\mathbf{x}_k'\right)^{-1}\left(\sum_{k \in s}\mathbf{x}_k\lambda_k^q\mathbf{y}_k'\right) \end{aligned} \]

where \(\mathbf{D}_{\lambda}=diag(\lambda_1^q,\ldots,\lambda_n^q)\) and \(\lambda_k^q\) are real-valued functions of the inclusion probabilities and the auxiliary information. Note also that the model \(\xi\) serves as a vehicle for finding an appropriate general regression estimator. Once it is found or defined, the model will not be useful for any other sampling purpose. The properties of the multiple general regression estimator, expectation and variance, are also defined from a design-based inference perspective.

14.4.2.1 Some Special Cases

The following scenarios are stated under a general framework and turn out to be special cases of the multiple general regression estimator. In most cases, their particularity is induced by the choice of the values of \(\lambda_k\).

• If \(P=1\), \(\mathbf{x}_k=x_k\), and \(\lambda_k^q=(\pi_kx_k)^{-1}\), then we have the ratio estimator for each characteristic of interest.
• If \(P=2\), \(\mathbf{x}_k=(1, x_k)'\), and \(\lambda_k^q=(\pi_k)^{-1}\), then we have the classical regression estimator.
• If \(P=M (\text{number of post-strata})\), \(\mathbf{x}_k=\delta_k=(0,\ldots0,1,0,\ldots,0)'\), and \(\lambda_k^q=(\pi_k)^{-1}\), where \(\delta_k\) represents \(M\) indicator variables, each indicator representing the membership of the population element in the post-stratum in question, then we have the post-stratification estimator.

Note that the multiple general regression estimator can also be written as follows: \[ \widehat{\mathbf{t}}_{Mgreg}=(\mathbf{W}'\odot\mathbf{Y}_s')\mathbf{1}_n, \]

where \[ \mathbf{W}=\left( \begin{array}{cccc} w_{1}^1 & w_1^2 & \ldots & w_1^Q \\ \vdots & \vdots & \ddots & \vdots \\ w_k^1 & w_k^2 & \ldots & w_k^Q \\ \vdots & \vdots & \ddots & \vdots \\ w_n^1 & w_n^2 & \ldots & w_n^Q \\ \end{array} \right) =\left( \begin{array}{ccccc} \mathbf{w}^1 & \mathbf{w}^2 & \ldots & \mathbf{w}^Q \\ \end{array} \right). \]

We have that \(\mathbf{w}^q=(w_1^q,\ldots,w_k^q,\ldots,w_n^q)'\) is a vector of weights such that \[ w_k^q=\frac{1}{\pi_k}\left(1+\lambda_k^q \mathbf{x}_k'\left(\sum_{k \in s}\mathbf{x}_k\lambda_k^q\mathbf{x}_k'\right)^{-1} (\mathbf{t}_{\mathbf{x}}-\widehat{\mathbf{t}}_{\mathbf{x}\pi}) \right). \]

These weights, as studied in previous chapters, are known as calibration weights, and they exactly reproduce the vector of totals \(\mathbf{t}_{\mathbf{x}}\) when applied to the available auxiliary information. Therefore, \(\mathbf{W}\) is called the calibration matrix. It is not difficult to show that the following relationship \[ \sum_{k\in S}w_k^q\mathbf{x}_k=\mathbf{X}_s'\mathbf{w}^q=\mathbf{t}_{\mathbf{x}}, \]

is satisfied for each \(q=1,\ldots,Q\). It is interesting to observe that \(\mathbf{t}_{\mathbf{x}}\) is calibrated under different choices of the weights \(w^q\). On the other hand, note that \[ \mathbf{w}^q=\boldsymbol{\Pi}^{-1}\mathbf{1}_n+\mathbf{A}_q\mathbf{X}_s\left(\mathbf{X}_s'\mathbf{A}_q\mathbf{X}_s\right)^{-1} (\mathbf{t}_{\mathbf{x}}-\widehat{\mathbf{t}}_{\mathbf{x}\pi}) \]

When dealing with post-stratified estimation, a generalized inverse must be used, relying on the property that the multiple general regression estimator is invariant to any inverse.

14.4.3 Joint Auxiliary Information

The least squares method is not the only path to obtaining a multiple general regression estimator. In this section, the existence of a joint information matrix is assumed, whose algebraic structure is defined by the following expression: \[ \mathbf{V}=\left( \begin{array}{cccccccc} y_{11} & y_{12} & \ldots & y_{1Q} & x_{11} & x_{12} & \ldots & x_{1P}\\ y_{21} & y_{22} & \ldots & y_{2Q} & x_{21} & x_{22} & \ldots & x_{2P}\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots\\ y_{n1} & y_{n2} & \ldots & y_{nQ} & x_{n1} & x_{n2} & \ldots & x_{nP}\\ \end{array} \right). \]

The estimator of the vector of totals of the characteristics of interest and the auxiliary information characteristics is given by \(\widehat{\mathbf{t}}_{\mathbf{v}\pi}\), which is defined as \[ \widehat{\mathbf{t}}_{\mathbf{v}\pi}=\mathbf{V}'\boldsymbol{\Pi}^{-1}\mathbf{1}_n. \]

Thus, suppose that \(\widehat{\mathbf{t}}_{\mathbf{v}\pi}\) follows a multivariate normal distribution with mean \[ E\left(\widehat{\mathbf{t}}_{\mathbf{v}\pi}\right) =(\mathbf{t}'_{\mathbf{Y}\pi},\mathbf{t}'_{\mathbf{X}\pi})'=\mathbf{t}_{\mathbf{v}}, \]

and variance matrix defined as \[ V\left(\widehat{\mathbf{t}}_{\mathbf{v}\pi}\right)=\left( \begin{array}{cc} V(\mathbf{\widehat{t}}_{\mathbf{y}\pi}) & C(\mathbf{\widehat{t}}_{\mathbf{y}\pi},\mathbf{\widehat{t}}_{\mathbf{x}\pi})\\ C(\mathbf{\widehat{t}}_{\mathbf{y}\pi},\mathbf{\widehat{t}}_{\mathbf{x}\pi}) & V(\mathbf{\widehat{t}}_{\mathbf{x}\pi})\\ \end{array}, \right), \]

where \(V(\mathbf{\widehat{t}}_{\mathbf{y}\pi})\) is considered a symmetric matrix such that the \(j\)-th element of its diagonal is given by the variance of \(\widehat{t}_{y_j\pi}\) \[ V(\widehat{t}_{y_j\pi})=\sum\sum_U\Delta_{kl}\frac{y_{jk}}{\pi_k}\frac{y_{jl}}{\pi_l}, \]

and the \(ij\) element, outside its diagonal, is given by the covariance of \(\widehat{t}_{y_i\pi}\) and \(\widehat{t}_{y_j\pi}\), \[ C(\widehat{t}_{y_i\pi},\widehat{t}_{y_j\pi})=\sum\sum_U\Delta_{kl}\frac{y_{ik}}{\pi_k}\frac{y_{jl}}{\pi_l}. \]

\(V(\mathbf{\widehat{t}}_{\mathbf{x}\pi})\) is defined analogously, and \(C(\mathbf{\widehat{t}}_{\mathbf{y}\pi},\mathbf{\widehat{t}}_{\mathbf{x}\pi})\), not necessarily symmetric, is such that the \(ij\) element is given by the covariance of \(\widehat{t}_{y_i\pi}\) and \(\widehat{t}_{x_j\pi}\) \[ C(\widehat{t}_{y_i}\pi,\widehat{t}_{x_j}\pi)=\sum\sum_U\Delta_{kl}\frac{y_{ik}}{\pi_k}\frac{x_{jl}}{\pi_l}. \]

Following the results of multivariate inference for populations with normal distribution, the conditional distribution of \(\widehat{\mathbf{t}}_{\mathbf{Y}\pi}\) given \(\widehat{\mathbf{t}}_{\mathbf{X}\pi}\) also follows a multivariate normal distribution with conditional mean given by \[ E(\widehat{\mathbf{t}}_{\mathbf{y}\pi}|\widehat{\mathbf{t}}_{\mathbf{x}\pi})=\mathbf{t}_{\mathbf{y}\pi}+ C(\mathbf{\widehat{t}}_{\mathbf{y}\pi},\mathbf{\widehat{t}}_{\mathbf{x}\pi})(V(\mathbf{\widehat{t}}_{\mathbf{x}\pi}))^{-1}(\mathbf{t}_{\mathbf{x}} -{\widehat{\mathbf{t}}}_{\mathbf{x}\pi}), \tag{14.2}\]

and conditional variance given by \[ V(\widehat{\mathbf{t}}_{\mathbf{y}\pi}|\widehat{\mathbf{t}}_{\mathbf{x}\pi})=V(\mathbf{\widehat{t}}_{\mathbf{y}\pi})- C(\mathbf{\widehat{t}}_{\mathbf{y}\pi},\mathbf{\widehat{t}}_{\mathbf{x}\pi}) (V(\mathbf{\widehat{t}}_{\mathbf{x}\pi}))^{-1}C(\mathbf{\widehat{t}}_{\mathbf{x}\pi},\mathbf{\widehat{t}}_{\mathbf{y}\pi}). \tag{14.3}\]

Note that (equation 14.2) and (equation 14.3) are estimated unbiasedly by \[ \begin{aligned} \widehat{\mathbf{t}}_{\mathbf{y}}&=\widehat{\mathbf{t}}_{\mathbf{y}\pi}+ \widehat{C}(\mathbf{\widehat{t}}_{\mathbf{y}\pi},\mathbf{\widehat{t}}_{\mathbf{x}\pi}) (\widehat{V}(\mathbf{\widehat{t}}_{\mathbf{x}\pi}))^{-1}(\mathbf{t}_{\mathbf{x}}-{\widehat{\mathbf{t}}}_{\mathbf{x}\pi})\\ &=\widehat{\mathbf{t}}_{\mathbf{y}\pi}+\widehat{\mathbf{B}}(\mathbf{t}_{\mathbf{x}}-{\widehat{\mathbf{t}}}_{\mathbf{x}\pi}) \end{aligned} \]

y, \[ \widehat{V}(\widehat{\mathbf{t}}_{\mathbf{y}})=V(\mathbf{\widehat{t}}_{\mathbf{Y}\pi})-C(\mathbf{\widehat{t}}_{\mathbf{Y}\pi}, \mathbf{\widehat{t}}_{\mathbf{X}\pi})(V(\mathbf{\widehat{t}}_{\mathbf{X}\pi}))^{-1} C(\mathbf{\widehat{t}}_{\mathbf{X}\pi},\mathbf{\widehat{t}}_{\mathbf{Y}\pi}), \]

respectively. On the other hand, observe that (13.4.22) looks like the multiple general regression estimator. However, its slope, \(\widehat{\mathbf{B}}\), would be different: while the slope of the general regression estimator is given by the least squares method, the slope of this latter estimator corresponds, according to the results of multivariate statistical inference, to a set of multiple regressions of \(\mathbf{X}\) on \(\mathbf{Y}\). This estimator of the vector of totals of the characteristic of interest should be called the optimal general regression estimator and has been studied by Cassady and Valliant (1993) in the context of model-based inference for estimating the total of a single characteristic of interest.

14.5 Optimal Sampling Designs

This section addresses the problem of sample selection under a unified criterion that accounts for the structural behavior of each characteristic of interest. That is, in the design stage of a multipurpose survey, a comprehensive sampling design must be chosen, and for this the Holmberg approach will be considered. Thus, it can be assumed that, in the planning stage of the strategy, auxiliary information characteristics may be involved, making it possible to take some positions on the validity of the statistical relationships between the characteristics of interest and the auxiliary information variables.

14.5.1 Holmberg’s Sampling Design

Suppose that the characteristics of interest involved in the survey all have the same importance¹. Under this approach, a brief summary of the Holmberg design used in multipurpose surveys is presented below:

1. For each characteristic of interest, the statistician, researcher, or end user must propose a sampling design, \(p_q(\cdot)\) \((q=1,\ldots,Q)\), that is optimal and such that the expected sample size is \(E(n(S))=n_q\). Of course, note that each of the \(Q\) sampling designs may be different. Moreover, the sample sizes in each proposed design do not necessarily have to be equivalent. Recall that in the traditional approach, which is not concerned with including several characteristics of interest, the statistician must propose a single sampling design, which is assumed to be optimal for all parameters to be estimated.

2. Each sampling design \(p_q(\cdot)\) induces a vector of inclusion probabilities of size \(N\) for each of the elements belonging to the finite population. These inclusion probabilities must take the following form (Holmberg 2002b, eq. 6): \[ \pi_{qk}=n_q\frac{\sigma_{qk}}{\sum_{k\in S}\sigma_{qk}}, \]

with \(\sigma_{qk}\) size measures, usually though not necessarily linked to a linear regression model. The “optimal sampling design” feature is obtained if \(\pi_{qk}\propto\sigma_{qk}\). Note that if the optimal sampling design for the \(q\)-th characteristic of interest is simple random sampling without replacement, then \(\sigma_{qk}=1\) for all \(k\in U\). On the other hand, with the choice \(\sigma_{qk}^2=\sigma^2_qx_{qk}^{\gamma_q}\), where \(\sigma^2_q\) is a constant and \(x_{qk}\) corresponds to the value of the \(k\)-th element for some auxiliary variable, or to a function of many auxiliary information variables, then the optimal sampling design must be proportional to the size of \(\sigma_{qk}\) (\(\pi\)PS). That is, \(\pi_{qk}\propto x_{qk}^{\gamma_q/2}\).

1. Based on the minimum average loss of relative efficiency criterion, ANOREL, the optimal sample size for the multipurpose survey will be given by \[ n^*\geq \frac{(\sum_{k\in U}\sqrt{a_{qk}})^2}{(1+c)Q+\sum_{k\in U}a_{qk}}, \]

where \[ a_{qk}= \sum_{q=1}^Q \frac{\sigma^2_{qk}}{\sum_{k\in U}\left( \frac{1}{\pi_{qk}}-1\right)\sigma^2_{qk}}, \]

and \(c\) is the maximum error allowed under the ANOREL criterion, on a scale from zero to one. Note that in practice, \(\sigma^2_{qk}\) is unknown and must be written as a function of the auxiliary information variables. Holmberg (2002b) states that subjective knowledge, experience, or external sources can be used to obtain approximations to the exact value of this quantity.

1. Once the sample size has been calculated, a single vector of inclusion probabilities that is optimal for all characteristics of interest must be created. This vector is induced by Holmberg’s sampling design, which minimizes the average loss of relative efficiency, and is given by the following expression: \[ \pi_{(opt)k}=\frac{n^*\sqrt{a_{qk}}}{\sum_{k\in U}\sqrt{a_{qk}}} \]

2. In most cases, the resulting inclusion probability vector, \(\boldsymbol{\pi}_{(opt)}=(\pi_{(opt)1},\ldots,\pi_{(opt)N})'\), is a vector of unequal inclusion probabilities. In this situation, a \(\pi\)PT sample selection scheme must be used.

14.5.2 A Numerical Example

In this section, an example of the multipurpose approach is considered. At the design stage, an optimal sampling design is chosen using the Holmberg approach (Holmberg 2002b), and at the estimation stage the matrix approach (Gutiérrez 2009) is implemented. Both stages are carried out using the R software. In particular, the sampling package is introduced for sample selection and estimation in several domains of interest.

For this purpose, a real population is considered, the population of Swiss municipalities MU281 available in Appendix B of Särndal et al. (1992). Thus, it is possible to plan a multipurpose survey in which the characteristics of interest and the domains of interest are provided in advance, and in which it is possible to have some beliefs about the structural behavior of the population and the relationship between the characteristics of interest and the auxiliary information characteristics. Note that the goal is not to present a perfect sampling design, but rather one that illustrates the practical development of the theory in a multipurpose survey. The characteristics of interest are:

\(y_1\)	=	P85	(Population in 1985)
\(y_2\)	=	RMT85	(Taxes accrued by municipalities in 1985)
\(y_3\)	=	REV84	(Real estate values in 1984)

The auxiliary information characteristics are:

\(x_1\)	=	P75	(Population in 1975)
\(x_2\)	=	S82	(Number of seats in the municipal council in 1982)

For domain estimation, the following variable is used:

\(z\)	=	REG	(geographic region indicator)

The following computational code was used to specify the survey characteristics.

library(sampling)
data(MU284)
MU281 <- MU284[MU284$RMT85 <= 3000, ]
attach(MU281)

Y1 <- P85
Y2 <- RMT85
Y3 <- REV84
X1 <- P75
X2 <- S82
Z <- REG

To have some degree of certainty about the quality of the estimation, the totals of the characteristics of interest and the auxiliary information characteristics are available.

Ty <- c(sum(Y1), sum(Y2), sum(Y3))
Tx <- c(281, sum(X1), sum(X2))

Ty

[1]   7033  53151 757246

Tx

[1]   281  6818 13257

Now, assuming that the importance of the three characteristics of interest is the same, the Holmberg approach for this particular case is described below:

1. In population MU281, the population size is \(N=281\). Suppose that the statistician considers that, for each of the three characteristics of interest, the sample size should be equal to 100.

N <- 281
n <- 100

1. Assume that, using knowledge from external sources, the statistician assumes that the best sampling designs, in the optimal sense, are: for \(y_1\), a \(\pi\)PT sampling design with \(\pi_{1k}\propto x_{1k}^{0.7}\); for \(y_2\), a \(\pi\)PT sampling design with \(\pi_{2k}\propto x_{1k}\); and finally, for \(y_3\), a simple random sampling design.

sigy1 <- sqrt(X1^(1.4))
sigy2 <- sqrt(X1^(2))
sigy3 <- rep(1, N)

pik1 <- n * sigy1 / (sum(sigy1))
pik2 <- n * sigy2 / (sum(sigy2))
pik3 <- n * sigy3 / (sum(sigy3))

1. The optimal sample size based on the ANOREL criterion for this multiparameter case would be \(n^*=108\). The following code verifies this.

a1 <- sigy1^2 / (sum(((1 / pik1) - 1) * sigy1^2))
a2 <- sigy2^2 / (sum(((1 / pik2) - 1) * sigy2^2))
a3 <- sigy3^2 / (sum(((1 / pik3) - 1) * sigy3^2))
aqk <- a1 + a2 + a3

n.st <- ((sum(sqrt(aqk)))^2) / ((1 + 0.03) * 3 + (sum(aqk)))
n.st <- as.integer(n.st)
n.st

[1] 108

1. The vector of optimal inclusion probabilities for the three characteristics of interest is given by the following code. Note that their sum in the population is equal to the sample size.

pikopt <- n.st * sqrt(aqk) / sum(sqrt(aqk))
sum(pikopt) == n.st

[1] TRUE

1. Because the entries of the resulting inclusion probability vector are unequal, the sample must be selected with some ordered sampling design, with unequal inclusion probabilities and fixed sample size. The UPopips function from the sampling package selects a sample with these characteristics. Once the sample is selected, the getdata function is used to extract the observed data.

sam <- UPopips(pikopt, "exponential")
sel <- sam == 1
head(getdata(MU281, sam))

   LABEL P85 P75 RMT85 CS82 SS82 S82 ME84 REV84 REG CL
5      5  56  52   536   20   27  61 3951  5183   1  1
6      6  16  15   134   16   12  41  918  2157   1  2
7      7  70  62   623   18   27  61 4367  7072   1  2
8      8  66  54   517   15   32  61 4345  5246   1  2
11    11  32  29   277   14   20  45 1993  3264   1  3
12    12  20  14   155   10   21  41 1312  1899   1  3

When the sample is selected, the statistician faces the problem of multiparameter estimation for the characteristics of interest. It is possible to write computational code to estimate the parameters of interest in the traditional way, or to write the computational code once, based on the matrix approach. For the example of population MU281, for which optimal inclusion probabilities \(\pi_{(opt)k}\) were obtained, the Horvitz-Thompson estimator for the vector of totals of the characteristics of interest, for the vector of totals of the auxiliary information characteristics, and for the population size is calculated using the following code.

Ys <- cbind(Y1, Y2, Y3)[sel, ]
Xs <- cbind(1, X1, X2)[sel, ]
PI <- diag(pikopt[sel])
ones <- rep(1, sum(sel))

TyHT <- t(Ys) %*% solve(PI) %*% ones
TxHT <- t(Xs) %*% solve(PI) %*% ones
NHT <- t(ones) %*% solve(PI) %*% ones

The result of running the previous code is a vector of estimated totals. In particular, the estimate of the totals of the characteristics of interest is given by

TyHT

     [,1]
Y1   5802
Y2  43001
Y3 691476

If one or more domains of interest are involved in the estimation stage, the matrix approach provides a simple, comprehensive, and effective estimation method. The domain of interest for this particular case corresponds to the variable REG, which contains 8 geographic categories. Thus, it is possible to obtain estimates of the parameters of interest disaggregated by these regions. Using the disjunctive function from the sampling package, it is possible to create the indicator matrix for the domains given in (13.3.7) and obtain the estimates corresponding to (13.3.11) and (13.3.13).

Z <- disjunctive(Z)[sel, ]
NdHT <- t(Z) %*% solve(PI) %*% ones
Ty1d <- t(Ys[, 1] * Z) %*% solve(PI) %*% ones
Ty2d <- t(Ys[, 2] * Z) %*% solve(PI) %*% ones
Ty3d <- t(Ys[, 3] * Z) %*% solve(PI) %*% ones

It is also possible to gather the estimation results using a simple data table given by:

TydHT <- data.frame(NdHT, Ty1d, Ty2d, Ty3d)
TydHT

  NdHT Ty1d Ty2d   Ty3d
1 26.1 1004 8475 107244
2 37.9  712 5207  89820
3 42.6  720 4828  77898
4 38.3  904 6251  94350
5 42.1  889 6226 115950
6 47.7  949 7152  97457
7  3.7  141 1009  13492
8 29.4  483 3853  95264

If the statistician suspects that a superpopulation-model-assisted inference approach can be used, then the relationships between the auxiliary information characteristics and the characteristics of interest must be established through a model. In this particular example, there are three models, \(\xi_q\) \((q=1,2,3)\), involved in a general model \(\xi\). The relationship is as dictated by the following expression: \[ Y_q=\beta_{q0}+\beta_{q1}X_1+\beta_{q2}X_2+\boldsymbol{\varepsilon}_i \ \ \ \ \ \ \ q=1,2,3. \]

Note that \(E_{\xi_i}(\boldsymbol{\varepsilon}_i)=0\) and that the variance structure of the previous models is induced by step number two of Holmberg’s design, which in particular is given by \[ \begin{aligned} \Sigma_1&=\sigma^2_1 diag(x_{11},x_{12},\ldots,x_{1N})^{1.4}\\ \Sigma_2&=\sigma^2_2 diag(x_{11},x_{12},\ldots,x_{1N})^{2}\\ \Sigma_3&=\sigma^2_3 \mathbf{I}_{N\times N} \end{aligned} \]

Then, the general model takes the following form:

\[ \begin{aligned} \left( \begin{array}{ccc} Y_{11} & Y_{21} & Y_{31} \\ Y_{12} & Y_{22} & Y_{32} \\ \vdots & \vdots & \vdots \\ Y_{1N} & Y_{2N} & Y_{3N} \\ \end{array} \right) &= \left( \begin{array}{ccc} 1 & X_{11} & X_{21} \\ 1 & X_{12} & X_{22} \\ \vdots & \vdots & \vdots \\ 1 & X_{1N} & X_{2N} \\ \end{array} \right) \left( \begin{array}{ccc} \beta_{10} & \beta_{20} & \beta_{30} \\ \beta_{11} & \beta_{21} & \beta_{31} \\ \beta_{12} & \beta_{22} & \beta_{32} \\ \end{array} \right)\notag\\ &+ \left( \begin{array}{ccc} \varepsilon_{11} & \varepsilon_{21} & \varepsilon_{31} \\ \varepsilon_{12} & \varepsilon_{22} & \varepsilon_{32} \\ \vdots & \vdots & \vdots \\ \varepsilon_{1N} & \varepsilon_{2N} & \varepsilon_{3N} \\ \end{array} \right) \end{aligned} \]

Thus, the estimation of the matrix of regression coefficients in the finite population, which involves the variance structure of each model given in (13.4.7), is calculated using the following code:

A1 <- diag(pikopt[sel] * Xs[, 2]^(1.4))
B1 <- (solve(t(Xs) %*% A1 %*% Xs)) %*% (t(Xs) %*% A1 %*% Ys[, 1])
A2 <- diag(pikopt[sel] * Xs[, 2]^(2))
B2 <- (solve(t(Xs) %*% A2 %*% Xs)) %*% (t(Xs) %*% A2 %*% Ys[, 2])
A3 <- diag(pikopt[sel])
B3 <- (solve(t(Xs) %*% A3 %*% Xs)) %*% (t(Xs) %*% A3 %*% Ys[, 3])

B <- matrix(c(B1, B2, B3), ncol = 3, nrow = 3)
B

       [,1]  [,2] [,3]
[1,]  2.562 148.5 1717
[2,]  1.085   9.7  105
[3,] -0.076  -3.7  -28

The next step is to implement the multiple general regression estimator for the totals of interest given in (13.4.10). The computational code requires only one line to perform the calculation, as shown below.

TyMgreg <- TyHT + t(B) %*% (Tx - TxHT)
TyMgreg

     [,1]
Y1   7090
Y2  53003
Y3 814446

This estimator can take different forms. Among others, it can be rewritten in simplified form as in (13.4.13). However, it is first necessary to calculate the calibration matrix given by expression (13.4.14). The following code shows the implementation of the theory.

w1 <- solve(PI) %*% ones + (A1 %*% Xs) %*% (solve(t(Xs) %*% A1 %*% Xs)) %*% (Tx - TxHT)
w2 <- solve(PI) %*% ones + (A2 %*% Xs) %*% (solve(t(Xs) %*% A2 %*% Xs)) %*% (Tx - TxHT)
w3 <- solve(PI) %*% ones + (A3 %*% Xs) %*% (solve(t(Xs) %*% A3 %*% Xs)) %*% (Tx - TxHT)
W <- cbind(w1, w2, w3)
TyMgreg <- t(W * Ys) %*% ones
TyMgreg

     [,1]
Y1   7090
Y2  53003
Y3 814446

The calibration principle shown in (13.4.16) can be easily verified for each column of the calibration matrix. In particular, for the second column, the result holds.

t(w2) %*% Xs

           X1    X2
[1,] 281 6818 13257

In this way, we have shown how to plan and develop a multipurpose survey: first, by using Holmberg’s sampling design at the design stage, and then the matrix approach at the estimation stage.

14.6 Frame and Lucy

When planning a survey, the traditional approach focuses on a single variable of interest, which is insufficient for the statistician who must estimate several parameters of interest. In this chapter, and indirectly throughout the entire book, a useful approach was proposed for the simultaneous estimation of several parameters of interest. In addition to its computational advantages, this matrix approach serves as a vehicle for introducing advanced sampling topics such as the general weighting system proposed by Lavallé and Caron (2008).

Of course, this chapter closes with Frame and Lucy, who, indirectly throughout the book, have shown that the matrix approach for simultaneous estimation should be used by both the theoretical and practical statistician. Suppose that the sampling frame has the quality of providing, in addition to the identification and location of each firm, an auxiliary information characteristic such as each firm’s Income. Along these lines, the reader who has followed the book directly up to this point will know that the population size is \(N=2396\) and that excellent results have been obtained with sampling designs proportional to the size of firm Income for the characteristics of interest Employees and Income. Moreover, these good results have been obtained with a sample size of \(n=2000\).

On the other hand, suppose that the relationship between the auxiliary information characteristic Income is linear for the characteristic of interest Employees but quadratic for the characteristic of interest Taxes. These characteristics must be defined in the computational environment, thereby determining the quantities \(\sigma_{qk}\) in expression (13.5.1), as follows.

data(BigLucy)
attach(BigLucy)
N <- nrow(BigLucy)
n <- c(2000, 2000)
sigy1 <- sqrt(Income^(1))
sigy2 <- sqrt(Income^(2))
sigma <- cbind(sigy1, sigy2)

Using the PikHol function from the TeachingSampling package, which contains three computational parameters: the first, n, is a vector of sample sizes according to the optimality of each design for each variable of interest involved in the survey; sigma is a matrix with \(N\) rows and as many columns as characteristics of interest, in which each of the quantities \(\sigma_{qk}\) determining the relationships between the characteristics of interest and the auxiliary information are stored; and finally, e, which corresponds to the maximum error allowed under the ANOREL criterion. The function returns a vector of optimal inclusion probabilities for all individuals in the finite population, whose sum gives the optimal sample size under this ANOREL criterion.

pis <- PikHol(n, sigma, e = 0.03)
n.Hol <- ceiling(sum(pis))
n.Hol

[1] 1984

From here on, everything becomes familiar because the S.piPS function from the TeachingSampling package is used to select a random sample of firms. The result of this function is, on one hand, a vector containing the realized sample and, on the other, a vector of inclusion probabilities for the selected firms. After data collection, the E.piPS function is used to obtain the estimates that are optimal under the ANOREL criterion.

res <- S.piPS(n.Hol, pis)
sam <- res[, 1]
Pik.s <- res[, 2]
sample_data <- BigLucy[sam, ]
attach(sample_data)

target_variables <- data.frame(Income, Employees, Taxes)
E.piPS(target_variables, Pik.s)

                     N      Income Employees     Taxes
Estimation     85479.5 36593858.47 5434258.1 1009075.2
Standard Error   932.0   121199.96   60380.9   16704.7
CVE                1.1        0.33       1.1       1.7
DEFF               Inf        0.04       0.7       0.1

estimates <- E.piPS(target_variables, Pik.s)
knitr::kable(estimates)

Table 14.1: Holmberg sampling: estimation of the totals of the characteristics of interest

	N	Income	Employees	Taxes
Estimation	85479.5	36593858.47	5434258.2	1009075.2
Standard Error	932.0	121199.96	60380.9	16704.7
CVE	1.1	0.33	1.1	1.7
DEFF	Inf	0.04	0.7	0.1

14.7 Exercises

1. Prove the following equality: \[ \begin{aligned} Cov(\hat{t}_{y,\pi}, \hat{t}_{x,\pi})=\sum_U\sum_U\Delta_{kl}\frac{y_k}{\pi_k}\frac{x_l}{\pi_l} \end{aligned} \]

2. Prove that, for a simple random sampling design, the following relationship holds: \[ \begin{aligned} Cov(\hat{t}_{y,\pi}, \hat{t}_{x,\pi})=\frac{N^2}{n}\left(1-\frac{n}{N}\right)Cov_S(y,x) \end{aligned} \] where \[ \begin{aligned} Cov_S(y,x)=\frac{1}{n-1}\sum_S(y_k-\bar{y}_S)(x_k-\bar{x}_S) \end{aligned} \]

3. Obtain an expression for Holmberg’s inclusion probabilities when all size measures are constant in a multipurpose study.

4. Prove expression (13.3.13).

Cassady, R., and R. Valliant. 1993. “Condicional Properties of Poststratified Estimators Under Normal Theory.” Survey Methodology 19: 183–92.

Gupta, A. K., and D. K. Nagar. 1999. Matriz Variate Distributions. Chapman; Hall.

Gutiérrez, H. A. 2009. “Estimation Stage in Survey Sampling: A Multiparameter Approach.” Revista Colombiana de Estadística 32: 79–97.

Holmberg, A. 2002a. “A Multiparameter Perspective on the Choice of Sampling Design in Surveys.” Statistics in Transition 5: 969–94.

Holmberg, A. 2002b. On the Choice of Sampling Design Under GREG Estimation in Multiparameter Surveys. No. 1. RD Department, Statistics Sweden.

Lavallé, P., and P. Caron. 2008. “Estimation Using the Generalized Weigth Share Method: The Use of Record Linkage.” Survey Methodology 27: 155–69.

Särndal, C. E., B. Swensson, and J. Wretman. 1992. Model Assisted Survey Sampling. Springer, New York.

Smith, T. M. F. 1976. “The Foundations of Survey Sampling: A Review (with Discussion).” Journal of the Royal Statistical Society 139: 183–204.

---

1. Of course, variants of this assumption are possible and can be consulted in Holmberg (2002b). However, in this chapter it will be assumed that the survey gives equal importance to all characteristics of interest.