Survey Sampling and Weighting

Survey statistics provides the mathematical framework for drawing inferences about populations from samples selected through known, probabilistic mechanisms. Unlike model-based inference, design-based inference treats population values as fixed and randomness arises solely from the sampling process. This lesson develops the theory and practice of survey sampling, weighting, and variance estimation for complex survey designs.

Design-Based Inference

Horvitz-Thompson Estimator

Design Effects

Weighting Procedures

Survey weights adjust for unequal selection probabilities and non-response. The final analysis weight is typically a product of several components.

Raking (Iterative Proportional Fitting)

Variance Estimation for Complex Surveys

Complex survey designs violate the independence assumption underlying simple variance formulas. Specialized variance estimation methods account for the design.

Python Implementation

import numpy as np
import pandas as pd
from scipy import stats

np.random.seed(42)

# --- Simulate Complex Survey Data ---
N = 10000  # population size
n_strata = 4  # number of strata
n_psu_per_stratum = 25  # PSUs per stratum

# Population
strata = np.random.choice(n_strata, N, p=[0.3, 0.3, 0.25, 0.15])
y = 50 + 10 * strata + np.random.randn(N) * 15  # outcome correlated with stratum
x = 2 * strata + np.random.randn(N) * 5  # auxiliary variable

# Inverse probability sampling (higher prob in larger strata)
probs = np.where(strata < 2, 0.1, 0.2)  # oversample small strata
sample_idx = np.random.choice(N, size=1000, replace=False, p=probs / probs.sum())
design_weights = 1.0 / probs[sample_idx]

# --- Horvitz-Thompson Estimator ---
y_sample = y[sample_idx]
w = design_weights

y_ht = np.sum(w * y_sample) / np.sum(w)
y_srs = np.mean(y_sample)

print("=== Horvitz-Thompson Estimator ===")
print(f"Population mean (true): {np.mean(y):.2f}")
print(f"HT estimate: {y_ht:.2f}")
print(f"SRS mean: {y_srs:.2f}")

# --- Design Effect ---
var_srs = np.var(y_sample, ddof=1) / len(y_sample)
# Approximate variance with Taylor linearization
z_i = w * (y_sample - y_ht) / np.sum(w)
var_ht = np.sum(z_i**2) * len(y_sample) / (len(y_sample) - 1)
deff = var_ht / var_srs
n_eff = len(y_sample) / deff

print(f"\n=== Design Effects ===")
print(f"DEFF: {deff:.2f}")
print(f"Effective sample size: {n_eff:.0f}")

# --- Raking / Iterative Proportional Fitting ---
def rake_weights(weights, sample_data, population_margins, max_iter=50, tol=1e-6):
    """Rake weights to match population margins for multiple variables."""
    w = weights.copy()
    for iteration in range(max_iter):
        w_old = w.copy()
        for var_name, margin in population_margins.items():
            categories = sample_data[var_name]
            unique_cats = np.unique(categories)
            for cat in unique_cats:
                mask = categories == cat
                w[mask] *= margin[cat] / np.sum(w[mask])
        if np.max(np.abs(w - w_old)) < tol:
            print(f"  Raking converged after {iteration + 1} iterations")
            break
    return w

# Post-stratification variables
age_group = np.random.choice(['18-34', '35-54', '55+'], len(sample_idx), p=[0.3, 0.4, 0.3])
gender = np.random.choice(['M', 'F'], len(sample_idx), p=[0.48, 0.52])

# Known population margins
pop_margins = {
    'age': {'18-34': 0.28, '35-54': 0.38, '55+': 0.34},
    'gender': {'M': 0.49, 'F': 0.51}
}

sample_df = pd.DataFrame({'age': age_group, 'gender': gender})
w_raked = rake_weights(w, sample_df, pop_margins)
y_raked = np.sum(w_raked * y_sample) / np.sum(w_raked)

print(f"\n=== Raked Estimate ===")
print(f"Raked mean: {y_raked:.2f}")
print(f"Weight range: [{w_raked.min():.2f}, {w_raked.max():.2f}]")

# --- Jackknife Variance Estimation ---
# Simulate PSU structure
n_strata_sample = 4
psu_per_stratum = 25
total_sample = n_strata_sample * psu_per_stratum

# Create fake PSU assignments
psu_ids = np.repeat(np.arange(total_sample), len(y_sample) // total_sample + 1)[:len(y_sample)]
stratum_ids = np.repeat(np.arange(n_strata_sample), psu_per_stratum)[:len(y_sample)]

def theta_hat(y_vals, w_vals):
    return np.sum(w_vals * y_vals) / np.sum(w_vals)

# Delete-one-PSU jackknife
estimates = []
for h in range(n_strata_sample):
    for j in range(psu_per_stratum):
        mask = ~((stratum_ids == h) & (psu_ids == j))
        if mask.sum() > 10:
            theta_jk = theta_hat(y_sample[mask], w[mask])
            estimates.append(theta_jk)

estimates = np.array(estimates)
var_jack = np.var(estimates, ddof=1)
se_jack = np.sqrt(var_jack)

print(f"\n=== Jackknife Variance ===")
print(f"Estimate: {y_ht:.2f}")
print(f"SE (jackknife): {se_jack:.2f}")
print(f"95% CI: [{y_ht - 1.96*se_jack:.2f}, {y_ht + 1.96*se_jack:.2f}]")

Summary