Ratio, Regression and Calibration Estimation

In survey sampling theory, auxiliary information consists of data regarding one or more supplementary characteristics $\mathbf{x} = (x_1, x_2, \dots, x_p)^\top$ available for the elements of the population. In Units 3, 4, and 5, we investigated how auxiliary information is deployed at the design stage prior to data collection—using auxiliary variables to partition the universe into homogeneous strata, to delineate compact clusters, or to assign unequal selection probabilities proportional to size (PPS).

However, in many practical survey investigations, auxiliary information cannot be utilized during sample selection:

1. Absence of Unit-Level Frame Records: Detailed auxiliary values $x_i$ for every single individual unit in the frame may not be recorded or accessible when the sample is drawn. For example, a national statistical office may know the exact aggregate population counts of males and females across five-year age brackets from a recent national census, but the operational sampling frame of telephone numbers or residential addresses does not record the age or sex of the occupants.

2. Multi-Dimensional Auxiliary Vectors: When multiple auxiliary characteristics are available simultaneously (such as age, sex, educational attainment, and regional geography), constructing a multidimensional stratified or PPS design that accommodates all variables simultaneously leads to extreme stratum fragmentation, empty sampling cells, and operational paralysis.

3. Multi-Purpose Surveys: Large-scale national surveys gather data on hundreds of distinct substantive variables (e.g., employment, health expenditures, nutrition, energy consumption). An optimal design for one survey variable may be disastrously inefficient for another.

To overcome these structural constraints, modern survey theory leverages auxiliary information at the estimation stage, after the sample has been gathered. By benchmarking sample totals against known, reliable population totals derived from censuses or pristine administrative registers, the statistician can adjust the estimation formula or calibrate the survey weights. This post-survey adjustment achieves two goals:

• Substantial Variance Reduction: Correlated auxiliary variables absorb large fractions of sampling variability, producing estimators whose standard errors are significantly smaller than the unadjusted Horvitz–Thompson estimator.
• Consistent Benchmarking: Calibration guarantees that when the survey weights are applied to the auxiliary variables, the resulting survey estimates reproduce the external demographic benchmarks exactly, eliminating discrepancies between national statistical publications.

The ratio estimator is one of the oldest, most elegant, and most widely applied techniques in survey practice. It is particularly effective when the survey variable $y$ is approximately proportional to a single positive auxiliary variable $x$, such that the relationship between $y$ and $x$ passes approximately through the origin.

Formulation of the Estimator

Let $U = \{1, 2, \dots, N\}$ denote a finite population of $N$ units. Associated with each unit $i$ is the target variable $y_i$ and a strictly positive auxiliary variable $x_i > 0$. Suppose the population total of the auxiliary variable is known with certainty from an independent source:

A simple random sample without replacement (SRSWOR) of size $n$ is selected. Let $\bar{y} = \frac{1}{n} \sum_{i \in s} y_i$ and $\bar{x} = \frac{1}{n} \sum_{i \in s} x_i$ denote the unweighted sample means.

The finite-population ratio of the two characteristics is:

The sample estimate of the population ratio is:

The Ratio Estimator of the Population Mean $\bar{Y}$ is defined as:

Similarly, the Ratio Estimator of the Population Total $Y$ is:

Superpopulation Model Interpretation

Under the model-based superpopulation perspective (Brewer, 1963; Royall, 1970), the ratio estimator arises naturally as the best linear unbiased predictor under the following regression model:

where:

Under this model, the conditional expectation of $y$ is a straight line through the origin, and the residual variance is proportional to $x_i$ (heteroscedasticity). The weighted least-squares estimator of $\beta$ under this model is precisely $\hat{\beta} = \frac{\sum y_i}{\sum x_i} = \hat{R}$.

Finite-Sample Bias of the Ratio Estimator

Because the ratio estimator is a non-linear function of two random variables ($\bar{y}$ and $\bar{x}$), it is technically biased in finite samples. To quantify this bias, we deploy a Taylor series expansion.

Define the relative sample deviations:

Under SRSWOR, the expectations of these relative deviations are zero:

Their variances and covariance under the randomization distribution are:

where $C_y = S_y / \bar{Y}$ and $C_x = S_x / \bar{X}$ are the population coefficients of variation, $\rho = \frac{S_{xy}}{S_y S_x}$ is the finite-population Pearson correlation coefficient, and $S_{xy} = \frac{1}{N - 1} \sum_{i=1}^N (y_i - \bar{Y})(x_i - \bar{X})$ is the population covariance.

Now express the sample ratio $\hat{R}$ in terms of $e_0$ and $e_1$:

Assuming $|e_1| < 1$, we expand $(1 + e_1)^{-1}$ using a second-order Taylor series:

Multiplying through:

Taking expectations:

Substituting the expectations:

Therefore, the approximate finite-sample bias of the ratio estimator of the mean $\hat{\bar{Y}}_R = \hat{R} \bar{X}$ is:

Properties of the Bias:

1. Order of Magnitude $O(1/n)$: The bias is proportional to $1/n$. As sample size increases, the bias decays to zero much faster than the standard error (which decays at rate $O(1/\sqrt{n})$).

2. Zero Bias Condition: The bias is identically zero if $\rho = C_x / C_y$.

3. Practical Negligibility: In empirical practice with sample sizes $n \ge 30$ and coefficients of variation $C_x \le 0.10$, the ratio of the bias to the standard error ($\text{Bias} / \text{SE}$) is typically well below $0.10$. Cochran (1977) established that whenever $\text{Bias} / \text{SE} < 0.10$, the bias exerts a negligible effect on nominal confidence interval coverage rates.

Mean Squared Error and Asymptotic Variance

To derive the Mean Squared Error ($\text{MSE}$) of the ratio estimator, we use the first-order linear approximation:

Squaring both sides:

Taking the design expectation:

Substitute the variance and covariance expressions:

Translating coefficients of variation back to population variances:

This yields the fundamental Mean Squared Error Formula for the Ratio Estimator:

We can write this compactly in terms of the population residuals $E_i = y_i - R x_i$:

Therefore:

where $S_d^2 = \frac{1}{N - 1} \sum_{i=1}^N (y_i - R x_i)^2$ is the finite-population variance of the residuals from the line $y = Rx$.

Sample Variance Estimator

From sample data, the residual variance $S_d^2$ is estimated using the sample residuals $e_i = y_i - \hat{R} x_i$:

The standard approximately unbiased estimator of the variance of the ratio estimator is:

The Fundamental Efficiency Condition: When Does Ratio Estimation Beat Simple Random Sampling?

A practitioner must determine whether applying the ratio estimator improves precision over the unadjusted sample mean $\bar{y}$. The variance of the simple sample mean under SRSWOR is:

The ratio estimator is more efficient than the simple sample mean if and only if:

Subtracting $S_y^2$ from both sides:

Dividing by $R S_x > 0$:

Substituting $R = \bar{Y} / \bar{X}$:

This delivers the celebrated Efficiency Criterion for the Ratio Estimator:

Critical Implications of the Efficiency Threshold:

1. Positive Correlation is Insufficient: It is not enough for $\rho$ to be positive ($\rho > 0$). If $C_x$ is substantially larger than $C_y$ (the auxiliary variable is much more dispersed than the target variable), the threshold $\frac{C_x}{2 C_y}$ can be large (e.g., $0.60$ or $0.75$). If $\rho = 0.40$, applying the ratio estimator will increase the variance rather than decrease it!

2. Equal Variances: When the coefficients of variation are equal ($C_x \approx C_y$), the condition becomes simply $\rho > 0.50$. The correlation between $y$ and $x$ must exceed $0.50$ for the ratio estimator to outperform simple random sampling.

3. Diagnostic Rule: A survey statistician should never blindly apply ratio estimation without first verifying that the sample correlation satisfies $r > \frac{c_x}{2 c_y}$.

While the ratio estimator is highly effective when the relationship between $y$ and $x$ is linear and passes through the origin, it loses efficiency when the relationship has a non-zero intercept. The linear regression estimator relaxes this constraint, accommodating any linear relationship $y = \alpha + \beta x + \varepsilon$ with an arbitrary intercept.

The Difference Estimator as a Theoretical Precursor

To understand the regression estimator, we first examine the difference estimator. Suppose $k$ is a pre-determined mathematical constant chosen before observing the sample. The difference estimator of the population mean is:

Properties:

1. Strict Design Unbiasedness:

The difference estimator is strictly unbiased for any constant $k$.

2. Exact Sampling Variance:

3. Optimal Multiplier:

Differentiating $V(\bar{y}_d)$ with respect to $k$ and setting to zero:

where $B$ is the population linear regression slope of $y$ on $x$. Substituting $k = B$ into the variance formula:

Since $\rho^2 = \frac{S_{xy}^2}{S_y^2 S_x^2}$, we have:

The Linear Regression Estimator

In practice, the true population regression slope $B = S_{xy} / S_x^2$ is unknown. We estimate it from the sample using the ordinary least squares slope:

The Linear Regression Estimator of the population mean is:

The corresponding estimator of the population total is:

Finite-Sample Bias and Asymptotic Variance

Because $b$ is a stochastic random variable that covaries with $\bar{x}$, $\bar{y}_{lr}$ has a slight finite-sample bias:

Under bivariate normality, $\text{Cov}(b, \bar{x}) = 0$, making $\bar{y}_{lr}$ strictly unbiased. In general finite populations, the bias is of order $O(1/n)$ and vanishes rapidly as $n$ grows.

To first order of approximation, the variance of the linear regression estimator is equal to the optimal difference variance:

An approximately unbiased sample estimator of this variance is:

where $r = \frac{s_{xy}}{s_y s_x}$ is the sample correlation coefficient.

Universal Superiority of the Regression Estimator

Comparing the asymptotic variance of the regression estimator with the simple random sampling mean:

Because $\rho^2 \ge 0$, the regression estimator never increases variance relative to simple random sampling (to order $O(1/n)$). Even if the correlation is weak (e.g., $\rho = 0.20$), the regression estimator reduces variance by $1 - 0.20^2 = 4\%$. If the correlation is strong (e.g., $\rho = 0.90$), it reduces variance by $1 - 0.90^2 = 81\%$.

Now compare the regression estimator with the ratio estimator:

This yields an essential mathematical theorem:

The linear regression estimator is always more precise than or equal to the ratio estimator. The two estimators achieve identical precision if and only if $R = B$, which occurs when the population regression line passes strictly through the origin ($\alpha = \bar{Y} - B \bar{X} = 0$).

When the population is partitioned into $L$ strata and auxiliary totals $X_h$ are known for each stratum $h$, two distinct approaches exist for constructing ratio and regression estimators:

1. The Combined Estimator

In combined estimation, a single overall ratio or regression slope is computed across all strata:

where $\bar{y}_{st} = \sum_{h=1}^L W_h \bar{y}_h$ and $\bar{x}_{st} = \sum_{h=1}^L W_h \bar{x}_h$.

The bias of the combined estimator is of order $O(1/n)$, where $n = \sum n_h$ is the total sample size across the entire survey.

2. The Separate Estimator

In separate estimation, independent ratio or regression estimators are computed within each individual stratum $h$, and then aggregated:

The Bias vs. Precision Trade-off:

• Precision Gain: If the true regression slopes $B_h$ differ markedly across strata, the separate estimator captures these stratum-specific slopes, achieving a lower variance than the combined estimator:

• Cumulative Bias Hazard: The separate estimator incurs a bias in every single stratum of order $O(1/n_h)$. Summing across $L$ strata, the total bias can accumulate to $L \times O(1/n_h)$. If the within-stratum sample sizes $n_h$ are small (e.g., $n_h = 10$ or $15$), this cumulative bias can become severe.

Standard Survey Rule:

• If sample sizes within strata are large ($n_h \ge 30$ or $50$ for all $h$), use the separate estimator.
• If sample sizes within strata are small ($n_h < 20$) or if slopes $B_h$ are relatively homogeneous across strata, use the combined estimator to protect against cumulative bias.

In modern official statistics, surveys collect data on dozens of demographic, geographic, and economic variables simultaneously. The classical univariate ratio and regression estimators cannot readily accommodate multi-dimensional benchmark constraints across dozens of overlapping margins.

In 1992, Jean-Claude Deville and Carl-Erik Särndal formulated a universal optimization framework known as Calibration Weighting. Calibration transforms survey weighting into a constrained mathematical optimization problem.

Formal Formulation of the Calibration Problem

Let $s \subset U$ denote a probability sample drawn with design weights:

Associated with each sampled unit $k$ is an auxiliary vector of $p$ characteristics:

The population totals of these auxiliary variables are known from external censuses or administrative registers:

The objective of calibration is to find a set of new, calibrated weights $\{w_k : k \in s\}$ that:

1. Are as close as possible to the original design weights $d_k$, according to a pre-specified distance metric $G(w_k / d_k)$.
2. Satisfy the Calibration Equations (benchmarking constraints) exactly:

The Optimization Problem

Mathematically, calibration solves the constrained minimization problem:

where:

• $G(u)$ is a strictly convex distance function satisfying $G(1) = 0$, $G'(1) = 0$, and $G''(1) = 1$.
• $q_k > 0$ are known scaling factors (typically $q_k = 1$).

Classical Distance Functions

The choice of distance function $G(u)$ determines the operational behavior and range of the calibrated weights:

1. Linear (Chi-Square) Distance:

This leads directly to the Generalized Regression Estimator (GREG). It admits an exact closed-form analytical solution. However, linear calibration can occasionally produce negative weights ($w_k < 0$), which are physically impossible in survey estimation.

2. Exponential (Raking / Entropy) Distance:

This leads to Iterative Proportional Fitting (Multiplicative Raking). Because $u \log u$ is defined only for $u > 0$, the calibrated weights are guaranteed to be strictly positive ($w_k > 0$ for all $k$).

3. Logit (Bounded) Distance:

This restricts the weight adjustment factor to a pre-specified bounded interval:

Bounded calibration prevents extreme weight inflation, protecting estimators against excessive variance while maintaining exact benchmarking.

The Generalized Regression Estimator (GREG)

Under the linear distance function $G(u) = \frac{1}{2}(u - 1)^2$, the optimization problem is solved using Lagrange multipliers. The Lagrangian is:

Differentiating with respect to $w_k$ and setting to zero:

Substitute $w_k$ into the calibration equations:

Let $\hat{\mathbf{X}}_\pi = \sum_{k \in s} d_k \mathbf{x}_k$ denote the uncalibrated Horvitz–Thompson estimator of the auxiliary totals. Then:

Assuming the sample covariance matrix $\hat{\mathbf{T}}_{xx} = \sum_{k \in s} d_k \mathbf{x}_k \mathbf{x}_k^\top$ is invertible:

Substituting $\boldsymbol{\lambda}$ back into the weight formula yields the Exact GREG Weight Formula:

The calibrated GREG estimator of the population total of any survey variable $y$ is:

where $\hat{\mathbf{B}} = \left( \sum_{j \in s} d_j \mathbf{x}_j \mathbf{x}_j^\top \right)^{-1} \sum_{j \in s} d_j \mathbf{x}_j y_j$ is the weighted least-squares regression coefficient vector.

Asymptotic Equivalence Theorem

Deville and Särndal (1992) proved a remarkable theoretical result: For any valid distance function $G(u)$, the resulting calibrated estimator $\hat{Y}_{cal} = \sum_{k \in s} w_k y_k$ is asymptotically equivalent to the Generalized Regression Estimator (GREG) to first order of approximation:

Consequently, all calibration estimators (whether linear, raking, or logit) share the exact same asymptotic variance:

where $E_i = y_i - \mathbf{x}_i^\top \mathbf{B}$ are the population regression residuals.

In modern survey field operations, nonresponse is inevitable. When field operations encounter non-respondents, simply replacing non-respondents with adjacent households introduces fatal selection bias. In 1946, Morris Hansen and William Hurwitz formulated an ingenious application of Two-Phase (Double) Sampling for Nonresponse, which converts nonresponse from an uncontrolled bias hazard into a controlled probability design.

The Two-Phase Nonresponse Protocol

Consider a target population of $N$ units. The population conceptually consists of two strata:

1. Response Stratum ($U_1$): Units who would respond if contacted under the initial standard survey mode (e.g. mail questionnaire or online web survey). Size $N_1$, proportion $W_1 = N_1 / N$.

2. Nonresponse Stratum ($U_2$): Units who do not respond to the initial survey mode. Size $N_2$, proportion $W_2 = N_2 / N$.

Neither $N_1$ nor $N_2$ is known prior to data collection. The survey protocol proceeds in two sequential phases:

1. Phase 1: An initial probability sample of $n$ units is selected from $N$ via SRSWOR. The questionnaire is administered via a low-cost survey mode (e.g., electronic mail or self-administered web portal at unit cost $c_1$).

• $n_1$ units respond successfully, providing observations $y_{11}, \dots, y_{1n_1}$.
• $n_2 = n - n_1$ units fail to respond.

2. Phase 2 (Subsampling Nonrespondents): From the $n_2$ non-respondents, a subsample of size $m = n_2 / k$ ($k \ge 1$) is selected at random via SRSWOR.

• Field interviewers deploy an intensive, high-cost tracking protocol (e.g., personal face-to-face visits with monetary incentives at unit cost $c_2 \gg c_1$) to achieve a $100\%$ response rate among the $m$ selected nonrespondents, obtaining measurements $y_{21}, \dots, y_{2m}$.

Hansen–Hurwitz Unbiased Estimator

The Hansen–Hurwitz Estimator of the population mean $\bar{Y}$ is:

where:

• $\bar{y}_1 = \frac{1}{n_1} \sum_{i=1}^{n_1} y_{1i}$ is the sample mean of the initial respondents.
• $\bar{y}_{2m} = \frac{1}{m} \sum_{j=1}^m y_{2j}$ is the sample mean of the subsampled nonrespondents.
• $w_1 = n_1 / n$ and $w_2 = n_2 / n$ are the empirical sample response and nonresponse proportions.

Proof of Strict Unbiasedness:

Using the law of total expectation, condition on the realized Phase 1 sample $s$:

In Phase 2, the $m$ units are drawn via SRSWOR from the $n_2$ nonrespondents. Therefore, $E_2(\bar{y}_{2m} \mid s) = \bar{y}_2$, where $\bar{y}_2 = \frac{1}{n_2} \sum_{i \in s_2} y_{2i}$ is the mean of all nonrespondents in the initial sample $s$:

where $\bar{y}$ is the unobserved overall sample mean of the full sample of size $n$. Now take expectation over Phase 1 SRSWOR selection from the population $U$:

Thus:

The Hansen–Hurwitz estimator is strictly design-unbiased for the true population mean, completely eliminating nonresponse bias through probability subsampling.

Exact Sampling Variance and Optimal Allocation

Using the conditional variance decomposition $V(\hat{\bar{Y}}_{HH}) = V_1(E_2(\hat{\bar{Y}}_{HH} \mid s)) + E_1(V_2(\hat{\bar{Y}}_{HH} \mid s))$:

1. $V_1(E_2(\hat{\bar{Y}}_{HH} \mid s)) = V_1(\bar{y}) = \left( \frac{1 - f}{n} \right) S^2$
2. $E_1(V_2(\hat{\bar{Y}}_{HH} \mid s)) = E_1 \left[ \left(\frac{n_2}{n}\right)^2 \left(\frac{1 - m/n_2}{m}\right) s_{2}^2 \right] = \left( \frac{k - 1}{n} \right) W_2 S_2^2$

Summing both components yields the exact variance:

where:

• $S^2$ is the overall population variance.
• $S_2^2$ is the variance among units in the nonresponse stratum $U_2$.
• $W_2$ is the population nonresponse rate.

Cost Optimization:

The total expected survey budget is:

where $c_0$ is overhead fixed cost, $c_1$ is initial mail/web contact cost, and $c_2$ is intensive face-to-face field follow-up cost.

Minimizing the product of variable cost and variance via Cauchy–Schwarz inequality yields the Optimal Subsampling Interval ($k^*$) and Optimal Initial Sample Size ($n^*$):

When follow-up costs are high ($c_2 / c_1 \gg 1$), $k^*$ increases, directing field interviewers to intensively trace a smaller, highly targeted subsample of nonrespondents while preserving design unbiasedness and protecting survey budgets.

In our derivation of the classical ratio estimator $\hat{\bar{Y}}_R = \frac{\bar{y}}{\bar{x}} \bar{X}$, we proved that the estimator possesses a finite-sample bias of order $O(1/n)$:

While this bias is negligible when $n \ge 50$ and $C_x \le 0.10$, in small-scale inquiries, environmental monitoring plots ($n \le 15$), or highly fragmented sub-strata, the ratio bias can distort point estimates and degrade nominal confidence interval coverage. Survey theory provides two celebrated mathematical techniques to eliminate or dramatically reduce this bias: the Quenouille Jackknife and the Hartley–Ross exactly unbiased estimator.

1. Quenouille's Jackknife Bias-Reduced Ratio Estimator

In 1956, Maurice Quenouille introduced the Jackknife as a universal nonparametric device to eliminate first-order bias from non-linear estimators.

Let a sample of $n$ pairs $(x_i, y_i)$ be partitioned randomly into $g$ mutually exclusive groups of size $m = n / g$ ($g \ge 2$). Let $\bar{y}_{(-j)}$ and $\bar{x}_{(-j)}$ denote the sample means of the $n - m$ observations remaining after deleting the $j$-th group ($j = 1, 2, \dots, g$). The partial ratio estimate computed on the remaining $g - 1$ groups is:

The average of the $g$ partial ratio estimates is:

The full-sample ratio estimate is $\hat{R} = \bar{y} / \bar{x}$. Quenouille's Jackknife Ratio Estimator is defined as:

The corresponding estimator of the population mean is $\hat{\bar{Y}}_Q = \hat{R}_Q \bar{X}$.

Proof of Bias Reduction:

Recall that the bias of the full-sample ratio estimator of size $n$ has the asymptotic expansion:

where $a_1$ and $a_2$ are population constants independent of sample size.

Each partial ratio estimate $\hat{R}_{(-j)}$ is based on a sample of size $n - m = n(1 - 1/g) = n \left(\frac{g - 1}{g}\right)$. Therefore:

Now take the expectation of Quenouille's linear combination:

Notice what happens to the first-order bias term:

The leading $O(1/n)$ bias term cancels out identically!

The remaining expectation is:

Therefore:

Quenouille's Jackknife reduces the bias from order $O(1/n)$ down to order $O(1/n^2)$. In a sample of $n = 20$, the bias is reduced by a factor of approximately 20, rendering it completely negligible.

2. The Hartley–Ross Exactly Unbiased Ratio Estimator

In 1954, H. O. Hartley and A. Goodman Ross derived a mathematically strictly design-unbiased ratio-type estimator of the population mean $\bar{Y}$.

Let $r_i = \frac{y_i}{x_i}$ denote the individual ratio for each sampled unit ($i = 1, \dots, n$). Let $\bar{r} = \frac{1}{n} \sum_{i \in s} r_i$ denote the unweighted average of the individual ratios. Let $s_{rx} = \frac{1}{n - 1} \sum_{i \in s} (r_i - \bar{r})(x_i - \bar{x})$ denote the sample covariance between the individual ratios and the auxiliary values $x_i$.

The Hartley–Ross Estimator of the population mean $\bar{Y}$ is:

Proof of Strict Design Unbiasedness:

For each unit in the population, $y_i = r_i x_i$. The population mean of $y$ is:

Using the algebraic covariance identity:

where $\bar{R}_u = \frac{1}{N} \sum_{i=1}^N r_i$ is the population mean of individual ratios, and $S_{rx} = \frac{1}{N - 1} \sum_{i=1}^N (r_i - \bar{R}_u)(x_i - \bar{X})$ is the population covariance.

Under SRSWOR:

Taking expectations of the Hartley–Ross estimator:

This proves that $\hat{\bar{Y}}_{HR}$ is strictly design-unbiased for any finite sample size $n \ge 2$, with identically zero bias.

When multiple auxiliary characteristics $\mathbf{x} = (x_1, x_2, \dots, x_p)^\top$ are available with known population totals $\mathbf{X} = (X_1, X_2, \dots, X_p)^\top$, survey precision can be dramatically enhanced by exploiting all auxiliary variables simultaneously.

1. Olkin's Multivariate Ratio Estimator

In 1958, Ingram Olkin extended the ratio estimator to $p$ auxiliary variables. For each auxiliary variable $j \in \{1, 2, \dots, p\}$, compute the univariate ratio estimate of the mean:

Olkin's Multivariate Ratio Estimator is a weighted convex combination:

The sampling variance of $\hat{\bar{Y}}_{MR}$ is:

where $\boldsymbol{\Sigma}_R$ is the $p \times p$ covariance matrix among the individual ratio estimators:

Applying Lagrange multipliers to minimize $\mathbf{W}^\top \boldsymbol{\Sigma}_R \mathbf{W}$ subject to $\mathbf{1}^\top \mathbf{W} = 1$ yields the optimal weight vector:

Olkin's multivariate ratio estimator extracts the joint predictive power of all $p$ auxiliary variables, achieving lower variance than any single univariate ratio estimator.

2. The Multiple Linear Regression Estimator

When the relationship between $y$ and $\mathbf{x}$ involves an intercept, the Multiple Linear Regression Estimator provides optimal precision.

Let $\mathbf{x}_i = (x_{i1}, \dots, x_{ip})^\top$ denote the vector of auxiliary characteristics for unit $i$, and let $\bar{\mathbf{X}} = \frac{1}{N} \sum_{i=1}^N \mathbf{x}_i$ denote the vector of known population means. From the sample of size $n$, compute the sample means $\bar{y}$ and $\bar{\mathbf{x}}$, and compute the sample OLS slope vector:

The Multiple Linear Regression Estimator of the population mean is:

To first order of approximation, the asymptotic variance of the multiple regression estimator is:

where $R_{y \cdot x_1 \dots x_p}^2$ is the Multiple Coefficient of Determination (the proportion of total variance in $y$ explained by the linear combination of the $p$ auxiliary variables).

Notice the immense inferential power of this result: If an investigator includes three auxiliary variables that together explain $85\%$ of the variation in $y$ ($R^2 = 0.85$), the multiple regression estimator reduces the sampling variance by an astounding $85\%$, cutting the required sample size by a factor of nearly 7 compared to simple random sampling!

In Section 6.5, we introduced the exponential distance function $G(u) = u \log u - u + 1$, noting that it generates Iterative Proportional Fitting (IPF), commonly called Multiplicative Raking. Raking was invented by W. Edwards Deming and Frederick F. Stephan in 1940 to adjust cross-tabulations of census data so that marginal totals match known external census benchmarks.

Information-Theoretic Foundations

The mathematical elegance of raking lies in its connection to information theory. Let $d_k$ denote the initial design weights, and let $w_k$ denote the calibrated weights. The Kullback–Leibler Information Divergence (relative entropy) from the design distribution to the calibrated distribution is defined as:

Minimizing relative entropy subject to the calibration constraints $\sum_{k \in s} w_k \mathbf{x}_k = \mathbf{X}$ yields the Lagrangian:

Differentiating with respect to $w_k$:

Solving for $w_k$:

Why Raking Prevents Negative Weights:

Because the exponential function $\exp(z)$ is strictly positive for all real numbers $z \in (-\infty, \infty)$:

Raking weights are guaranteed to be strictly positive! The pathology of negative weights produced by linear GREG is mathematically impossible under raking.

The Deming–Stephan Raking Algorithm for Two-Dimensional Margins

Consider a survey sample cross-classified by two categorical auxiliary variables:

• Row variable $A$ with $R$ categories (e.g., age groups).
• Column variable $B$ with $C$ categories (e.g., educational attainment).

Let $N_{r\cdot}$ denote the known population total for row category $r$ ($r = 1, \dots, R$). Let $N_{\cdot c}$ denote the known population total for column category $c$ ($c = 1, \dots, C$). Let $w_{rc}^{(0)} = \sum_{i \in s_{rc}} d_i$ denote the initial design-weighted total in cell $(r, c)$.

The iterative algorithm cycles between rows and columns:

1. Row Step: Adjust weights so that row totals match the row benchmarks:

2. Column Step: Adjust weights so that column totals match the column benchmarks:

3. Iteration: Repeat steps 1 and 2 sequentially for $t = 0, 1, 2, \dots$ until the maximum discrepancy between weighted totals and marginal benchmarks across all rows and columns is less than a tolerance $\epsilon = 10^{-6}$.

Deming and Stephan (1940) and Ireland and Kullback (1968) proved that if the contingency table cells are strictly positive, the raking algorithm converges geometrically to the unique global minimum of the relative entropy function.

In classical GREG estimation, the assistance model is strictly linear: $y_i = \mathbf{x}_i^\top \boldsymbol{\beta} + \varepsilon_i$. In environmental, ecological, and complex socioeconomic phenomena, the relationship between $y$ and $\mathbf{x}$ is frequently non-linear (e.g., agricultural crop yields saturating with fertilizer application, or energy consumption exhibiting quadratic temperature response).

When the linear model is misspecified, GREG remains design-consistent, but it loses substantial efficiency. To overcome this limitation, F. Jay Breidt and Jean Opsomer (2000) formulated Nonparametric Model-Assisted Calibration using local polynomial smoothing.

Local Polynomial Calibration

Assume the superpopulation model:

where $m(\cdot)$ is an unknown smooth, non-linear function, and $E_\xi(\varepsilon_i \mid x_i) = 0$.

For each sampled unit $k \in s$, the non-linear curve $m(x_k)$ is estimated via kernel-weighted local polynomial regression:

where $\mathbf{W}_{k,s} = \text{diag}\left\{ d_i K\left( \frac{x_i - x_k}{h} \right) \right\}$ is a matrix of kernel weights governed by bandwidth $h$.

The nonparametric regression estimator of the population total is:

Breidt and Opsomer proved that $\hat{Y}_{np}$ is asymptotically design-unbiased for any underlying population structure, while achieving the theoretical lower bound on sampling variance when the true relationship is non-linear.

In Section 6.5, we derived the Generalized Regression Estimator (GREG) under single-stage sampling with arbitrary inclusion probabilities. In large-scale national surveys, the sampling design is virtually always a stratified multistage clustered design. Establishing how GREG behaves in the presence of hierarchical clustering, stratification, and multi-stage inclusion weights is essential for modern survey practice.

The Hierarchical GREG Formulation

Consider a population partitioned into $H$ strata. In stratum $h$, there are $N_h$ Primary Sampling Units (PSUs), and PSU $(h, i)$ contains $M_{hi}$ secondary units. Let $w_{hij}$ denote the base design weight of observation unit $j$ within PSU $i$ of stratum $h$:

where $\pi_{hi}$ is the first-stage inclusion probability of the PSU and $\pi_{j|hi}$ is the second-stage conditional inclusion probability.

Associated with each observation unit is a $p$-dimensional vector of auxiliary benchmark variables:

The true population total vector is known from external censuses or administrative registers:

Let $\hat{\mathbf{X}}_\pi$ denote the uncalibrated Horvitz–Thompson estimator of the auxiliary totals:

We specify a linear assisting model:

where $q_{hij} > 0$ are known heteroscedasticity weights.

The design-weighted generalized least-squares estimator of the regression parameter vector $\boldsymbol{\beta}$ is:

The Complex GREG Estimator of the Population Total is:

where $\hat{Y}_\pi = \sum_{h,i,j} w_{hij} y_{hij}$.

The Calibrated $g$-Weights

A remarkable feature of the GREG estimator is that it can be expressed directly as a weighted sum of the observed survey values:

where $w_{hij}^*$ is the calibrated analytic weight:

and the calibration adjustment factor (or $g$-factor) for unit $(h, i, j)$ is:

Fundamental Invariance Property:

When the calibrated weights $w_{hij}^*$ are applied to the auxiliary vector $\mathbf{x}_{hij}$, they reproduce the known benchmark totals exactly:

This mathematical guarantee holds for any complex multistage design, ensuring that published survey tables match external demographic census totals across all benchmark margins.

Design-Consistent Variance Estimation Using Linearized Residuals

To evaluate the sampling variance of $\hat{Y}_{GREG}$ under the complex multistage design, we compute the sample regression residuals:

Under the Ultimate Cluster principle, the design variance of $\hat{Y}_{GREG}$ is computed directly by aggregating the weighted residuals at the Primary Sampling Unit level:

The mean weighted residual per PSU in stratum $h$ is:

The design-consistent $g$-Weighted Residual Variance Estimator is:

Notice the inclusion of the $g$-factor inside $e_{hi}^* = \sum w_{hij} g_{hij} e_{hij}$: Särndal, Swensson, and Wretman (1992) proved that incorporating $g_{hij}$ into the residual variance estimator accounts for the conditional variance of the estimator given the realized auxiliary sample mean, delivering superior confidence interval coverage in moderately sized samples.

A recurring operational hazard in complex survey calibration occurs when the auxiliary vector $\mathbf{x}$ contains a large number of categorical indicator variables crossing multiple demographic dimensions (e.g., age, sex, education, marital status, race, and geographic district).

When some demographic cross-classification cells in the sample are sparsely populated or collinear, the sample cross-product matrix:

becomes ill-conditioned (its condition number explodes, and its determinant approaches zero).

Under standard linear GREG calibration, inverting an ill-conditioned matrix $\hat{\mathbf{T}}_{xx}$ causes the Lagrange multipliers $\boldsymbol{\lambda}$ to take extreme values. Consequently:

• Many calibrated weights $w_k^*$ become negative ($w_k^* < 0$), which violates physical reality.
• Other weights become massive outliers (e.g., $w_k^* / d_k > 50$), inflating sampling variance and destabilizing domain estimates.

The Ridge Calibration Formulation

To resolve ill-conditioning without dropping essential demographic benchmarks, P. Bardsley and R. L. Chambers (1984) and Raymond Chambers (1996) formulated Ridge Calibration.

In ridge calibration, the statistician recognizes that satisfying calibration equations with absolute mathematical perfection on noisy sample data is sub-optimal. Instead, the calibration problem is relaxed into a penalized optimization problem that balances two competing objectives:

1. Minimizing the distance between the calibrated weights $w_k$ and the design weights $d_k$.
2. Minimizing a weighted quadratic penalty on the calibration discrepancy $\mathbf{C} = \sum_{k \in s} w_k \mathbf{x}_k - \mathbf{X}$.

Mathematically, ridge calibration solves:

where:

• $\lambda_{ridge} > 0$ is a scalar ridge tuning parameter governing the penalty weight.
• $\mathbf{C}$ is a diagonal cost matrix specifying the relative importance of matching specific benchmarks (e.g., total population count must match tightly, while a detailed age-by-education interaction can tolerate a minor discrepancy).

The Closed-Form Ridge Weight Formula

Differentiating the penalized objective function with respect to $w_k$ and setting to zero yields the exact Ridge Calibration Weight Formula:

where:

• $\mathbf{d} = (d_1, \dots, d_n)^\top$ is the vector of initial design weights, and $\mathbf{D} = \text{diag}(\mathbf{d})$.
• $\mathbf{Q} = \text{diag}(q_1, \dots, q_n)$ is the diagonal matrix of model variance weights.
• $\mathbf{X}_s$ is the $p \times n$ matrix of sample auxiliary values.

The Magic of the Ridge Regularizer:

Notice what happens algebraically inside the inverted matrix:

The addition of the diagonal ridge term $\lambda_{ridge} \mathbf{C}$ guarantees that the eigenvalues of the matrix are bounded strictly away from zero!

• Even if $\mathbf{X}_s \mathbf{D} \mathbf{Q} \mathbf{X}_s^\top$ is strictly singular or severely ill-conditioned, the regularized matrix is strictly positive definite and easily invertible.
• The tuning parameter $\lambda_{ridge}$ pulls extreme weights back toward the initial design weights $d_k$, completely eliminating negative weights and extreme leverage outliers.
• By accepting a tiny, negligible calibration discrepancy on collinear margins, ridge calibration achieves a massive reduction in sampling variance, producing stable, robust weights across all survey indicators.

In agricultural economics, consumer price monitoring, and longitudinal health research, the parameter of interest is frequently not a simple ratio $R = Y / X$, but a ratio of two ratios—known as the Double Ratio Estimator.

The Double Ratio Architecture

Consider a two-wave longitudinal inquiry or an economic index comparing two market sectors. Let:

• $\hat{R}_1 = \bar{y}_1 / \bar{x}_1$ denote the estimated ratio in Sector 1 (e.g., crop yield per hectare in Province 1).
• $\hat{R}_2 = \bar{y}_2 / \bar{x}_2$ denote the estimated ratio in Sector 2 (e.g., crop yield per hectare in Province 2).

The target parameter is the Relative Comparative Efficiency:

The sample estimator is:

Taylor Series Linearization of the Double Ratio

To derive the sampling variance of $\hat{\theta}$, we express $\log \hat{\theta}$ in terms of the four sample means:

Taking differential variations:

Squaring both sides and taking design expectations yields the asymptotic relative variance:

Substituting the Taylor expansions for the individual ratios:

Methodological Interpretation:

When Sector 1 and Sector 2 share common macro-environmental conditions, the covariance between the two ratios $\text{Cov}(\hat{R}_1, \hat{R}_2)$ is strongly positive. The subtraction of this positive covariance in the double ratio variance formula results in substantial variance cancellation, allowing economists to measure relative efficiency shifts and price elasticities with extraordinary empirical precision.

In modern survey software suites (such as R's `survey` package or Python's `samplics`), calibration weighting is the standard production step connecting raw field data to published official estimates.

Best-Practice Calibration Workflow:

1. Formulate Auxiliary Vector $\mathbf{x}_k$: Include an intercept term $1$ to preserve the population size $N$, followed by indicator dummy variables for categorical margins (age categories, sex, educational levels, geographic districts). To prevent collinearity, drop one redundant category per marginal partition.

2. Execute Linear Calibration (GREG): Compute initial weights. Inspect the resulting weight distribution:

• Check the minimum weight $\min w_k$. If any $w_k < 0$, linear calibration is discarded.
• Check the ratio of maximum to minimum weight $\max w_k / \min w_k$. If the ratio exceeds $10.0$ or $15.0$, variance inflation will occur.

3. Switch to Logit or Raking Calibration if Needed: If linear weights are negative or extreme, rerun calibration using the Logit distance function with bounds $[L, U] = [0.4, 2.5]$, ensuring that no weight is reduced below $40\%$ or inflated above $250\%$ of its design weight.

4. Variance Estimation via Linearized Residuals: To estimate standard errors of calibrated estimators, compute the survey residuals $e_k = y_k - \mathbf{x}_k^\top \hat{\mathbf{B}}$. Apply the complex survey design variance formula (incorporating strata and clusters) directly to the weighted residuals $w_k e_k$, guaranteeing honest standard errors.

In the next unit, we address what happens when sampled units cannot be observed—examining Survey Nonresponse, Imputation, and the Total Survey Error paradigm.