U <- c("Yves", "Ken", "Erik", "Sharon", "Leslie")
U[1][1] "Yves"U[2][1] "Ken"2 Probability Samples and Estimators
The mathematical basis for developing the sampling model lies in the theory of statistical inference and, more directly, in the application of the basic principles of probability theory. The results of the sampling model are valid only if one starts from the certainty of having a sample that satisfies the conditions required by statistical inference.
Bautista (1998)
2.1 Population and Random Sample
The process of estimation and inference in finite populations, which are ultimately the populations we easily encounter in reality and on which sampling focuses, is very different from the inferential process of classical statistics. Classical statistics treats observed values as realizations of a random variable. In contrast, sampling assumes that the observed values correspond to fixed population parameters. Starting from this fact, let us formalize some concepts that are of vital importance in the study and analysis of sampling.
2.1.1 Finite Population
ImportantDefinition
A finite population is a set of \(N\) elements \(\{e_1, e_2, ..., e_N\}\). Each unit can be identified unambiguously by a set of labels. Let \(U=\{1,2,...,N\}\) be the set of labels of the finite population. The population size is not necessarily known.
It is the set of \(N\), where \(N<\infty\), units that make up the universe of study. \(N\) is commonly called the population size. Each element belonging to the population can be identified by a label. Let \(U\) be the set of labels, such that
\[\begin{equation*} U=\{1,...,k,...,N\}. \end{equation*}\]
The subscript \(k\) will be used to denote the physical existence of the \(k\)-th element. Note that the population size, \(N\), is not always known, and in some cases the objective of the investigation is to estimate it.
2.1.2 Random Sample
It is a subset of the population that has been drawn by means of a statistical selection mechanism. We will denote the random sample with an uppercase letter \(S\) and one of its realizations with a lowercase letter \(s\). Thus, without ambiguity, a selected (realized) sample is the set of units belonging to
\[\begin{equation*} s=\{1,...,k,...,n(S)\}. \end{equation*}\]
The number of components of \(s\) is called the sample size and is not always fixed. That is, in some cases \(n(S)\) is a random quantity. The set of all possible samples is known as the support. Drawing an analogy with classical statistical inference, the support generated by a random sample corresponds to the sample space generated by a random variable.
The preceding definition of a sample, in which the included elements are listed within a set, corresponds to the classical notation. However, a sample can also be denoted as a vector of size \(N\). In this way, the \(k\)-th entry of the vector will denote the number of times the element was included or selected; if the value is zero, it indicates that the element was not included in the selected sample; if the value is different from zero, it indicates that the element was selected. Although both forms of notation have the same interpretation, to avoid confusion the sample in vector form will be denoted by a bold \(\mathbf{s}\), while the sample in set form will be denoted by a plain, nonbold \(s\). More precise definitions of a random sample with or without replacement are given below.
2.1.2.1 Random Sample Without Replacement
ImportantDefinition
A sample without replacement is denoted by a column vector \[\begin{equation} \mathbf{s}=(I_1,I_2,...,I_N)' \in \{0,1\}^N \end{equation}\] where \[\begin{equation} I_k= \begin{cases} 1 & \text{if the $k$th element belongs to the sample,}\\ 0 & \text{otherwise} \end{cases} \end{equation}\]
A random sample is said to be without replacement if each element is included from among the elements that have not yet been chosen; in this way the set \(s\) will never have repeated elements. The sample size corresponds to the cardinality of \(s\). \[\begin{equation} n(S)=\sum_{k \in U}I_k. \end{equation}\]
Because \(n(S)\) is not a fixed quantity, one of the following scenarios may occur: a) the sample contains no elements, in which case the sample is said to be empty; b) the sample contains all elements of the population, in which case it is known as a census.
2.1.2.2 Random Sample With Replacement
ImportantDefinition
A sample with replacement is denoted by a column vector \[\begin{equation} \mathbf{s}=(n_1,n_2,...,n_N)' \in \mathbb{N}^N \end{equation}\] where \(n_k\) is the number of times element \(k\) appears in the sample
In some cases, for the convenience of the selection mechanism, the user prefers to draw a random sample with replacement if the inclusion of each element considers all elements, whether or not they have already been chosen to belong to the sample. In this way, the user can select a sample whose selection process includes an individual \(m\) times (note that \(m\) may be greater than \(N\)). However, in a random sample with replacement, two or more components may be identical. An element that is included more than once in \(s\) is called a repeated element.
In principle, the sample size is given by \[\begin{equation} n(S)=m=\sum_{k \in U}n_k. \end{equation}\]
The number of distinct elements in a random sample \(S\) with replacement is called the effective sample size and, with probability one, is less than or equal to \(N\).
2.1.3 Sampling Supports
In the following chapters, the specific treatment of particular sampling strategies will begin; that is, sampling designs that fit certain situations and estimators that improve the efficiency of the strategy. However, before proceeding, it is necessary for the reader to understand that sampling strategies are defined in terms of the type of sampling used for sample selection. In general, there are two basic distinctions.
- Type of sampling: selection of units with replacement or without replacement.
- Sample size: fixed or random sample size.
As will be seen in later chapters, the sampling strategy, the theoretical treatment for parameter estimation, and the type of support are defined according to the preceding conditions. This section deals specifically with the different forms that the support of a sampling design can take depending on the two basic distinctions. To begin, it is necessary to state the following definitions.
ImportantDefinition
A support \(Q\) is a set of samples.
ImportantDefinition
A support is called symmetric if, for any \(s \in Q\), all permutations of \(s\) also belong to \(Q\).
In the following chapters, unless otherwise stated, the term support will refer to a symmetric support. Some particular symmetric supports are:
The symmetric support without replacement defined as \[\begin{equation*} \mathcal{S}=\{0,1\}^N \end{equation*}\] Note that \[\begin{equation*} \#(\mathcal{S})=2^N \end{equation*}\] For example, if \(N=3\), then \(\mathcal{S}\) is defined by the following samples: { \[\begin{equation*} \mathcal{S}=\{(0,0,0)',(1,0,0)',(0,0,1)',(1,0,1)',(0,1,0)',(1,1,0)',(0,1,1)',(1,1,1)'\} \end{equation*}\] }
The fixed-size symmetric support without replacement defined as \[\begin{equation*} \mathcal{S}_n=\left\{ \textbf{s} \in \mathcal{S}| \sum_{k\in U}s_k=n \right\} \end{equation*}\] Note that \[\begin{equation*} \#(\mathcal{S}_n)=\binom{N}{n} \end{equation*}\] For example, if \(N=3\) and \(n=2\), then \(\mathcal{S}_n\) is defined by the following samples: \[\begin{equation*} \mathcal{S}_n=\{(1,0,1)',(1,1,0)',(0,1,1)'\} \end{equation*}\]
The symmetric support with replacement defined as \[\begin{equation*} \mathcal{R}=\mathbb{N}^N \end{equation*}\] where \(\mathbb{N}\) is the set of natural numbers. Note that this support is a countable but infinite set; therefore \[\begin{equation*} \#(\mathcal{R})=\infty \end{equation*}\]
The fixed-size symmetric support with replacement defined as \[\begin{equation*} \mathcal{R}_m=\left\{ \textbf{s} \in \mathcal{R} | \sum_{k\in U}n_k=m \right\} \end{equation*}\] Note that \[\begin{equation*} \#(\mathcal{R}_m)=\binom{N+m-1}{m} \end{equation*}\] For example, if \(N=3\) and \(m=2\), then \(\mathcal{R}_m\) is defined by the following samples: \[\begin{equation*} \mathcal{R}_m=\{(2,0,0)',(0,0,2)',(0,2,0)',(1,1,0)',(1,0,1)',(0,1,1)'\} \end{equation*}\]
Till’e (2006) states that, geometrically, each vector \(\mathbf{s}\) represents the vertex of an \(N\)-cube. In addition, we have the following result:
TipResult
For the supports defined above, the following properties hold:
- \(\mathcal{S},\mathcal{S}_n,\mathcal{R},\mathcal{R}_m\) are symmetric supports.
- \(\mathcal{S}\subset\mathcal{R}\).
- The set \(\{\mathcal{S}_0,\mathcal{S}_1,\ldots,\mathcal{S}_N\}\) is a partition of \(\mathcal{S}\).
- The set \(\{\mathcal{R}_0,\mathcal{R}_1,\ldots,\mathcal{S}_N,\ldots\}\) is an infinite partition of \(\mathcal{R}\).
- \(\mathcal{S}\subset\mathcal{R}\) for all \(n=0,1,\ldots,N\).
2.1.3.1 Probability Samples
Not all random samples are probability samples. A sample (with or without replacement) is a probability sample if:
- It is possible to construct (or at least define theoretically) a support \(Q\), such that \(Q=\{s_1,\ldots,s_q,\ldots,s_Q\}\), of all possible samples obtained by a selection method. Here \(s_q\), \(q=1,\ldots,Q\), is a sample belonging to the support \(Q\).
- The selection probabilities assigned by the random process to each possible sample belonging to the support are known before the final sample is selected.
Note that a random sample is not necessarily a probability sample. In poor practice, some researchers use random methods for including elements without having a sampling frame and without satisfying the two preceding conditions; in this way, although the elements are chosen randomly or by chance, the resulting sample cannot be classified as a probability sample. From this point forward, unless otherwise stated, the term sample refers to a probability sample. Some comments of interest are:
- The universe \(U\) is finite.
- The probability sample \(s\) may contain repeated objects. This occurs when the sampling procedure is with replacement.
- The sample \(s\) with repetitions may have a size greater than the population.
- The sample \(s\) without repetitions may have a maximum size equal to \(N\).
- If there is no sampling frame, it is impossible to carry out a probability sampling procedure, except when a census is conducted.
- If the selected sample is not probabilistic, then no statistical estimate can be constructed.
- The statistician must be accountable for deceptions or frauds committed, whether through ignorance, bad faith, or the convenience of keeping a job or business for which they are not qualified, against clients, cities, and countries that trust the figures resulting from their analyses.
NoteExample
Suppose a finite population of size \(N=5\), in which the members of the population are each identified by name. The population is composed of the following elements:
In R, a character vector is used to index the population. Note that the elements belonging to the vector are specified by using quotation marks. In this case, the identifiers of each population element are assigned to the object U.
To obtain the support \(Q\) of all possible samples of size \(n=2\) from this population of size \(N=5\), the Support function from the TeachingSampling package is used. This function has three arguments: the population size N, the fixed size of each possible sample n, and finally a characteristic y, which may be numeric or may be a set of labels. The output of the function will be a data set containing all possible fixed-size samples. When the argument y is different from FALSE, the result of the function will be the population characteristic for each individual. In the following example, the function Support(N,n,y=FALSE) is used to obtain the set of possible samples of size two from the population \(U\), while the function Support(N,n,U) returns the set of labels in each of the 10 possible samples.
N <- length(U)
N[1] 5n <- 2
Support(N, n) [,1] [,2]
[1,] 1 2
[2,] 1 3
[3,] 1 4
[4,] 1 5
[5,] 2 3
[6,] 2 4
[7,] 2 5
[8,] 3 4
[9,] 3 5
[10,] 4 5Support(N, n, U) [,1] [,2]
[1,] "Yves" "Ken"
[2,] "Yves" "Erik"
[3,] "Yves" "Sharon"
[4,] "Yves" "Leslie"
[5,] "Ken" "Erik"
[6,] "Ken" "Sharon"
[7,] "Ken" "Leslie"
[8,] "Erik" "Sharon"
[9,] "Erik" "Leslie"
[10,] "Sharon" "Leslie"
ImportantDefinition
A sampling design \(p(\cdot)\) is a multivariate probability distribution defined on a support \(Q\); that is, \(p(\cdot)\) is a function that maps from \(Q\) to \((0,1]\) such that \(p(s)>0\) for all \(s \in Q\) and \[\begin{equation} \sum_{s \in Q}p(s)=1 \end{equation}\]
Given the support \(Q\), a sampling design is a function \(p(\cdot)\) such that \(p(s)\) gives the selection probability of the realized sample \(s\) under a particular selection scheme. In other words, if \(S\) is a random sample that takes the value \(s\) with probability \(p(s)\), such that
\[\begin{equation} Pr(S=s)=p(s) \ \ \ \ \ \ \ \ \ \ \ \ \text{for all } s\in Q . \end{equation}\]
Then \(p(\cdot)\) is called the sampling design.
The sampling design is a function that maps from the support \(Q\) to the interval \(]0,1]\). Because it is a probability distribution, \(p(\cdot)\) satisfies:
- \(p(s)\geq0\) for all \(s\in Q\)
- \(\sum_{s\in Q}p(s)=1\)
Note that the sampling design does not refer to an algorithm or procedure that enables sample selection. Given a sampling design, the statistician’s job is to find an algorithm that allows the selection of samples whose selection probability corresponds to the probability induced by the sampling design. For making inferences about the parameters of interest, the sampling design plays a very important role because the statistical properties (expectation, variance, and others) of the random quantities calculated from a sample are determined by it.
Given a support \(Q\), a sampling design may be:
- Without replacement if all possible samples in \(Q\) are without replacement.
- With replacement if all possible samples in \(Q\) are with replacement.
- Fixed-size if all possible samples in \(Q\) have the same sample size \(n(S)=n\).
Cassel et al. (1976) explain that the possibility of identifying each and every possible sample belonging to the support \(Q\) is a crucial factor that makes it possible to:
- designate a set of samples to which a positive selection probability is assigned, and
- distribute the entire probability mass among the members of \(Q\).
The most important feature of probability sampling is that it makes it possible to know, at least theoretically, the selection probability of all possible samples in the support \(Q\). However, a sampling design also makes it possible to know the inclusion probability of element \(k\) in the sample \(S\).
2.1.3.2 Selection Algorithm
A sampling design is a probability distribution over a support \(Q\); however, it is in no way a procedure that selects the sample per se.
ImportantDefinition
A selection algorithm is a procedure used to select a probability sample.
Till’e (2006) states that one way to select a sample is to list all possible samples, generate a random variable with a uniform distribution on the interval \([0,1]\), and then make the corresponding selection. Algorithms of this type, which list all possible samples, are known as enumerative selection algorithms; however, these algorithms are computationally inefficient and can only be implemented when the sampling design is known and the population size \(N\) is small. Throughout the book, various selection algorithms specific to each sampling design will be included, allowing the selection of a probability sample.
2.1.4 Inclusion Probability
The inclusion of the \(k\)-th element in a particular sample \(s\) is a random event defined by the indicator function \(I_k(s)\), which is given by
\[ I_k(s)= \begin{cases} 1 & \text{if $k \in s$}\\ 0 & \text{if $k \notin s$}. \end{cases} \tag{2.1}\]
Note that the function \(I_k(s)\) is a function of the random variable \(S\). To shorten the notation, we will write \(I_k=I_k(s)\), understanding that \(I_k\) is the indicator function for the \(k\)-th element. Under a sampling design \(p(\cdot)\), an inclusion probability is assigned to each element of the population to indicate the probability that the element belongs to the sample. For the \(k\)-th element of the population, the inclusion probability is denoted by \(\pi_k\) and is known as the first-order inclusion probability, given by
\[\begin{equation} \pi_k=Pr(k \in S)=Pr(I_k=1)=\sum_{s \ni k} p(s). \end{equation}\]
Here the subscript \(s \ni k\) refers to the sum over all samples that contain the \(k\)-th element. Note that, from the preceding definition, for a sample to be considered probabilistic, all elements in the population must have inclusion probability strictly greater than zero.
ImportantDefinition
The expectation of a random sample, in the sense of Definitions 2.1.2 and 2.1.3, is given by \[\begin{equation} \boldsymbol{\mu}=E(\mathbf{s})=\sum_{\mathbf{s}\in Q}p(\mathbf{s})\mathbf{s} \end{equation}\]
If the sampling design is without replacement, then \(\boldsymbol{\mu}=\boldsymbol{\pi}\), where \(\boldsymbol{\pi}=(\pi_1,\ldots,\pi_N)'\) is the vector of inclusion probabilities induced by the sampling design. The following result provides a simple way to compute the \(N\) inclusion probabilities.
TipResult
Given a support \(Q\), the inclusion probability \(\pi_k\) is the probability that the \(k\)-th element belongs to the random sample \(S\) and can be written as follows: \[\begin{equation} \pi_k=E(I_k(S))=\sum_{s \in Q}I_k(s)p(s) \end{equation}\]
Proof.
\(I_k(S)\) is a function of the random sample \(S\); the proof follows from the definition of the expectation of a function of a random variable. On the other hand, \(I_k(S)\) can take only the two values 1 and 0; therefore \[\begin{align*} E(I_k(S))&=(1)Pr(I_k(S)=1)+(0)Pr(I_k(S)=0)\\ &=Pr(I_k(S)=1)=Pr(k\in S)=\pi_k \end{align*}\]
Analogously, \(\pi_{kl}\) is known as the second-order inclusion probability and denotes the probability that elements \(k\) and \(l\) belong to the sample. It is denoted by \(\pi_{kl}\) and is given by
\[\begin{equation} \pi_{kl}=Pr(k\in S\text{ and }l\in S)=Pr(I_kI_l=1)=\sum_{s \ni \text{ $k$ and $l$}} p(s). \end{equation}\]
Here the subscript \(s \ni \text{ $k$ and $l$}\) refers to the sum over all samples that contain the \(k\)-th and \(l\)-th elements.
NoteExample
Consider the following sampling design \(p(\cdot)\), which assigns the following selection probabilities to each of the 10 possible samples of size 2 from the support \(Q\) of the population \(U\).
p <- c(0.13, 0.2, 0.15, 0.1, 0.15, 0.04, 0.02, 0.06, 0.07, 0.08)
p [1] 0.13 0.20 0.15 0.10 0.15 0.04 0.02 0.06 0.07 0.08That is, the first sample has a selection probability of 0.13, the second sample has a selection probability of 0.15, and so on until the tenth, whose selection probability is 0.08. With the following instructions, we verify that the properties of the sampling design are satisfied.
sum(p)[1] 1p < 0 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEBy using the Ik function from the TeachingSampling package, it is possible to create the \(N=5\) indicator functions of the population elements for each of the 10 possible fixed-size samples without replacement. This function has two arguments: the population size N and the fixed size of each possible sample n. A data table is created from the labels, the selection probability, and the 5 indicator functions of the possible samples contained in the support \(Q\).
Ind <- Ik(N, n)
Q <- Support(N, n, U)
data.frame(Q, p, Ind) X1 X2 p X1.1 X2.1 X3 X4 X5
1 Yves Ken 0.13 1 1 0 0 0
2 Yves Erik 0.20 1 0 1 0 0
3 Yves Sharon 0.15 1 0 0 1 0
4 Yves Leslie 0.10 1 0 0 0 1
5 Ken Erik 0.15 0 1 1 0 0
6 Ken Sharon 0.04 0 1 0 1 0
7 Ken Leslie 0.02 0 1 0 0 1
8 Erik Sharon 0.06 0 0 1 1 0
9 Erik Leslie 0.07 0 0 1 0 1
10 Sharon Leslie 0.08 0 0 0 1 1Once the indicator variables have been calculated for each element and in each possible sample, calculating the inclusion probabilities is very simple: multiply the selection probabilities by each of the indicator variables. The result is summed by columns, and the output is a vector of size \(N=5\) of inclusion probabilities.
multip <- p * Ind
colSums(multip)[1] 0.58 0.34 0.48 0.33 0.27The Pik function from the TeachingSampling package returns the vector of inclusion probabilities for all population elements. It has two arguments: a vector p of selection probabilities for all possible samples and a matrix Ind of \(N\) indicator variables. Note that the sum of inclusion probabilities is the expected sample size, in this case equal to 2.
pik <- Pik(p, Ind)
pik [,1] [,2] [,3] [,4] [,5]
[1,] 0.58 0.34 0.48 0.33 0.27Thus, the population element with the highest probability of being included is Yves, while the element with the lowest inclusion probability is Sharon. On the other hand, by using the Pikl function from the TeachingSampling package, it is possible to calculate the matrix of second-order inclusion probabilities for the design p under consideration. This function has only three arguments: N, the population size; n, the fixed sample size; and p, the sampling design used. The output of this function is a square and symmetric matrix of size \(N \times N\) whose entries correspond to second-order inclusion probabilities. In this particular case, the function is run as follows.
pikl <- Pikl(N, n, p)
pikl [,1] [,2] [,3] [,4] [,5]
[1,] 0.58 0.13 0.20 0.15 0.10
[2,] 0.13 0.34 0.15 0.04 0.02
[3,] 0.20 0.15 0.48 0.06 0.07
[4,] 0.15 0.04 0.06 0.33 0.08
[5,] 0.10 0.02 0.07 0.08 0.27Note that, under this sampling design, Yves and Erik correspond to the pair of elements with the highest inclusion probability.
2.1.5 Characteristic of Interest and Parameters of Interest
The purpose of any sampling study is to study a characteristic of interest \(y\) associated with each unit of the population. That is, the characteristic of interest takes the value \(y_k\) for unit \(k\). It is important to note that the \(y_k\)s are not considered random variables but fixed quantities; therefore, their notation uses the lowercase letter \(y\). The objective of sampling research is to estimate a function of interest \(T\), called a parameter, of the characteristic of interest in the population.
\[\begin{equation*} T=f\{y_1,\ldots,y_k,\ldots,y_N\}. \end{equation*}\]
Some of the most common parameters of interest are:
- The population total, \[ t_y=\sum_{k \in U}y_k \tag{2.2}\]
- The population mean, \[ \bar{y}_U=\frac{\sum_{k \in U}y_k}{N}=\frac{t_y}{N} \tag{2.3}\]
- The population variance, \[\begin{equation} S^2_{yU}=\frac{\sum_{k \in U}(y_k-\bar{y}_U)^2}{N-1} \end{equation}\]
There are other parameters of interest, such as the population median, population percentiles, the ratio between two population totals, or, as mentioned earlier, the size of a population, in which case we would be interested in \(N\). Among others, some examples of sampling studies interested in the preceding parameters are:
- Total number of people belonging to the labor force.
- Percentage of people who would use a product.
Obviously, these population quantities are unknown, and this is why a sampling study is required: through it, these population parameters can be estimated from a selected sample.
NoteExample
Suppose that, in our example population, we want to estimate the total of the variable \(y\). The value for each population element is as follows:
y <- c(32, 34, 46, 89, 35)
y[1] 32 34 46 89 35The data.frame function creates the data set containing the names (labels) and the value of the characteristic of interest for each population element.
data.frame(U, y) U y
1 Yves 32
2 Ken 34
3 Erik 46
4 Sharon 89
5 Leslie 35Some population parameters of interest for the characteristic \(y\) are the population total and mean, given by \(t_y\) and \(\bar{y}_U\), respectively.
ty <- sum(y)
ty[1] 236ybar <- ty / N
ybar[1] 472.1.6 Statistic and Estimator
A statistic is a real-valued function \(G\) of the random sample \(S\) and depends only on the elements belonging to \(S\). When a statistic is used to estimate a parameter, it is called an estimator, and the realizations of the estimator in a selected sample \(s\) are called estimates.
Since \(G\) is a statistic, its statistical properties are determined by the sampling design. That is, given the selection probability of each sample \(s \in Q\), the expectation, variance, and other properties of interest are defined from \(p(s)\).
The expectation of a statistic \(G\) is
\[\begin{equation} E(G)=\sum_{s\in Q}p(s)G(s). \end{equation}\]
The variance of the statistic \(G\) is defined as
\[\begin{align} Var(G)&=E[G-E(G)]^2\\ &=\sum_{s\in Q}p(s)[G(s)-E(G)]^2. \end{align}\]
Where \(G(s)\) is the real value taken by the statistic \(G\) in the selected (realized) sample \(s\), and \(Q\) is the support induced by the sampling design. Note that the properties of statistics and, consequently, of estimators are defined with sums because the sampling design induces a discrete probability distribution over all possible samples \(s\) belonging to the support \(Q\).
2.1.6.1 The Statistic \(I_k\)
The quantity \(I_k\) given by equation 2.1 is a statistic that takes values randomly depending on the sampling design used.
TipResult
The most important properties of this statistic are:
- \(E(I_k)=\pi_k\)
- \(Var(I_k)=\pi_k(1-\pi_k)\)
- \(Cov(I_k,I_l)=\pi_{kl}-\pi_k\pi_l\) for all \(k \neq l\)
Proof.
By Result 2.1.2, the first property follows immediately; now, from the definition of variance, we have \[\begin{align*} Var(I_k(S))&=E[I_k(S)-E(I_k(S))]^2\\ &=Pr(I_k(S)=1)[1-\pi_k]^2+Pr(I_k(S)=0)[0-\pi_k]^2\\ &=\pi_k(1-\pi_k) \end{align*}\] and finally, from the definition of covariance we have \[\begin{align*} Cov(I_k(S),I_l(S))&=E[I_k(S)I_l(S)]-E[I_k(S)]E[I_k(S)]\\ &=(1)Pr(I_k(S)I_l(S)=1)+(0)Pr(I_k(S)I_l(S)=0)-\pi_k\pi_l\\ &=\pi_{kl}-\pi_k\pi_l \end{align*}\]
The covariance of the indicator statistics for elements \(k\) and \(l\), \(Cov(I_k,I_l)\), is known as \(\Delta_{kl}\). This quantity, depending on the design, may take positive, negative, or even zero values.
2.1.6.2 The Statistic \(n(S)\) or Sample Size
As already seen, the sample size is a random quantity, depending on the design. Note that this value can be expressed as a function of the inclusion statistics.
\[\begin{align} n(S)=\sum_UI_k. \end{align}\]
TipResult
Some properties of interest are:
- \(E(n(S))=\sum_U\pi_k\)
- \(Var(n(S))=\sum_U\pi_k-(\sum_U\pi_k)^2+\sum\sum_{k \neq l}\pi_{kl}\).
Proof.
For the first property, we have \[\begin{align*} E[n(S)]=E\left[\sum_UI_k\right]=\sum_UE[I_k]=\sum_U\pi_k \end{align*}\] Recalling the properties of the variance of a sum, we have \[\begin{align*} Var[n(S)]&=Var\left[\sum_UI_k\right]\\ &=\sum_UVar[I_k]+\sum\sum_{k\neq l}Cov[I_k,I_l]\\ &=\sum_U\pi_k-\sum_U\pi_k^2-\sum\sum_{k \neq l}\pi_k\pi_l+\sum\sum_{k \neq l}\pi_{kl}\\ &=\sum_U\pi_k-\left(\sum_U\pi_k\right)^2+\sum\sum_{k \neq l}\pi_{kl} \end{align*}\]
In addition, when the variation of the sample size is zero because a fixed sample-size design has been used, the following properties hold.
TipResult
If the sampling design is fixed-size and equal to \(n\),
- \(E(n(S))=\sum_U\pi_k=n\)
- \(\sum_U\pi_{kl}=n\pi_l\)
- \(\sum_U\Delta_{kl}=0\)
- \(\pi_k(1-\pi_k)=\sum_{l\neq k}(\pi_k\pi_l-\pi_{kl})\)
Proof.
The first property follows by recalling that the expectation of a constant is the constant itself. Note that \(\pi_{kl}= E[I_k(S)I_l(S)]\), so \[\begin{align*} \sum_{l\in U}\pi_{kl}=\sum_{l\in U}E[I_k(S)I_l(S)]&=\sum_{l\in U}\sum_{s\in Q}p(s)I_k(s)I_l(s)\\ &=\sum_{s\in Q}p(s)I_k(s)\sum_{l\in U}I_l(s)\\ &=n(S)\sum_{s\in Q}p(s)I_k(s)=n\pi_k \end{align*}\] The third property holds because \[\begin{align*} \sum_U\Delta_{kl}&=\sum_U(\pi_{kl}-\pi_k\pi_l)\\ &=\sum_U\pi_{kl}-\pi_k\sum_U\pi_l\\ &=n\pi_k-n\pi_k=0 \end{align*}\] To prove the last property, it is necessary to redefine the sample size so that \(n=\sum_{l\neq k}I_l(S)+I_k(S)\). Then, \[\begin{align*} \pi_k(1-\pi_k)&=Var(I_k(S))\\ &=Cov(I_k(S),I_k(S))\\ &=Cov\left(I_k(S),n-\sum_{l\neq k}I_l(S)\right)\\ &=-\sum_{l\neq k}Cov(I_k(S),I_l(S))\\ &=\sum_{l\neq k}(\pi_k\pi_l-\pi_{kl}) \end{align*}\]
NoteExample
Continuing with the development of Example 2.1.3, we will now use the vector of inclusion probabilities and the matrix of second-order probabilities to verify Results 2.1.4 and 2.1.5. First, note that the expectation of the sample size, which corresponds to 2 because the design is fixed-size, is obtained as follows.
A <- sum(pik)
A[1] 2Now, the square of the sum of the inclusion probabilities is obtained as follows.
B <- (sum(pik))^2
B[1] 4And the sum of the distinct elements of the matrix of second-order inclusion probabilities is
C <- sum(pikl) - sum(diag(pikl))
C[1] 2To verify the second part of Result 2.1.4, it is enough to carry out the operation A-B+C. This sum is zero and indeed corresponds to the variance of the sample size in this sampling design; since, in this particular case, the sample size was always fixed and equal to 2, the variance must be zero.
The next step in this example consists of verifying the second part of Result 2.1.5. In short, this section says that the row (or column) sums of the matrix of second-order inclusion probabilities must correspond exactly to the product of the sample size and the vector of first-order inclusion probabilities. This is easily corroborated with the following code.
n * pik [,1] [,2] [,3] [,4] [,5]
[1,] 1.2 0.68 0.96 0.66 0.54colSums(pikl)[1] 1.16 0.68 0.96 0.66 0.54rowSums(pikl)[1] 1.16 0.68 0.96 0.66 0.54Note that the row and column sums coincide perfectly with \(n\times \pi_k\) for all \(k\in U\). On the other hand, we will verify the third property, which states that the row (or column) sums of the variance-covariance matrix of the sample membership indicator variables must yield a vector of zeros of size five. For this, the Deltakl function from the TeachingSampling package is used. This function has three arguments: N, the population size; n, the fixed sample size; and p, the sampling design used. The output of this function is a square and symmetric matrix of size \(N \times N\) whose entries correspond to the variances and covariances of the sample membership indicator variables. For this example, implementing the following code makes it possible to obtain the desired matrix and verify the result.
Delta <- Deltakl(N, n, p)
Delta [,1] [,2] [,3] [,4] [,5]
[1,] 0.244 -0.067 -0.078 -0.0414 -0.0566
[2,] -0.067 0.224 -0.013 -0.0722 -0.0718
[3,] -0.078 -0.013 0.250 -0.0984 -0.0596
[4,] -0.041 -0.072 -0.098 0.2211 -0.0091
[5,] -0.057 -0.072 -0.060 -0.0091 0.1971rowSums(Delta)[1] -0.000000000000000139 -0.000000000000000083 -0.000000000000000056
[4] -0.000000000000000069 -0.000000000000000014colSums(Delta)[1] -0.000000000000000139 -0.000000000000000083 -0.000000000000000056
[4] -0.000000000000000069 -0.000000000000000014In this way, the row (or column) sum of the variance-covariance matrix of the sample membership indicator variables is zero in each column (or row).
When a statistic is constructed with the intention of estimating a parameter, it is called an estimator. Thus, the most commonly used properties of an estimator \(\hat{T}\) of a parameter of interest \(T\) are the bias, defined by
\[\begin{equation} B(\hat{T})=E(\hat{T})-T \end{equation}\]
and the mean squared error, given by
\[\begin{align} MSE(\hat{T})&=E[\hat{T}-T]^2\\ &=Var(\hat{T})+B^2(\hat{T}). \end{align}\]
If the bias of an estimator is zero, the estimator is said to be unbiased, and when this occurs the mean squared error becomes the variance of the estimator.
Särndal et al. (1992) state that the objective in a sampling study is to estimate one or more population parameters. The most important decisions when addressing a sampling estimation problem are
- The choice of a sampling design and a selection algorithm that allows the design to be implemented.
- The choice of a mathematical formula or estimator that calculates an estimate of the parameter of interest in the selected sample.
The preceding decisions are not independent. That is, the choice of an estimator will usually depend on the sampling design used.
ImportantDefinition
Let \(\hat{T}\) be an estimator of a parameter \(T\) and let \(p(\cdot)\) be a sampling design defined on a support \(Q\). A sampling strategy is defined as the pair \((p(\cdot),\hat{T})\).
As its name indicates, this book focuses on the search for the best combination of sampling design and estimator; this problem has been considered throughout the development of sampling theory. The sampling strategy is chosen in two stages, namely: the design stage, referring to the period during which the sampling design to be used is decided, together with the sampling algorithm that enables sample selection, and finally the probability sample is selected. Once the information is collected and recorded, the estimation stage begins, in which estimates for the characteristic of interest are calculated using the estimator associated with the chosen sampling strategy.
2.2 Sampling Estimators
Each element belonging to the population has an associated characteristic of interest \(y\). For the \(k\)-th element, the value taken by this characteristic of interest is \(y_k\). The objective of sampling research is to estimate a parameter \(T\) that is of interest. The statistician’s objective is to infer about \(T\) based on a sample \(s\). One indicator of the precision of an estimator is the estimated coefficient of variation, given by
\[\begin{equation} cve(\hat{T})=\frac{\sqrt{\widehat{Var}(\hat{T})}}{\hat{T}} \end{equation}\]
where \(\widehat{Var}(\hat{T})\) is the variance estimator based on the selected sample \(s\). The estimated coefficient of variation is a commonly used measure for expressing the error incurred by selecting a sample and not using the entire population in measuring the variable of interest. If a census were conducted and the estimator reproduced the population parameter, then \(\widehat{Var}(\hat{T})\) would be zero and, therefore, the \(cve\) would also be zero.
Next, some of the most widely used estimators in the history of sampling are reviewed. As the reading of the book progresses, new estimators will emerge and, consequently, new sampling strategies that make it possible to obtain results with almost clinical precision. Most of the estimators presented in this book are estimators of totals or of functions of totals.
2.2.1 The Horvitz-Thompson Estimator
2.2.1.1 Estimator of the Population Total
Narain (1951) discovered this estimator, although his article was edited and published by an Indian journal with limited circulation. Later, Horvitz and Thompson (1952) published similar results in the most important statistics journal at that time, JASA (Journal of the American Statistical Society). Since then, this estimator has been known as the Horvitz-Thompson estimator or \(\pi\) estimator, although rigorously it should be called the Narain-Horvitz-Thompson estimator. In this book we will follow the international and classical notation.
For a universe \(U\), we want to estimate the population total \(t_y\) of the characteristic of interest \(y\) given by equation 2.2. The Horvitz-Thompson (HT) estimator for \(t_y\) is defined as:
\[ \hat{t}_{y,\pi}=\sum_S\frac{y_k}{\pi_k}=\sum_Sd_ky_k \tag{2.4}\]
Where \(\pi_k\) is the inclusion probability for the \(k\)-th element, and \(d_k\) is known as the expansion factor and corresponds to the inverse of the inclusion probability. Note that the Horvitz-Thompson estimator is random because it is constructed from a sum over the random sample \(S\). The motivation behind this estimator, as Brewer (2002) indicates, rests on the principle of representativeness, which states that each element included in a sample represents itself and a group of units that do not belong to the selected sample, whose characteristics are close to those of the element included in the sample. The expansion factor is nothing other than the number of population elements minus one (not included in the sample) represented by the included element.
TipResult
If all first-order inclusion probabilities are greater than zero (\(\pi_{k}>0\) for all \(k\)), the Horvitz-Thompson estimator is unbiased for the population total. Therefore, \[\begin{equation} E(\hat{t}_{y,\pi})=t_y \end{equation}\]
Proof.
Rewriting the Horvitz-Thompson estimator as \(\hat{t}_{y,\pi}=\sum_SI_k(S)\frac{y_k}{\pi_k}\), we have \[\begin{align*} E(\hat{t}_{y,\pi})=E\left(\sum_UI_k(S)\frac{y_k}{\pi_k}\right) =\sum_U\frac{y_k}{\pi_k}E\left(I_k(S)\right)=\sum_U\pi_k\frac{y_k}{\pi_k}=t_y \end{align*}\]
If the sampling design is such that the first-order inclusion probabilities maintain a good positive correlation with the measurement of the characteristic of interest; in other words, if \(\pi_k \propto y_k\), the Horvitz-Thompson estimator reduces to a constant and therefore has zero variance. In practice, an optimal sampling strategy (Cassel et al. 1976) is one that uses the Horvitz-Thompson estimator together with a sampling design that induces a good correlation between the vector of inclusion probabilities and the vector of values of the characteristic of interest. However, in multipurpose surveys, where the goal is to estimate parameters for several characteristics of interest among which there is no good correlation, when using the Horvitz-Thompson estimator it is difficult to avoid the weak, and even negative, correlation that exists between the characteristics of interest and the vector of inclusion probabilities. Nevertheless, this fact can be mitigated by including auxiliary information in the construction of the estimator.
2.2.1.2 Variance of the Horvitz-Thompson Estimator
TipResult
The variance of the Horvitz-Thompson estimator is given by the following expression \[ Var_1(\hat{t}_{y,\pi})=\sum\sum_U\Delta_{kl}\frac{y_k}{\pi_k}\frac{y_l}{\pi_l}. \tag{2.5}\]
Proof.
From the definition of variance, we obtain the following \[\begin{align*} Var_1(\hat{t}_{y,\pi})&=Var\left(\sum_UI_k(S)\frac{y_k}{\pi_k}\right)\\ &=\sum_U\frac{y_k^2}{\pi_k^2}Var(I_k(S))+\sum\sum_{k\neq l}\frac{y_k}{\pi_k}\frac{y_l}{\pi_l}Cov(I_k(S),I_l(S))\\ &=\sum_U\frac{y_k^2}{\pi_k^2}(\pi_k-\pi_k^2) +\sum\sum_{k\neq l}\frac{y_k}{\pi_k}\frac{y_l}{\pi_l}(\pi_{kl}-\pi_k\pi_l)\\ &=\sum\sum_U\frac{y_k}{\pi_k}\frac{y_l}{\pi_l}(\pi_{kl}-\pi_k\pi_l)\\ &=\sum\sum_U\Delta_{kl}\frac{y_k}{\pi_k}\frac{y_l}{\pi_l} \end{align*}\]
Sen (1953) and Yates and Grundy (1953) derived the following result when the sampling design is fixed-size.
TipResult
If the design \(p(\cdot)\) has fixed sample size, then the variance of the Horvitz-Thompson estimator is written as \[ Var_2(\hat{t}_{y,\pi})=-\frac{1}{2}\sum\sum_U\Delta_{kl}\left(\frac{y_k}{\pi_k}-\frac{y_l}{\pi_l}\right)^2 \tag{2.6}\]
Proof.
Using the properties of Result 2.1.5, we have \[\begin{align*} Var_2(\hat{t}_{y,\pi})&= -\frac{1}{2}\sum\sum_{U}\Delta_{kl}\left(\frac{y_k}{\pi_k}-\frac{y_l}{\pi_l}\right)^2\\ &=-\frac{1}{2}\sum\sum_{U}\Delta_{kl}\left(\frac{y_k^2}{\pi_k^2}+\frac{y_l^2}{\pi_l^2}- 2\frac{y_k}{\pi_k}\frac{y_l}{\pi_l}\right)\\ &=-\frac{1}{2}\left[\sum\sum_{U}\Delta_{kl}\frac{y_k^2}{\pi_k^2}+\sum\sum_{U}\Delta_{kl} \frac{y_l^2}{\pi_l^2}-2\sum\sum_U\Delta_{kl}\frac{y_l}{\pi_k}\frac{y_k}{\pi_l}\right]\\ &=-\frac{1}{2}\left[2\sum\sum_{U}\Delta_{kl}\frac{y_k^2}{\pi_k^2}- 2\sum\sum_U\Delta_{kl}\frac{y_l}{\pi_k}\frac{y_k}{\pi_l}\right]\\ &=-\sum_{U}\frac{y_k^2}{\pi_k^2}\sum_U\Delta_{kl}+\sum\sum_U\Delta_{kl}\frac{y_l}{\pi_k}\frac{y_k}{\pi_l}\\ &=\sum\sum_U\Delta_{kl}\frac{y_l}{\pi_k}\frac{y_k}{\pi_l} =Var_1(\hat{t}_{y,\pi}) \end{align*}\] since \(\sum_U\Delta_{kl}=0\) for fixed-size designs. Therefore, in the case of fixed-size sampling designs, the variance of the Horvitz-Thompson estimator can be calculated by means of \(Var_2(\hat{t}_{y,\pi})\).
2.2.1.3 Variance Estimation
It is possible to construct two unbiased estimators for equation 2.5 and equation 2.6. For this, all second-order inclusion probabilities must be strictly positive (\(\pi_{kl}>0\) for all \(k\)). Under the preceding assumption, the following results hold.
TipResult
An unbiased estimator for equation 2.5 is given by \[ \widehat{Var}_1(\hat{t}_{y,\pi})=\sum\sum_S \dfrac{\Delta_{kl}}{\pi_{kl}}\frac{y_k}{\pi_k}\frac{y_l}{\pi_l} \tag{2.7}\]
TipResult
If the design has fixed sample size, an unbiased estimator for equation 2.6 is given by \[\begin{equation} \widehat{Var}_2(\hat{t}_{y,\pi})=-\frac{1}{2}\sum\sum_S\frac{\Delta_{kl}}{\pi_{kl}}\left(\frac{y_k}{\pi_k}-\frac{y_l}{\pi_l}\right)^2 \end{equation}\]
Proof.
The preceding results are immediate after rewriting the estimators \(\widehat{Var}_1(\hat{t}_{y,\pi})\) and \(\widehat{Var}_2(\hat{t}_{y,\pi})\) in terms of \(U\) and multiplying by the product of the indicator functions \(I_k(S)I_l(S)\). Applying expectation gives \(E[I_k(S)I_l(S)]=\pi_{kl}\), which proves the result.
Bautista (1998) highlights the following three important comments about the estimates produced by the preceding expressions.
- If the second-order inclusion probabilities are greater than zero for all elements in the sample, but not for the remaining elements that were not included in the sample, the unbiasedness of the preceding expressions cannot be guaranteed.
- Variance estimates may yield negative results, which cannot be used or interpreted. To avoid this situation, it is necessary to ensure that the covariance between the inclusion statistics for each pair of elements in the population is negative (\(\Delta_{kl}<0\) \(\forall\) $k l $).
- The estimates produced by the preceding expressions do not necessarily coincide in all cases.
For his part, Till’e (2006) adds that, in practice, using the expressions for the variance estimators is very difficult to implement because the double sum makes the computational calculation process very long and inefficient. Therefore, for each sampling design used, expressions that can be simplified must be created, or in some cases approximations must be used.
2.2.1.4 Confidence Interval for the Horvitz-Thompson Estimator
H’ajek (1960) proves the asymptotic convergence of the Horvitz-Thompson estimator to a normal distribution. When the sample size is sufficiently large (depending on the behavior of the population, a few dozen individuals may be enough), a confidence interval of level \((1-\alpha)\) for the population total \(t_y\) can be constructed as follows:
\[ IC(1-\alpha)=\left[\hat{t}_{y,\pi}-z_{1-\alpha / 2}\sqrt{ Var(\hat{t}_{y,\pi})},\hat{t}_{y,\pi}+z_{1-\alpha / 2}\sqrt{Var(\hat{t}_{y,\pi})}\right] \tag{2.8}\]
where \(z_{1-\alpha / 2}\) refers to the quantile \((1-\alpha / 2)\) of a random variable with a standard normal distribution. Note that
\[\begin{equation*} 1-\alpha=\sum_{Q_0 \supset s}p(s), \end{equation*}\]
where \(Q_0\) is the set of all possible samples whose confidence interval contains the population total \(t_y\). In practice, the variance of the estimator is very rarely known; therefore, the estimated confidence interval of level \((1-\alpha)\) can be obtained from the data in the selected sample by replacing, in equation 2.8, the variance of the estimator with its corresponding estimate, yielding the following expression
\[\begin{equation} IC_s(1-\alpha)=\left[\hat{t}_{y,\pi}-z_{1-\alpha / 2}\sqrt{ \widehat{Var}(\hat{t}_{y,\pi})},\hat{t}_{y,\pi}+z_{1-\alpha / 2}\sqrt{ \widehat{Var}(\hat{t}_{y,\pi})}\right]. \end{equation}\]
When using a sampling strategy to estimate a parameter in finite populations, the properties of the strategy are studied in terms of:
- Reliability: defined as the sum of the probabilities of the samples whose confidence interval contains the parameter.
- Precision: defined as the length of the confidence interval.
Note that the preceding properties are functions of the confidence interval. To determine reliability, the parameter \(T\) (unknown) must be known; therefore, in practical terms, reliability cannot be calculated. To determine precision and reliability, it is necessary to know the design-based variance of the estimator used, say \(\hat{T}\); however, calculating this variance \(Var(\hat{T})\) almost always requires knowing the values \(y_k\) for all \(k=1,...,N\). Thus, precision cannot be calculated either. However, an estimator of \(Var(\hat{T})\) should be proposed (preferably unbiased) which, together with \(\hat{T}\), provides a bound for bias and precision.
2.2.1.5 Estimation of Other Parameters
Although equation 2.4 is an estimator of the population total of the characteristic of interest, it can be used to estimate other population quantities of interest. If the population size \(N\) is known, the population mean defined in equation 2.3 can be estimated with the Horvitz-Thompson estimator.
TipResult
The population mean is estimated unbiasedly using the following expression \[\begin{align} \hat{\bar{y}}_{\pi}&=\dfrac{1}{N}\left(\hat{t}_{y,\pi}\right)=\dfrac{1}{N}\sum_s\frac{y_k}{\pi_k} \end{align}\] The variance and estimated variance of the estimator of the population mean are given by \[\begin{equation} Var(\hat{\bar{y}}_{\pi})=\dfrac{1}{N^2}Var(\hat{t}_{y,\pi}) \end{equation}\] \[\begin{equation} \hat{Var}(\hat{\bar{y}}_{\pi})=\dfrac{1}{N^2}\hat{Var}(\hat{t}_{y,\pi}) \end{equation}\] respectively.
However, it is the rule rather than the exception that, in most cases in which the user faces an investigation whose objectives depend on carrying out a sampling study, the population size is unknown. In that case, we can use the Horvitz-Thompson estimator to estimate it, since \(N\) can be written as follows
\[\begin{equation} N=\sum_U 1, \end{equation}\]
taking the familiar form of a population total. Thus, we have the following result.
TipResult
The population size is estimated unbiasedly using the following expression \[\begin{equation} \hat{N}_{\pi}=\sum_S\frac{1}{\pi_k}. \end{equation}\]
When the population total of a characteristic of interest and the population size have been estimated using the Horvitz-Thompson estimator, an estimator for the population mean arises, given by
\[\begin{align} \widetilde{y}_S&=\dfrac{\hat{t}_{y,\pi}}{\hat{N}_{\pi}}\\ &=\sum_S\dfrac{y_k}{\pi_k} \ \Bigl/ \ \sum_S\dfrac{1}{\pi_k}. \end{align}\]
The preceding expression is a ratio, or a quotient between two population totals. The statistical properties of the preceding estimators will be treated later in the relevant sections of the book.
Till’e (2006) notes that, even when \(N\) is known, a poor property of the Horvitz-Thompson estimator for the population mean appears when it is used and the characteristic of interest is constant for all elements of the population (\(y_k=C\) \(\forall k \in U\)). Of course, under the preceding conditions it is clear that the population mean is equal to the constant (\(\bar{y}_U=C\)). However, the estimator \(\hat{\bar{y}}_{\pi}\) takes the following form
\[\begin{equation} \hat{\bar{y}}_{\pi}=\dfrac{1}{N}\sum_s\frac{y_k}{\pi_k}=\dfrac{1}{N}\sum_s\frac{C}{\pi_k}=\dfrac{C}{N}\sum_s\frac{1}{\pi_k}=C\frac{\hat{N}_{\pi}}{N}. \end{equation}\]
In this regard, Bautista (1998) states that in cases where the value of \(N\) is known, it is preferable to ignore it and use the estimator \(\widetilde{y}_S\), since its variation is smaller and, when \(y_k=C\) \(\forall k \in U\), it reproduces the population mean with zero variance because
\[\begin{equation*} \widetilde{y}_S=\dfrac{\hat{t}_{y,\pi}}{\hat{\bar{y}}_{\pi}}=\dfrac{C\hat{\bar{y}}_{\pi}}{\hat{\bar{y}}_{\pi}}=C. \end{equation*}\]
When the population size is known and, as will be seen later, for some sampling designs without replacement, a new alternative estimator of the population total can be created, inspired by the following argument: if \(\widetilde{y}_S\) estimates the population mean, then \(N\widetilde{y}_S\) will estimate the population total. Therefore, the alternative estimator is given by the following expression
\[\begin{equation} \hat{t}_{y,alt}=N\widetilde{y}_S=\hat{t}_{y,\pi}\dfrac{N}{\hat{N}_{\pi}} \end{equation}\]
which can be seen as a correction of the Horvitz-Thompson estimator through estimation of the population size. The variance and variance estimation will be topics of later chapters.
NoteExample
The HT function from the TeachingSampling package returns the estimate of the population total for one or more characteristics of interest. This function has two arguments: the vector of size \(n\) of inclusion probabilities pik and the set of values of the characteristic or characteristics of interest for the individuals belonging to the sample. y may be a vector in the case of a single characteristic of interest or a matrix in the case of several.
Thus, if the first sample (whose elements are Yves and Ken) had been selected, and given that the inclusion probabilities of these two elements are 0.58 and 0.34, respectively, and the values of the characteristic of interest are 32 and 34, respectively, the Horvitz-Thompson estimator would produce the following estimate:
y.s <- c(32, 34)
pik.s <- c(0.58, 0.34)
HT(y.s, pik.s) [,1]
[1,] 155Note that the population total for the variable of interest \(y\) is equal to 236. On the other hand, the calculation or estimation of the variance of the Horvitz-Thompson estimator is not implemented because the double sum makes the computational processes very long and slow. Therefore, if these values are needed, the process must be carried out manually. The variance estimation is performed by taking into account that \(\pi_{12}=0.13\). Thus,
\[\begin{align*} \frac{\Delta_{11}}{\pi_{11}}&=\frac{\pi_{11}-\pi_{1}\pi_{1}}{\pi_{11}}=\frac{0.58-0.58^2}{0.58}=0.42\\ \frac{\Delta_{12}}{\pi_{12}}&=\frac{\pi_{12}-\pi_{1}\pi_{2}}{\pi_{12}}=\frac{0.13-0.58*0.34}{0.13}=-0.52\\ \frac{\Delta_{21}}{\pi_{21}}&=\frac{\pi_{11}-\pi_{2}\pi_{1}}{\pi_{21}}=\frac{0.13-0.34*0.58}{0.13}=-0.52\\ \frac{\Delta_{22}}{\pi_{22}}&=\frac{\pi_{22}-\pi_{2}\pi_{2}}{\pi_{22}}=\frac{0.34-0.34^2}{0.34}=0.66 \end{align*}\]
Therefore, using equation 2.7, the variance estimator is
\[\begin{align*} \widehat{Var}(\hat{t}_{\pi})=\frac{\Delta_{11}}{\pi_{11}}\frac{y_1}{\pi_1}\frac{y_1}{\pi_1} +\frac{\Delta_{12}}{\pi_{12}}\frac{y_1}{\pi_1}\frac{y_2}{\pi_2} +\frac{\Delta_{21}}{\pi_{21}}\frac{y_2}{\pi_2}\frac{y_1}{\pi_1} +\frac{\Delta_{22}}{\pi_{22}}\frac{y_2}{\pi_2}\frac{y_2}{\pi_2} \end{align*}\]
and its respective estimate is
\[\begin{align*} 0.42\left(\frac{32}{0.58}\right)^2-2(0.52)\left(\frac{32}{0.58}\frac{34}{0.34}\right)+0.66\left(\frac{34}{0.34}\right)^2\cong2140 \end{align*}\]
The estimated coefficient of variation is
\[\begin{equation*} cve(\hat{t}_{\pi})=\frac{\sqrt{2140}}{155.1724}\cong0.3 \end{equation*}\]
And the estimated confidence interval with a confidence level of 95 percent for this estimate is the following:
\[\begin{align*} IC_s(0.95)&\cong \left[155-(1.96)\sqrt{2140},155+(1.96)\sqrt{2140}\right]\\ &\cong \left[64,246\right] \end{align*}\]
Continuing with the lexical-graphic exercise of estimating the population total \(t_y\) over all possible samples of size 10 from the population \(U\), we have table 2.1, which can be reproduced by running the following computational code.
all.pik <- Support(N, n, pik)
all.y <- Support(N, n, y)
all.HT <- rep(0, 10)
for (k in 1:10) {
all.HT[k] <- HT(all.y[k, ], all.pik[k, ])
}
all.HT [1] 155 151 325 185 196 370 230 366 225 399AllSamples <- data.frame(Q, p, all.pik, all.y, all.HT)| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|
| Yves | Ken | 0.13 | 0.58 | 0.34 | 32 | 34 | 155 |
| Yves | Erik | 0.20 | 0.58 | 0.48 | 32 | 46 | 151 |
| Yves | Sharon | 0.15 | 0.58 | 0.33 | 32 | 89 | 325 |
| Yves | Leslie | 0.10 | 0.58 | 0.27 | 32 | 35 | 185 |
| Ken | Erik | 0.15 | 0.34 | 0.48 | 34 | 46 | 196 |
| Ken | Sharon | 0.04 | 0.34 | 0.33 | 34 | 89 | 370 |
| Ken | Leslie | 0.02 | 0.34 | 0.27 | 34 | 35 | 230 |
| Erik | Sharon | 0.06 | 0.48 | 0.33 | 46 | 89 | 366 |
| Erik | Leslie | 0.07 | 0.48 | 0.27 | 46 | 35 | 225 |
| Sharon | Leslie | 0.08 | 0.33 | 0.27 | 89 | 35 | 399 |
The vector all.HT contains the Horvitz-Thompson estimates for each of the 10 possible samples; its expectation is calculated as
sum(p * all.HT)[1] 236Note that the expectation of the Horvitz-Thompson estimator exactly reproduces the population total. The variance is calculated as follows.
\[\begin{multline*} Var(\hat{t}_{\pi})=(0.13)(155.2-236)^2+(0.2)(151.0-236)^2+\cdots\\ +(0.08)(399.3-236)^2=7847.2 \end{multline*}\]
Using the VarHT function from the TeachignSampling package, it is possible to reproduce this same variance calculation. However, this function uses the theoretical variance expression \(Var_1(\hat{t}_{y,\pi})\) given by equation 2.5 for fixed-size sampling designs. It has four arguments: y, a vector containing the values of the characteristic of interest for each and every element of the population; N, the population size; n, the fixed sample size; and p, the sampling design used. The result of this function is the calculation of the theoretical variance value of the Horvitz-Thompson estimator for a particular sampling design and configuration of population values. Continuing with the sampling design given in Example 2.1.2 and the configuration of values of the characteristic of interest from Example 2.1.3, the variance calculation is exactly the same as that given by the lexical-graphic exercise.
VarHT(y, N, n, p)[1] 78472.2.2 The Hansen-Hurwitz Estimator
2.2.2.1 On Sampling With Replacement
Consider a finite population of \(N\) elements and a sampling design that allows the selection of a realized sample \(s\), with replacement, of size \(m\). As Lohr (2000) states, the most intuitive way to understand this type of with-replacement sampling design is to think of drawing \(m\) independent samples of size 1. An element is drawn from the population to be included in the sample with probability \(p_k\); however, that same element participates in the next random draw. This process is repeated \(m\) times; that is, there is a total of \(m\) random draws.
Under the preceding selection scheme, it is clear that an element can be selected into the sample more than once; therefore, although the size of the selected sample with replacement is \(m\), the effective sample size is not necessarily \(m\). Note that selecting an element that is repeated more than once provides no new information. This is why, in practice, sampling designs that allow the selection of samples without duplicates are preferred.
Särndal et al. (1992) state that the general framework of sampling with replacement has the following characteristics:
- Each element of the population is directly associated with a positive number \(p_k\) (\(k=1,\ldots,N\)), such that \[\begin{equation*} \sum_Up_k=1. \end{equation*}\] \(p_k\) is known as the selection probability of the \(k\)-th element. Note that these probabilities are not necessarily equal.
- To select the first element that will belong to the sample of size \(m\), a random draw is carried out so that \[\begin{equation*} Pr(\text{Select element }k)=p_k,\text{ $k \in U$}. \end{equation*}\]
- The selected element is replaced in the population and again becomes part of the next random draw with the same selection probability \(p_k\).
- The same set of probabilities is used to select the remaining elements. In total, \(m\) independent random draws are carried out.
Now, in sampling with replacement, the selection probability of an element is not the same as its inclusion probability. The following results hold.
ImportantDefinition
Under a with-replacement design, the random variable \(n_k(S)\) is defined as the number of times the \(k\)-th element is selected in the random sample \(S\).
TipResult
The random variable \(n_k(S)\) follows a binomial distribution such that \[\begin{equation*} E(n_k(S))=mp_k, \ \ \ \ \ \ \ \ \ Var(n_k(S))=mp_k(1-p_k) \end{equation*}\]
Proof.
Since each of the \(m\) draws induces independent statistical events, the selection of the \(k\)-th element in a particular draw follows a Bernoulli distribution with parameter \(p_k\). Since there are \(m\) draws, \(n_k(S)\) follows a binomial distribution and can take the values \(0,1,\ldots,m\); defining success as the selection of the \(k\)-th element in the sample proves the result.
ImportantDefinition
In general, a sampling design with replacement is defined as \[ p(s)= \begin{cases} \frac{m!}{n_1(s)!\ldots n_N(s)!}\prod_U(p_k)^{n_k(s)} &\text{if $\sum_Un_k(s)=m$}\\ 0 &\text{otherwise} \end{cases} \tag{2.9}\] Where \(n_k(s)\) is the number of times the \(k\)-th element is selected in the realized sample \(s\).
Note the difference (and at the same time the similarity) between the variable \(n_k(S)\) and the variable \(I_k(S)\). Also, by the preceding definition, the sampling design with replacement follows a multinomial distribution and therefore satisfies the conditions of a sampling design; that is, \(\sum_{s\in Q}p(s)=1\), where \(Q\) is the support containing all possible samples with replacement of size \(m\). The cardinality of \(Q\) is
\[\begin{equation} \#Q=\binom{N+m-1}{m} \end{equation}\]
TipResult
In sampling with replacement, the first-order inclusion probability of the \(k\)-th element is given by: \[\begin{equation} \pi_k=1-(1-p_k)^m \end{equation}\]
Proof.
Since these are independent events with an associated probability of success (success being equivalent to the event that the element \(k \in s\)) \(p_k\), each of these random draws is determined by a Bernoulli-type probability distribution. Consequently, when \(m\) independent trials are performed, the binomial probability distribution is used to find the first-order inclusion probabilities of each element in the population. \[\begin{align*} \pi_k=Pr(k\in S)&= 1-Pr(k\notin s)\\ &=1-\binom{m}{m}(1-p_k)^m(p_k)^{m-m}\\ &=1-(1-p_k)^m \end{align*}\]
TipResult
In sampling with replacement, the second-order inclusion probabilities \(\pi_{kl}\) are given by: \[\begin{equation} \pi_{kl}=1-(1-p_k)^m-(1-p_l)^m+(1-p_k-p_l)^m \ \ \ \ \ k\neq l=1 \ldots, N \end{equation}\]
Proof.
To find this probability, we must negate \((k\in S \text{ and }l\in s)\). This negation results in \((k\notin s \text{ or } l\notin s)\). Suppose we have two events, \(A=(k\notin s)\) and \(B=(l\notin s)\); therefore, \(Pr(A\cup B)=Pr(A)+Pr(B)-Pr(A\cap B)\). The preceding probabilities are governed by a binomial model, so: \[\begin{align*} \pi_{kl}&=Pr(k\in S \text{ and } l\in s)\\ &=1-Pr(k\notin s)-Pr(l\notin s)+Pr(k,l\notin s)\\ &=1-(1-p_k)^m-(1-p_l)^m+ \binom{m}{m}(1-p_k-p_l)^m(p_k+p_l)^{m-m}\\ &=1-(1-p_k)^m-(1-p_l)^m+(1-p_k-p_l)^m \end{align*}\] The fourth addend in the preceding equality is obtained by considering each trial as a Bernoulli process, where success is not choosing either \(k\) or \(l\). Therefore \[\begin{align*} Pr(\text{Success})&=1-Pr(\text{Failure})\\&=1-Pr(\text{Choosing $k$})-Pr(\text{Choosing $l$})+Pr(\text{Choosing both})\\ &=1-p_k-pl \end{align*}\] Since this is a single trial, the probability of choosing both is zero.
This is seen more clearly with the typical die example. If the event is rolling a die and success is not rolling 3 or 5, then the probability of success will be: \(1-Pr(\text{Failure})\), that is, \(1-Pr(\text{Rolling 5})-Pr(\text{Rolling 1})+Pr(\text{Rolling 5 and 1})\). It is obvious that the last addend is zero because this is a single roll.
NoteExample
The reader should not confuse the concept of sample with replacement with the concept of ordered draw. In our example population, the population size is \(N=5\). If a sampling design is used that induces fixed-size samples equal to \(m=2\), then there would be \(N^m=5^2=25\) possible ordered draws. However, there are only \(\binom{N+m-1}{m}=\binom{6}{2}=15\) possible samples with replacement. This scenario is easily shown with the help of the random variable \(n_k(S)\). The possible ordered draws are given as follows.
(1,1) (2,1) (3,1) (4,1) (5,1)
(1,2) (2,2) (3,2) (4,2) (5,2)
(1,3) (2,3) (3,3) (4,3) (5,3)
(1,4) (2,4) (3,4) (4,4) (5,4)
(1,5) (2,5) (3,5) (4,5) (5,5)However, although all possible ordered draws do not constitute the sampling support, they do help define it. In fact, the first step in constructing the sampling support with replacement is determining all possible draws. The OrderWR function from the TeachingSampling package makes it possible to know all possible fixed-size draws for a sampling design with replacement.
This function has three arguments: the first corresponds to the population size N; the second corresponds to the number of selections, m, which does not necessarily have to be less than the population size; and the last corresponds to a characteristic ID, which may be a set of labels or any other type of continuous identifier. The result of the OrderWR function will be a set of all possible ordered draws with fixed size m. When the argument ID is different from FALSE, the output of the function will correspond to the label or continuous identifier for each element of the population. In the following example, this function is used in our example population \(U\).
N <- length(U)
N[1] 5m <- 2
OrderWR(N, m, ID = FALSE) [,1] [,2]
[1,] 1 1
[2,] 1 2
[3,] 1 3
[4,] 1 4
[5,] 1 5
[6,] 2 1
[7,] 2 2
[8,] 2 3
[9,] 2 4
[10,] 2 5
[11,] 3 1
[12,] 3 2
[13,] 3 3
[14,] 3 4
[15,] 3 5
[16,] 4 1
[17,] 4 2
[18,] 4 3
[19,] 4 4
[20,] 4 5
[21,] 5 1
[22,] 5 2
[23,] 5 3
[24,] 5 4
[25,] 5 5OrderWR(N, m, ID = U) [,1] [,2]
[1,] "Yves" "Yves"
[2,] "Yves" "Ken"
[3,] "Yves" "Erik"
[4,] "Yves" "Sharon"
[5,] "Yves" "Leslie"
[6,] "Ken" "Yves"
[7,] "Ken" "Ken"
[8,] "Ken" "Erik"
[9,] "Ken" "Sharon"
[10,] "Ken" "Leslie"
[11,] "Erik" "Yves"
[12,] "Erik" "Ken"
[13,] "Erik" "Erik"
[14,] "Erik" "Sharon"
[15,] "Erik" "Leslie"
[16,] "Sharon" "Yves"
[17,] "Sharon" "Ken"
[18,] "Sharon" "Erik"
[19,] "Sharon" "Sharon"
[20,] "Sharon" "Leslie"
[21,] "Leslie" "Yves"
[22,] "Leslie" "Ken"
[23,] "Leslie" "Erik"
[24,] "Leslie" "Sharon"
[25,] "Leslie" "Leslie"Note that the set of ordered draws contains the sampling support with replacement. However, with the help of the SupportWR function from the TeachingSampling package, the true support induced by the with-replacement sampling design is defined. The arguments of this function are the same three arguments of the OrderWR function: N, m, and ID. The result of the function is the set of all possible fixed-size samples with replacement. For this particular example, the support is given by the following samples and not by all possible ordered draws.
SupportWR(N, m, ID = FALSE) [,1] [,2]
[1,] 1 1
[2,] 1 2
[3,] 1 3
[4,] 1 4
[5,] 1 5
[6,] 2 2
[7,] 2 3
[8,] 2 4
[9,] 2 5
[10,] 3 3
[11,] 3 4
[12,] 3 5
[13,] 4 4
[14,] 4 5
[15,] 5 5SupportWR(N, m, ID = U) [,1] [,2]
[1,] "Yves" "Yves"
[2,] "Yves" "Ken"
[3,] "Yves" "Erik"
[4,] "Yves" "Sharon"
[5,] "Yves" "Leslie"
[6,] "Ken" "Ken"
[7,] "Ken" "Erik"
[8,] "Ken" "Sharon"
[9,] "Ken" "Leslie"
[10,] "Erik" "Erik"
[11,] "Erik" "Sharon"
[12,] "Erik" "Leslie"
[13,] "Sharon" "Sharon"
[14,] "Sharon" "Leslie"
[15,] "Leslie" "Leslie"Of course, each possible sample with replacement that belongs to the support has different selection probabilities depending on the configuration of the individual selection probabilities for each element, \(p_k\). Suppose that each of the five elements of the population has selection probabilities given by
\[\begin{equation*} p_k= \begin{cases} 1/4, &\text{for $k=\textbf{Yves, Ken, Leslie}$},\\ 1/8, &\text{for $k=\textbf{Sharon, Erik}$} \end{cases} \end{equation*}\]
Note that \(\sum_U p_k=1\). For this particular configuration, and following equation 2.9, the selection probabilities \(p(s)\) of the samples in the support and the value of the variable \(n_k(S)\) would be given by the configuration shown in table 2.2, which is produced by the following code.
pk <- c(0.25, 0.25, 0.125, 0.125, 0.25)
QWR <- SupportWR(N, m, ID = U)
pWR <- p.WR(N, m, pk)
nkWR <- nk(N, m)
SamplesWR <- data.frame(QWR, pWR, nkWR)| 1 | 2 | 3 | n1 | n2 | n3 | n4 | n5 |
|---|---|---|---|---|---|---|---|
| Yves | Yves | 0.06 | 2 | 0 | 0 | 0 | 0 |
| Yves | Ken | 0.13 | 1 | 1 | 0 | 0 | 0 |
| Yves | Erik | 0.06 | 1 | 0 | 1 | 0 | 0 |
| Yves | Sharon | 0.06 | 1 | 0 | 0 | 1 | 0 |
| Yves | Leslie | 0.13 | 1 | 0 | 0 | 0 | 1 |
| Ken | Ken | 0.06 | 0 | 2 | 0 | 0 | 0 |
| Ken | Erik | 0.06 | 0 | 1 | 1 | 0 | 0 |
| Ken | Sharon | 0.06 | 0 | 1 | 0 | 1 | 0 |
| Ken | Leslie | 0.13 | 0 | 1 | 0 | 0 | 1 |
| Erik | Erik | 0.02 | 0 | 0 | 2 | 0 | 0 |
| Erik | Sharon | 0.03 | 0 | 0 | 1 | 1 | 0 |
| Erik | Leslie | 0.06 | 0 | 0 | 1 | 0 | 1 |
| Sharon | Sharon | 0.02 | 0 | 0 | 0 | 2 | 0 |
| Sharon | Leslie | 0.06 | 0 | 0 | 0 | 1 | 1 |
| Leslie | Leslie | 0.06 | 0 | 0 | 0 | 0 | 2 |
Note that the sum of the selection probabilities induced by the sampling design is equal to one and that each of them is greater than zero. The reader should notice that the sample belonging to the support is given in terms of \(n_k(S)\). Thus, if the seventh sample given by 1 0 1 0 0 has been selected, it actually does not matter whether Yves was selected before or after Erik, and the selection probability of this particular sample is 0.125 because
\[\begin{align*} p(s)&=\frac{2!}{1!0!1!0!0!}\left[ \left(\frac{1}{4}\right)^1\left(\frac{1}{4}\right)^0\left(\frac{1}{8}\right)^1 \left(\frac{1}{8}\right)^0\left(\frac{1}{4}\right)^0\right]\\ &=2\left(\frac{1}{32}\right)=0.0625 \end{align*}\]
2.2.2.2 Estimator of the Population Total
Hansen et al. (1953) propose a convenient estimator for the total of a population \(t_y\) when the sampling design is with replacement. The logic used in constructing this estimator is given below. Let the random event be:
This event defines the creation of random variables, which will be used later and whose behavior can be modeled through the following result.
TipResult
Let \(U_1,U_2,\ldots,U_m\) be a sequence of independent and identically distributed random variables with \(E(U_i)=\mu\) and \(Var(U_i)=\sigma^2\). Let \(\bar{U}=\sum_{i=1}^mU_i/m\). Then \(E(\bar{U})=\mu\), \(Var(\bar{U})=\sigma^2/m\), and an unbiased estimator of \(Var(\bar{U})\) is given by the following expression \[\begin{equation} \widehat{Var}(\bar{U})=\frac{1}{m(m-1)}\sum_{i=1}^m(U_i-\bar{U})^2 \end{equation}\] and consequently, an unbiased estimator for \(\sigma^2\) is given by \[\begin{equation} \hat{\sigma^2}=\frac{1}{m-1}\sum_{i=1}^m(U_i-\bar{U})^2. \end{equation}\]
Proof.
The expectation of \(\bar{U}\) is \[\begin{equation} E(\bar{U})=\frac{1}{m}\sum_{i=1}^mE(U_i)=\mu \end{equation}\] The variance is determined by \[\begin{equation} Var(\bar{U})=\frac{1}{m^2}\sum_{i=1}^mVar(U_i)=\sigma^2/m \end{equation}\] Note that the covariance terms are zero because the variables are independent of one another. Now, since \[\begin{equation} \sum_{i=1}^m(U_i-\bar{U})^2=\sum_{i=1}^mU_i^2-m\bar{U}^2 \end{equation}\] then, \[\begin{equation} E(\sum_{i=1}^m(U_i-\bar{U})^2)=\sum_{i=1}^mE(U_i^2)-mE(\bar{U}^2) \end{equation}\] On the other hand, \[\begin{align*} E(U_i^2)&=Var(U_i)+[E(U_i)]^2=\sigma^2+\mu^2\\ E(\bar{U}^2)&=Var(\bar{U})+[E(\bar{U})]^2=\sigma^2/m+\mu^2 \end{align*}\] This leads to the proof of the theorem because \[\begin{equation} E(\sum_{i=1}^m(U_i-\bar{U})^2)=(m-1)\sigma^2 \end{equation}\]
The preceding is a very powerful result that can be used for any type of independent and identically distributed random variables, and it will be the basis for proving results in parameter estimation using sampling designs with replacement. Continuing with the theoretical framework of sampling with replacement, we have the following definition.
ImportantDefinition
The random variable \(Z_i\) is defined such that \[\begin{equation} Z_i=y_{k_i}/p_{k_i} \ \ \ \ \ \ k\in U \ \ \ i=1,\ldots,m \end{equation}\] where the quantity \(y_{k_i}\) is the value of the characteristic of interest of the \(k\)-th element selected in the \(i\)-th draw. Analogously, \(p_{k_i}\) is the value of the selection probability of the \(k\)-th element selected in the \(i\)-th draw.
TipResult
The distribution of the random variable \(Z_i\) is given by \[\begin{equation} Pr\left(Z_i=\frac{y_k}{p_k}\right)=p_k, \end{equation}\] therefore the expectation and variance of the random variable \(Z_i\) are \[\begin{equation} E(Z_i)=t_y \end{equation}\] y \[\begin{equation} Var(Z_i)=\sum_Up_k\left(\frac{y_k}{p_k}-t_y\right)^2, \end{equation}\] respectively.
Proof.
Since these are \(m\) independent random draws, the random variable \(Z_i\) can take the following values \[\begin{equation*} \frac{y_1}{p_1},\frac{y_2}{p_2}\ldots,\frac{y_N}{p_N} \end{equation*}\] with probabilities \[\begin{equation*} p_1,p_2\ldots,p_N \end{equation*}\] respectively. Then, using the general definition of the expectation operator, we have \[\begin{align*} E(Z_i)=\sum_U\frac{y_k}{p_k}Pr\left(Z_i=\frac{y_k}{p_k}\right)=\sum_U\frac{y_k}{p_k}p_k=t_y \end{align*}\] and analogously we have the variance \[\begin{align*} Var(Z_i)=\sum_U\left(\frac{y_k}{p_k}-E(Z_i)\right)^2Pr\left(Z_i=\frac{y_k}{p_k}\right) =\sum_U\left(\frac{y_k}{p_k}-t_y\right)^2p_k \end{align*}\]
Since the \(m\) draws are independent events, the variables \(Z_i\) are also independent. Note that the quantity \(Z_i\) is an estimate of the population total with the \(i\)-th selected sample of size 1. Now, since there are \(m\) draws, there will be \(m\) estimates of the population total; therefore, as in many other statistical procedures, we use the average of these \(m\) estimates to obtain a unified estimate for \(t_y\). The Hansen-Hurwitz estimator takes the following form
\[ \hat{t}_{y,p}=\frac{1}{m}\sum_{i=1}^{m}\frac{y_{k_i}}{p_{k_i}} \tag{2.10}\]
To have a sampling strategy that is efficient in estimating \(t_y\), it is convenient to use the Hansen-Hurwitz estimator when the selection probabilities are proportional to the characteristic of interest; that is, when \(p_k\propto y_k\). If this occurs, the estimator will have nearly zero variance and the estimate will be very precise.
TipResult
If \(p_k>0\) for all \(k\in U\), the estimator \(\hat{t}_{y,p}\) is unbiased
Proof.
The random variables \(Z_i\) are independent (because each trial is independent) and their distribution is induced by \(Pr(Z_i=y_k/p_k)=p_k\), \(k \in U\); that is, they are identically distributed. Therefore, the Hansen-Hurwitz estimator can be written as: \[\begin{align*} \hat{t}_{y,p}=\frac{1}{m}\sum_{i=1}^{m}\frac{y_i}{p_i}=\frac{1}{m}\sum_{i=1}^{m}Z_i=\bar{Z} \end{align*}\] and thus, with \(p_k>0\) for all \(k\in U\), we have \[\begin{align*} E(\hat{t}_{y,p})=\frac{1}{m}\sum_{i=1}^{m}E(Z_i)=\frac{1}{m}\sum_{i=1}^{m}t_y=t_y \end{align*}\]
2.2.2.3 Variance of the Hansen-Hurwitz Estimator
One of the most important characteristics of the Hansen-Hurwitz estimator is the simplicity of its variance expression. This same feature means that, even though sampling is with replacement, the Hansen-Hurwitz estimator is frequently used by users of sampling studies.
TipResult
The variance of the Hansen-Hurwitz estimator is given by the following expression \[ Var(\hat{t}_{y,p})=\frac{1}{m}\sum_{k=1}^{N}p_k\left(\frac{y_k}{p_k}-t_y\right)^2 \tag{2.11}\]
Proof.
By the independence of the selections, we have \[\begin{align*} Var(\hat{t}_{y,p})&=Var\left(\frac{1}{m}\sum_{i=1}^{m}Z_i\right)\\ &=\frac{1}{m^2}\sum_{i=1}^{m}Var(Z_i)\\ &=\frac{1}{m}Var(Z_i)\\ &=\frac{1}{m}\sum_U\left(\frac{y_k}{p_k}-t_y\right)^2p_k \end{align*}\]
The preceding expression makes the computational calculation of the variance of the Hansen-Hurwitz estimator very simple. However, this variance can be written in several forms, some of them very useful for the theoretical development of the properties of the estimator.
TipResult
In general, the variance of the Hansen-Hurwitz estimator can be written as follows \[ Var(\hat{t}_{y,p})=\frac{1}{m}\left(\sum_{k=1}^{N}\frac{y_k^2}{p_k}-t_y^2\right) \tag{2.12}\]
Proof.
\[\begin{align*} Var(\hat{t}_{y,p})=&=\frac{1}{m}\sum_{k=1}^{N}p_k\left(\frac{y_k}{p_k}-t_y\right)^2\\ &=\frac{1}{m}\sum_{k=1}^{N}p_k\left(\frac{y_k^2}{p_k^2}-2t_y\frac{y_k}{p_k}+t_y^2\right)\\ &=\frac{1}{m}\sum_{k=1}^{N}\left(\frac{y_k^2}{p_k}-2t_yy_k+p_kt_y^2\right)\\ &=\frac{1}{m}\left(\sum_{k=1}^{N}\frac{y_k^2}{p_k}-2t_y\sum_{k=1}^{N}y_k+t_y^2\sum_{k=1}^{N}p_k\right)\\ &=\frac{1}{m}\left(\sum_{k=1}^{N}\frac{y_k^2}{p_k}-2t_y^2+t_y^2\right) =\frac{1}{m}\left(\sum_{k=1}^{N}\frac{y_k^2}{p_k}-t_y^2\right) \end{align*}\]
2.2.2.4 Variance Estimation
TipResult
An unbiased estimator of equation 2.11 is \[ \widehat{Var}(\hat{t}_{y,p})=\frac{1}{m(m-1)}\sum_{i=1}^{m}\left(\frac{y_i}{p_i}-\hat{t}_{y,p}\right)^2 \tag{2.13}\]
Proof.
By developing the variance of the estimator, we arrive at the fact that it is equal to \[\frac{1}{m}Var(Z_i).\] Now, using Result 2.2.11, since \(Z_1,\ldots,Z_m\) form a random sample of variables with expectation \(t_y\) and identical variance, a natural and unbiased estimator for the variance of \(Z_i\) is \[\frac{1}{m-1}\sum_{i=1}^{m}(Z_i-\bar{Z})^2=\frac{1}{m-1}\sum_{i=1}^{m}\left(\frac{y_i}{p_i}-\hat{t}_{y,p}\right)^2\] therefore, an unbiased estimator of the variance of the Hansen-Hurwitz estimator is \[\begin{align*} \widehat{Var}(\hat{t}_{y,p})=\frac{1}{m}\frac{1}{m-1}\sum_{i=1}^{m}\left(\frac{y_{k_i}}{p_{k_i}}-\hat{t}_{y,p}\right)^2 \end{align*}\]
TipResult
An alternative expression for variance estimation of the Hansen-Hurwitz estimator in sampling with replacement is \[\begin{align*} \widehat{Var}(\hat{t}_{y,p})=\frac{1}{m(m-1)}\sum_{i=1}^{m}\left(\frac{y_{k_i}}{p_{k_i}}\right)^2-m\hat{t}_{y,p}^2 \end{align*}\]
Proof.
Starting from the preceding result, we have \[\begin{align*} m(m-1)\widehat{Var}(\hat{t}_{y,p}) &=\sum_{i=1}^{m}\left(\frac{y_{k_i}}{p_{k_i}}-\hat{t}_{y,p}\right)^2\\ &=\sum_{i=1}^{m}\left(\frac{y_{k_i}^2}{p_{k_i}^2}-2\hat{t}_{y,p}\frac{y_{k_i}}{p_{k_i}}+\hat{t}_{y,p}^2\right)\\ &=\sum_{i=1}^{m}\left(\frac{y_{k_i}^2}{p_{k_i}^2}\right)-2\hat{t}_{y,p}\sum_{i=1}^{m}\frac{y_{k_i}}{p_{k_i}}+m\hat{t}_{y,p}\\ &=\sum_{i=1}^{m}\left(\frac{y_{k_i}^2}{p_{k_i}^2}\right)-2m\hat{t}_{y,p}^2+m\hat{t}_{y,p}\\ &=\sum_{i=1}^{m}\left(\frac{y_{k_i}}{p_{k_i}}\right)^2-m\hat{t}_{y,p}^2 \end{align*}\]
Even if the sampling design is with replacement, it is possible to use the Horvitz-Thompson estimator, since it preserves its unbiasedness. The comparison between the precision of the Horvitz-Thompson estimator and the Hansen-Hurwitz estimator in a design with repetition depends on the configuration of the values of the characteristic of interest in the population, \(y_k\) \(\forall k=1,2,...,N\). However, the Horvitz-Thompson estimator is generally more efficient than the Hansen-Hurwitz estimator, although the latter is easier to calculate. When the sampling design is fixed-size, the Horvitz-Thompson and Hansen-Hurwitz estimators coincide.
NoteExample
Continuing with the lexical-graphic exercise of estimating the population total \(t_y\) for all possible samples with replacement of size 2 from population U, we have the following table showing the sampling support with the help of the SupportWR function.
all.y <- SupportWR(N, n, y)
all.pk <- SupportWR(N, n, pk)
all.HH <- rep(0, 15)
for (k in 1:15) {
all.HH[k] <- HH(all.y[k, ], all.pk[k, ])
}
AllSamplesWR <- data.frame(QWR, all.pk, pWR, all.y, all.HH)| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|
| Yves | Yves | 0.25 | 0.25 | 0.06 | 32 | 32 | 128 |
| Yves | Ken | 0.25 | 0.25 | 0.13 | 32 | 34 | 132 |
| Yves | Erik | 0.25 | 0.12 | 0.06 | 32 | 46 | 248 |
| Yves | Sharon | 0.25 | 0.12 | 0.06 | 32 | 89 | 420 |
| Yves | Leslie | 0.25 | 0.25 | 0.13 | 32 | 35 | 134 |
| Ken | Ken | 0.25 | 0.25 | 0.06 | 34 | 34 | 136 |
| Ken | Erik | 0.25 | 0.12 | 0.06 | 34 | 46 | 252 |
| Ken | Sharon | 0.25 | 0.12 | 0.06 | 34 | 89 | 424 |
| Ken | Leslie | 0.25 | 0.25 | 0.13 | 34 | 35 | 138 |
| Erik | Erik | 0.12 | 0.12 | 0.02 | 46 | 46 | 368 |
| Erik | Sharon | 0.12 | 0.12 | 0.03 | 46 | 89 | 540 |
| Erik | Leslie | 0.12 | 0.25 | 0.06 | 46 | 35 | 254 |
| Sharon | Sharon | 0.12 | 0.12 | 0.02 | 89 | 89 | 712 |
| Sharon | Leslie | 0.12 | 0.25 | 0.06 | 89 | 35 | 426 |
| Leslie | Leslie | 0.25 | 0.25 | 0.06 | 35 | 35 | 140 |
The vector Est contains the Hansen-Hurwitz estimates for each of the 15 possible samples with replacement; its expectation is calculated as
sum(all.HH * pWR)[1] 236Note that the expectation of the estimator is equal to the total of the characteristic of interest, confirming its unbiasedness. On the other hand, to select a sample with replacement, R includes the sample function, whose main arguments are
x is the population size, and size is an integer that determines the sample size. To select a sample with replacement, the argument replace must take the value TRUE, that is, replace = TRUE. Each element belonging to the population must have an associated vector of selection probabilities whose sum is equal to one. In R, the argument prob contains this vector of probabilities; when this argument is omitted, the sample function assumes that the selection probabilities are identical for each individual in the population. Thus, for example, to select a sample with replacement from the sampling frame of \(U\) of size \(m=3\), with selection probabilities given by
pk[1] 0.25 0.25 0.12 0.12 0.25Note that the sum of the selection probabilities is equal to one and that the labels or names for each individual in the population are contained in the object U.
U[1] "Yves" "Ken" "Erik" "Sharon" "Leslie"To select a sample with replacement of size \(m=3\), the following code must be written.
sam <- sample(N, 3, replace = TRUE, prob = pk)
sam[1] 2 4 3For the preceding selection, the first element was chosen twice and the third element once. The indexing of the labels (names) and values of the characteristic of interest of the elements chosen in the sample is done using
pkm <- pk[sam]
pkm[1] 0.25 0.12 0.12ym <- y[sam]
ym[1] 34 89 46Note that the sample size is 3, but the effective sample size is \(n(S)=2\). Here pkm is the vector of selection probabilities for the individuals belonging to the sample, and ym is the vector of values of the characteristic of interest for the individuals belonging to the sample. The HH function from the TeachingSampling package estimates the population total for the characteristic of interest. This function has two arguments: y, the vector of values of the characteristic of interest for the individuals in the sample, and pk, their corresponding selection probabilities.
est <- HH(ym, pkm)[1]
est[1] 405To estimate the variance, a vector of differences dif is created between \(\frac{y_i}{p_i}\) and the estimate. Then it is squared, summed, and divided by \(m(m-1)\).
dif <- rep(0, 3)
dif[1] <- (ym[1] / pkm[1]) - est
dif[2] <- (ym[2] / pkm[2]) - est
dif[3] <- (ym[3] / pkm[3]) - est
dif[1] -269 307 -37Var <- (1 / 3) * (1 / 2) * sum(dif^2)
Var[1] 27996sqrt(Var)[1] 167Then, the corresponding estimated coefficient of variation is
\[\begin{equation*} cve(\hat{t}_{p})=\frac{167.32}{405.33}\cong 41\% \end{equation*}\]
Note that using the HH function gives the same result from the procedure.
HH(ym, pkm) y
Estimation 405
Standard Error 167
CVE 41We can think of the estimated coefficient of variation as a measure of precision. Thus, the preceding estimates could be considered unacceptable because this measure is very high.
The objective of this book is for the reader to be able to propose sampling strategies that allow precise and reliable estimates. That is, estimates whose coefficient of variation is acceptable and whose confidence interval length is short with a satisfactory confidence level.
2.2.3 The Horvitz-Thompson Estimator in With-Replacement Designs
2.3 Representative Samples
Sampling theory has been enriched in recent decades by valuable contributions worldwide; although the basis of sampling theory is probability theory, whose axiomatic development is several hundred years old, its practical development did not occur until the beginning of the twentieth century. However, in classical statistical inference, based on the thinking of Ronald Fisher and others, the population is assumed to be infinite. A fundamental aspect of sampling theory is that it is based on reality, where populations, no matter how large, are finite in nature.
Starting from this fact, it is possible to ground inference based on a random sample that comes from a finite population, and from this perspective the results of the inferences will differ significantly. In fact, this is a warning for people who make inferences with data from a sampling study to update their knowledge and avoid major mistakes when presenting inferential results (Chambers and Skinner 2003). For this reason, sampling theory covers fundamental aspects of statistics, because from a controlled experiment to a sample survey, one must think about the mechanism for collecting information and, from there, about inference.
A common example in classrooms is to describe the population on the board with a happy face; the teacher says that a representative sample of the population is a sample in which the same happy face can still be seen. That is, there is a belief that a representative sample is a reduced model of the population, and from this comes an argument about the validity of the sample: a good sample is one that resembles the population, so that categories appear in the same proportions as in the population. Nothing could be more false than this belief. In some cases, it is essential to overrepresent some categories or even select units with unequal probabilities.
Till’e (2006) cites the following example: suppose the objective is to estimate iron production in a country and we know that iron is produced by two giant companies with thousands of employees and by hundreds of small companies with few employees. Is the best way to select the sample to assign the same probability to each company? Of course not. First we find out the production of the large companies. Then we select a sample of the small companies.
The sample should not be a reduced model of the population; it should be a tool used to obtain estimates. Thus, the concept of a representative sample loses weight. Moreover, for H’ajek (1981), a sampling strategy is a pair: a sampling design (a probability distribution over all possible samples) and an estimator. Sampling theory has studied optimal strategies that make it possible to ensure the quality of estimates. Therefore, the concept of representativeness should be associated with sampling strategies and not only with samples.
Following Till’e (2006), a strategy is said to be representative if it makes it possible to estimate a population total exactly; that is, without bias and with zero variance. If, for example, the Horvitz-Thompson estimator is used together with an appropriate sampling design, this strategy is representative only if, together with the selected sample, the estimator reproduces some population totals; such samples are called balanced samples. There are also estimators that give the strategy the qualification of representative, some of which are known as calibration estimators.
2.4 Exercises
Prove that under a sampling design \(p(s)\), the mean squared error of any estimator \(\hat{T}(s)\) of a parameter \(T\) is equal to the variance \(Var(\hat{T})\) plus the squared bias \(B^2(\hat{T})\).
Hint: \(MSE\left(\hat{T}\right)=E_p\left(\hat{T}(s)-T\right)^2=\sum_{s\in Q}\left(\hat{T}(s)-T\right)^2p(s)\).
Show that \(\pi_{kl}=E_p \left( I_k(s) I_l(s) \right)\).
Suppose you have access to the finite population of size \(N=5\) from Example 2.2.1 and assume the following sampling design without replacement \[\begin{equation*} p(S=s)= \begin{cases} 0.2, & \text{for $s=\{Ken, Erik, Sharon\}$, $s=\{Ken, Leslie\}$},\\ 0.3, & \text{for $s=\{Yves, Erik, Leslie\}$, $s=\{Yves, Sharon\}$},\\ 0, & \text{otherwise}. \end{cases} \end{equation*}\]
- Calculate all first- and second-order inclusion probabilities.
- Is the preceding a fixed sample-size sampling design? Explain.
- List all values taken by the random variable \(n(S)\) and verify the relationships \(E_p(n(S))=\sum_U\pi_k\) and \(Var_p(n(S))=\sum_U\pi_k-\left(\sum_U\pi_k\right)^2+\sum\sum_{k\neq l}\pi_{kl}\).
- Suppose you have access to the finite population of size \(N=5\) from Example 2.2.1 and assume the following sampling design without replacement \[\begin{equation*} p(S=s)= \begin{cases} 0.1, & \text{if $n(S)=3$},\\ 0, & \text{otherwise}. \end{cases} \end{equation*}\]
- Define all possible samples that belong to the support induced by the preceding sampling design.
- Calculate all first- and second-order inclusion probabilities.
- Verify that \(\sum_U\pi_k=3\) and that \(\sum_U\pi_k-\left(\sum_U\pi_k\right)^2+\sum\sum_{k\neq l}\pi_{kl}=0\). Explain.
- Verify that \(\sum_U\pi_{k1}=3\times \pi_1\), \(\sum_U\pi_{k2}=3\times \pi_2\), through \(\sum_U\pi_{k5}=3\times \pi_5\).
- Calculate all possible covariances \(\Delta_{kl}\) and verify that \(\sum_U\Delta_{k1}=0\), through \(\sum_U\Delta_{k5}=0\).
- Prove or refute the following statement: “Under any sampling design, the population sum of the first-order inclusion probabilities is always equal to the sample size.”
- Prove or refute the following statement: “Under any sampling design, the Horvitz-Thompson estimator can be used to obtain an unbiased estimate of the population total.”
- Suppose you have access to the finite population of size \(N=5\) from Example 2.2.1 and that \(y_k\) denotes the value of the characteristic of interest for the \(k\)-th individual. Thus, we have: \[y_{Yves}=32, \ \ y_{Ken}=34, \ \ y_{Erik}=46, \ \ y_{Sharon}=89, \ \ y_{Leslie}=35\]
- For the sampling design in Exercise 2.3, for each possible sample calculate the Horvitz-Thompson estimate, the variance estimate, the \(cve\), and the 95% confidence interval estimate. Finally, show that the estimator is unbiased and calculate the variance of the estimator using equation 2.5.
- For the sampling design in Exercise 2.4, for each possible sample calculate the Horvitz-Thompson estimate, the variance estimate, the \(cve\), and the 95% confidence interval estimate. Finally, show that the estimator is unbiased and calculate the variance of the estimator using equation 2.5 and equation 2.6. Are these variances equal? Explain.
- For the sampling design in Exercise 2.3, for each possible sample calculate the Horvitz-Thompson estimate of the mean (expression 2.2.10), the estimate of the population size (expression 2.2.14), the alternative estimate of the mean (expression 2.2.15), and the alternative estimate of the total (expression 2.2.18).
- For the sampling design in Exercise 2.4, for each possible sample calculate the Horvitz-Thompson estimate of the mean (expression 2.2.10), the estimate of the population size (expression 2.2.14), the alternative estimate of the mean (expression 2.2.15), and the alternative estimate of the total (expression 2.2.18).
- Prove or refute the following statement: “Under any sampling design with replacement, the Hansen-Hurwitz estimator can be used to obtain an unbiased estimate of the population total.”
- Prove or refute the following statement: “The selection probability of an individual is always equal to its inclusion probability.”
- Prove or refute the following statement: “Any sampling design with replacement can be seen as a particular case of the multinomial distribution.”
- Prove or refute the following statement: “For a population of size \(N\), the number of possible samples with replacement of size \(m\) is \(N^m\).”
- Suppose you have access to the finite population of size \(N=5\) from the preceding exercises and assume the following selection probabilities \[\begin{equation*} p_k= \begin{cases} 0.3, & \text{for $k=Yves, Leslie$},\\ 0.2, & \text{for $k=Erik$},\\ 0.1, & \text{for $k=Ken, Sharon$}. \end{cases} \end{equation*}\]
- How many samples with replacement of size \(m=3\) can be selected? Explicitly specify the sampling design for these samples and verify that \(\sum_{s\in Q}p(s)=1\).
- For this sampling design, and taking into account the values of the characteristic of interest from Exercise 2.7, for each possible sample calculate the Hansen-Hurwitz estimate, the variance estimate, the \(cve\), and the 95% confidence interval estimate. Finally, show that the estimator is unbiased and calculate the variance of the estimator using equation 2.11.
- Is it possible to use another type of estimator to obtain unbiased estimates of the population total?
- Rigorously prove that the variance estimator of the Hansen-Hurwitz estimator corresponds to equation 2.13.
Bautista, J. L. 1998. Diseños de Muestreo Estad’istico. Universidad Nacional de Colombia.
Brewer, K. R. W. 2002. Combined Sampling Inference, Weighting Basu’s Elephants. London: Arnorld.
Cassel, C. M., C. E. Särndal, and J. Wretman. 1976. Foundations of Inference in Survey Sampling. Wiley.
Chambers, R. L., and C. J. Skinner, eds. 2003. Analysis of Survey Data. Wiley.
H’ajek, J. 1960. “Limiting Distributions in Simple Random Sampling from a Finite Poulation.” Publication of Mathematical Institute of the Hungarian Academy of Science 5: 361–74.
H’ajek, J. 1981. Sampling from a Finite Population. New York: Marcel Dekker.
Hansen, M. H., W. N. Hurwitz, and G. Madow W. 1953. Sample Survey Methods and Theory. Vols. I and II. John Wiley; Sons.
Horvitz, D. G., and D. J. Thompson. 1952. “A Generalization of Sampling Without Replacement from a Finite Universe.” Journal of the American Statistical Association 47: 663–85.
Lohr, S. 2000. Sampling: Design and Analysis. Thompson.
Narain, R. D. 1951. “On Sampling Without Replacement with Varying Probabilities.” Journal of Indian Society of Agricultural Statistics 3: 169–75.
Särndal, C. E., B. Swensson, and J. Wretman. 1992. Model Assisted Survey Sampling. Springer, New York.
Sen, A. R. 1953. “On the Estimate of the Variance in Sampling with Varying Probabilities.” Journal of the Indian Society of Agricultural Statistics 5: 119–27.
Till’e, Y. 2006. Sampling Algorithms. Springer.
Yates, F., and P. M. Grundy. 1953. “Selecting Withou Replacement from Within Estrata with Probability Proportional to Size.” Journal of the Royal Statitical Society B15: 235–61.