Also assumed, and not repeated: Unit 1 for the survey vocabulary, census against sample, and sampling against non-sampling error; Unit 2 for simple random sampling with and without replacement, \(\operatorname{Var}(\bar y) = \frac{N-n}{Nn}S^{2}\), and the finite population correction; and Unit 3 for stratification and allocation.
This unit begins with the four things that section does not do: the variance of the Hansen–Hurwitz estimator and how to estimate it; Lahiri's method actually worked, with a proof that it selects proportionally to size; sampling without replacement, where the inclusion probabilities \(\pi_i\) and \(\pi_{ij}\) have to be computed rather than assumed; and the Horvitz–Thompson estimator, which is the general answer to every unequal-probability design.
Quoting the estimator from the Foundation section, so that the derivation has something to work on: \(n\) units are drawn independently, unit \(i\) with probability \(p_i\) where \(\sum_i p_i = 1\), and
\[ \hat Y_{HH} = \frac{1}{n}\sum_{r=1}^{n}\frac{y_r}{p_r}, \]the sum being over the \(n\) draws, a unit counted as often as it is drawn.
Unbiasedness, in one line. That section asserts this without proof; the proof is short. Each draw is an independent trial, so for any one of them
\[ E\!\left(\frac{y}{p}\right) = \sum_{i=1}^{N} p_i \cdot \frac{Y_i}{p_i} = \sum_{i=1}^{N} Y_i = Y, \]and averaging \(n\) such terms leaves \(Y\). The division by \(p_i\) is doing exactly the work of undoing the unequal selection: a unit \(k\) times as likely to be chosen contributes \(1/k\) as much when it is.
The variance. Since the draws are independent and identically distributed,
\[ \operatorname{Var}\!\left(\hat Y_{HH}\right) = \frac{1}{n}\operatorname{Var}\!\left(\frac{y}{p}\right) = \frac{1}{n}\sum_{i=1}^{N} p_i\left(\frac{Y_i}{p_i} - Y\right)^{2}. \]Read the bracket: it is the amount by which unit \(i\)'s scaled value departs from the population total. If \(Y_i\) were exactly proportional to \(p_i\), every bracket would be zero and the variance would vanish. That is the whole design principle: choose the size measure \(X_i\) so that \(Y_i/X_i\) is as nearly constant as possible.
An unbiased estimator of the variance, from the sample alone:
\[ \widehat{\operatorname{Var}}\!\left(\hat Y_{HH}\right) = \frac{1}{n(n-1)}\sum_{r=1}^{n}\left(\frac{y_r}{p_r} - \hat Y_{HH}\right)^{2}, \]which is the ordinary sample variance of the \(n\) quantities \(y_r/p_r\), divided by \(n\) — because \(\hat Y_{HH}\) is their mean and the draws are independent. This is the one real convenience of sampling with replacement, and section 3 shows what is lost by accepting it.
Given. A population of \(N = 5\) units with a known size measure \(X_i\) and (for this demonstration only) known values \(Y_i\):
| Unit | 1 | 2 | 3 | 4 | 5 | Total |
|---|---|---|---|---|---|---|
| \(X_i\) | 10 | 20 | 30 | 25 | 15 | 100 |
| \(Y_i\) | 12 | 26 | 33 | 30 | 19 | 120 |
| \(p_i = X_i/100\) | 0.10 | 0.20 | 0.30 | 0.25 | 0.15 | 1.00 |
Step 1 — the scaled values.
\[ \frac{Y_i}{p_i}: \quad \frac{12}{0.10} = 120, \quad \frac{26}{0.20} = 130, \quad \frac{33}{0.30} = 110, \quad \frac{30}{0.25} = 120, \quad \frac{19}{0.15} = 126.666667. \]Step 2 — the departures from \(Y = 120\).
\[ 0, \qquad 10, \qquad -10, \qquad 0, \qquad 6.666667. \]Two of the five are exactly zero, because for those units \(Y_i/X_i = 1.2\) precisely.
Step 3 — the variance, at \(n = 2\).
\[ \sum_i p_i\left(\frac{Y_i}{p_i} - Y\right)^{2} = 0.20(100) + 0.30(100) + 0.15(44.444444) = 20 + 30 + 6.666667 = 56.666667, \] \[ \operatorname{Var}\!\left(\hat Y_{HH}\right) = \frac{1}{2}\cdot\frac{170}{3} = \frac{85}{3} = 28.333333. \]Step 4 — what equal probabilities would have cost. Under simple random sampling with replacement the estimator of the total is \(N\bar y\), with variance \(N^{2}\sigma^{2}/n\). Here \(\bar Y = 24\) and
\[ \sigma^{2} = \frac{(-12)^{2} + 2^{2} + 9^{2} + 6^{2} + (-5)^{2}}{5} = \frac{144 + 4 + 81 + 36 + 25}{5} = \frac{290}{5} = 58, \] \[ \operatorname{Var}(N\bar y) = \frac{25 \times 58}{2} = 725. \]Step 5 — compare.
\[ \frac{725}{85/3} = \frac{2175}{85} = 25.588235. \]Step 6 — an actual estimate. If the two draws fall on units 2 and 3,
\[ \hat Y_{HH} = \frac{1}{2}\left(130 + 110\right) = 120, \]the true total exactly; if they fall on units 3 and 5, \(\tfrac12(110 + 126.666667) = 118.333333\).
Interpretation. Selecting proportionally to size has made the estimator more than twenty-five times as precise as equal-probability sampling — equivalent to multiplying the sample size by twenty-five. The gain is that large because \(Y_i\) is nearly proportional to \(X_i\) here; it is not a general figure. The corresponding warning is also worth stating: if \(Y_i\) were inversely related to \(X_i\), PPS would be worse than equal-probability sampling, and the size measure would have been chosen badly.
The cumulative-total method of the Foundation section requires the running totals \(X_1, X_1+X_2, \ldots\) to be formed before any selection can be made. For \(N\) in the thousands that is a substantial piece of arithmetic, and it must be redone whenever the frame changes. Lahiri's method avoids it entirely, at the cost of some rejected trials.
The procedure. Let \(M = \max_i X_i\). Then repeat:
Proof that it selects with probability proportional to size. On any one trial the two draws are independent, so
\[ P(\text{trial accepts and gives unit } i) = \frac{1}{N}\cdot\frac{X_i}{M}, \]and the probability that a trial accepts at all is the sum of these,
\[ P(\text{accept}) = \sum_{k=1}^{N}\frac{X_k}{NM} = \frac{\sum_k X_k}{NM}. \]Given that some trial eventually accepts, the conditional probability that it delivered unit \(i\) is the ratio,
\[ P(\text{select } i) = \frac{X_i/(NM)}{\left(\sum_k X_k\right)/(NM)} = \frac{X_i}{\sum_k X_k} = p_i. \quad\blacksquare \]No cumulation appears anywhere in that argument, which is the point.
The cost. The number of trials until an acceptance is geometric with success probability \(\sum_k X_k/(NM)\), so on average
\[ E(\text{trials per selection}) = \frac{NM}{\sum_k X_k} = \frac{M}{\bar X}, \]the ratio of the largest size to the mean size. The method is efficient when the sizes are comparable and wasteful when one unit dwarfs the rest — if a single unit is a hundred times the average, about a hundred trials are needed per selection, and cumulating would have been cheaper.
Given. The five units of Example 1.1, with \(X = (10, 20, 30, 25, 15)\).
Step 1 — the two constants. \(N = 5\), \(M = \max X_i = 30\).
Step 2 — the acceptance probability.
\[ \frac{\sum_k X_k}{NM} = \frac{100}{5 \times 30} = \frac{100}{150} = \frac{2}{3} = 0.666667. \]Two trials in three succeed, so
\[ E(\text{trials per selection}) = \frac{3}{2} = 1.5. \]Step 3 — three trials.
| Trial | Unit \(i\) | \(R\) | \(X_i\) | \(R \le X_i\)? | Outcome |
|---|---|---|---|---|---|
| 1 | 3 | 22 | 30 | yes | select unit 3 |
| 2 | 1 | 26 | 10 | no | discard, repeat |
| 3 | 5 | 9 | 15 | yes | select unit 5 |
Two selections from three trials, against an expectation of \(1.5\) trials each.
Interpretation. Notice what unit 1 suffers: with \(X_1 = 10\) against \(M = 30\), two trials in three that land on it are thrown away. That is exactly right — unit 1 should be selected a third as often as unit 3 — and it is achieved without a single cumulative total being formed.
Sampling with replacement can select the same unit twice, which wastes an observation and is indefensible when \(n\) is an appreciable fraction of \(N\). Sampling without replacement fixes that, and immediately destroys the convenient structure of section 1: the draws are no longer independent, so there is no product to take expectations over and no sample variance to fall back on.
Everything is instead expressed through two sets of numbers, which are properties of the design and not of the data:
\[ \pi_i = P(\text{unit } i \text{ is in the sample}), \qquad \pi_{ij} = P(\text{units } i \text{ and } j \text{ are both in the sample}). \]Two identities hold for any design of fixed size \(n\), and both are worth checking on any design before using it:
\[ \sum_{i=1}^{N}\pi_i = n, \qquad \sum_{j \ne i}\pi_{ij} = (n-1)\,\pi_i, \qquad \sum_{i<j}\pi_{ij} = \frac{n(n-1)}{2}. \]Why the first. Write \(a_i\) for the indicator that unit \(i\) is in the sample. Then \(\sum_i a_i = n\) always, and taking expectations gives \(\sum_i \pi_i = n\) since \(E(a_i) = \pi_i\). Why the second. Condition on unit \(i\) being in the sample; the other \(n-1\) places are filled by other units, so \(\sum_{j\ne i} P(j \in s \mid i \in s) = n-1\), and multiplying by \(\pi_i\) gives the result.
Computing them is design-specific. For successive draws with probabilities \(p_i\), renormalised after each draw, \(n = 2\) gives
\[ \pi_i = p_i + \sum_{j \ne i}p_j\,\frac{p_i}{1 - p_j}, \qquad \pi_{ij} = \frac{p_ip_j}{1-p_i} + \frac{p_jp_i}{1-p_j}, \]the first term of \(\pi_i\) being selection on the first draw and the sum being selection on the second after some other unit took the first. Note what these are not: \(\pi_i \ne np_i\). The inclusion probabilities are only approximately proportional to size, and the discrepancy grows with \(n/N\) — which is why designs that make \(\pi_i = np_i\) exactly, such as Midzuno's and Brewer's, were invented.
For any design with known \(\pi_i > 0\), the Horvitz–Thompson estimator of the total is
\[ \hat Y_{HT} = \sum_{i \in s}\frac{y_i}{\pi_i}, \]the sum over the distinct units in the sample, each weighted by the reciprocal of its inclusion probability.
Unbiasedness. Write it as a sum over the whole population using the indicator \(a_i\):
\[ \hat Y_{HT} = \sum_{i=1}^{N} a_i\,\frac{Y_i}{\pi_i} \;\Longrightarrow\; E\!\left(\hat Y_{HT}\right) = \sum_{i=1}^{N}\pi_i\,\frac{Y_i}{\pi_i} = \sum_{i=1}^{N} Y_i = Y. \]The proof needs nothing about the design beyond \(\pi_i > 0\) for every unit — which is the condition that no unit is unreachable, and the only condition there is.
The variance, from \(\operatorname{Var}(a_i) = \pi_i(1-\pi_i)\) and \(\operatorname{Cov}(a_i, a_j) = \pi_{ij} - \pi_i\pi_j\):
\[ \operatorname{Var}\!\left(\hat Y_{HT}\right) = \sum_i \frac{1-\pi_i}{\pi_i}Y_i^{2} + 2\sum_{i<j}\frac{\pi_{ij} - \pi_i\pi_j}{\pi_i\pi_j}\,Y_iY_j. \]The Yates–Grundy form. For a design of fixed sample size the identities of section 3 allow this to be rewritten as a single sum of squares:
\[ \operatorname{Var}\!\left(\hat Y_{HT}\right) = \sum_{i<j}\left(\pi_i\pi_j - \pi_{ij}\right) \left(\frac{Y_i}{\pi_i} - \frac{Y_j}{\pi_j}\right)^{2}. \]It is much the better form, for three reasons. It is manifestly zero when every \(Y_i/\pi_i\) is equal, which is the design ideal. Its unbiased estimator
\[ \widehat{\operatorname{Var}}_{YG} = \sum_{i<j \in s}\frac{\pi_i\pi_j - \pi_{ij}}{\pi_{ij}} \left(\frac{y_i}{\pi_i} - \frac{y_j}{\pi_j}\right)^{2} \]is non-negative whenever \(\pi_i\pi_j \ge \pi_{ij}\) for all pairs, a condition most sensible designs satisfy — whereas the first form's estimator can come out negative, which is an embarrassment in a report. And it exposes the design requirement directly: \(\pi_{ij}\) must be strictly positive for every pair, or the variance cannot be estimated at all.
Given. \(N = 4\) units with \(p = (0.1,\ 0.2,\ 0.3,\ 0.4)\) and \(Y = (3,\ 6,\ 9,\ 12)\), total \(Y = 30\). Two units are drawn successively without replacement, the probabilities being renormalised after the first draw.
Step 1 — the inclusion probabilities. For unit 1,
\[ \pi_1 = 0.1 + 0.2\!\left(\frac{0.1}{0.8}\right) + 0.3\!\left(\frac{0.1}{0.7}\right) + 0.4\!\left(\frac{0.1}{0.6}\right) = \frac{197}{840} = 0.234524, \]and in the same way
\[ \pi_2 = \frac{139}{315} = 0.441270, \qquad \pi_3 = \frac{73}{120} = 0.608333, \qquad \pi_4 = \frac{451}{630} = 0.715873. \]Step 2 — check the first identity.
\[ \sum_i \pi_i = 0.234524 + 0.441270 + 0.608333 + 0.715873 = 2 = n. \checkmark \]Note also that \(\pi_1 = 0.234524\) is not \(np_1 = 0.2\): successive sampling does not give inclusion probabilities exactly proportional to size.
Step 3 — the pairwise probabilities. For the pair \((1,2)\),
\[ \pi_{12} = \frac{(0.1)(0.2)}{0.9} + \frac{(0.2)(0.1)}{0.8} = \frac{17}{360} = 0.047222, \]and the remaining five are
\[ \pi_{13} = \frac{8}{105} = 0.076190, \quad \pi_{14} = \frac{1}{9} = 0.111111, \quad \pi_{23} = \frac{9}{56} = 0.160714, \] \[ \pi_{24} = \frac{7}{30} = 0.233333, \quad \pi_{34} = \frac{13}{35} = 0.371429. \]Step 4 — check the other two identities.
\[ \sum_{i<j}\pi_{ij} = 1 = \frac{n(n-1)}{2}, \qquad \sum_{j \ne 1}\pi_{1j} = \frac{17}{360} + \frac{8}{105} + \frac{1}{9} = \frac{197}{840} = \pi_1 = (n-1)\pi_1. \checkmark \]Step 5 — an estimate. If the sample is \(\{2, 4\}\),
\[ \hat Y_{HT} = \frac{6}{0.441270} + \frac{12}{0.715873} = 13.597122 + 16.762749 = 30.3599, \]against the true total \(30\).
Step 6 — the Yates–Grundy variance. The six terms use \(Y_i/\pi_i = 12.791878,\ 13.597122,\ 14.794521,\ 16.762749\):
\[ \operatorname{Var}\!\left(\hat Y_{HT}\right) = \sum_{i<j}\left(\pi_i\pi_j - \pi_{ij}\right) \left(\frac{Y_i}{\pi_i} - \frac{Y_j}{\pi_j}\right)^{2} = 2.428334. \]Step 7 — verify it by brute force. There are only six possible samples, and their probabilities are the \(\pi_{ij}\), which sum to 1. Computing \(E\!\left(\hat Y_{HT}\right)\) and \(E\!\left(\hat Y_{HT}^{2}\right) - \left[E\!\left(\hat Y_{HT}\right)\right]^{2}\) directly over those six gives
\[ E\!\left(\hat Y_{HT}\right) = 30.000000 = Y, \qquad \operatorname{Var}\!\left(\hat Y_{HT}\right) = 2.428334. \checkmark \]The Yates–Grundy formula and the definition agree to every digit, which is the point of doing the enumeration: the formula is not an approximation.
Interpretation. The four ratios \(Y_i/\pi_i\) run from \(12.79\) to \(16.76\) — not constant, so the variance is not zero, but close enough that a standard error of \(\sqrt{2.428334} = 1.558\) on a total of \(30\) is \(5.2\%\) of it. Had the size measure been chosen so that \(Y_i/\pi_i\) were constant, the design would have estimated the total with no error at all from two units out of four.