Sampling Weights in Significance Tests

A condensed rewrite of Tim Bock's "The Correct Treatment of Sampling Weights in Statistical Tests" (Displayr), with the derivation the original omits.

Part 1 — The problem, in brief

Sampling weights correct for over- and under-representation of groups in a survey. If a population is 51% female and a sample is 40% female, weights restore the balance.

The critical asymmetry is this:

Weights make your point estimates right. They do not make your standard errors right.

Percentages, means, regression coefficients, cluster centroids — essentially every parameter — are computed correctly by standard software once a weight is applied. What breaks is the standard error (SE), and therefore every p-value, confidence interval, and significance flag built on top of it.

Why it breaks

Most statistical routines were written for frequency weights (a weight of 3 means "this row occurred 3 times"), not sampling weights (a weight of 3 means "this respondent represents 3 people in the population"). Given sampling weights, the software computes the weighted sample size ${\sum }_i{w}_i$ and treats it as the number of independent observations.

That is simply false. If you have 67 female respondents and weight them up to represent 183 people, you still have 67 pieces of information. Testing as though you had 183 is inventing evidence.

The consequence is a systematic underestimate of the SE, and therefore an inflated false-positive rate. In the worked example in Part 3, the naive weighted SE is 26% too small, turning a genuinely null result into an apparently significant one.

What specific software does

The four common workarounds, and why each fails

1. Use the weighted percentages with the unweighted sample size. Works only in the special case where the test variable is the weighting variable. If you weight on gender and then compare men with women, the weights cancel within the comparison and you get away with it. Change the weight so it targets the outcome variable instead and the method fails immediately.

2. Run the significance test on unweighted data. Same objection, and a worse one: you are now testing a different quantity from the one you are reporting. In the Part 3 example, the unweighted test returns $p=0.020$ while the correct weighted test returns $p=0.161$. You would publish a "significant difference" that the weighted estimates do not support.

3. Rescale the weights to a mean of 1. This fixes only the most visible symptom. Weights summing to 1 make nothing significant; weights grossed to a population of 60 million make everything significant. Rescaling to a mean of 1 puts you in a plausible-looking middle ground — but the standard error is still computed from the wrong formula. Every calculation in Part 3's "naive" row uses mean-1 weights and is still wrong by a quarter.

4. Use Kish's effective sample size. $n_{\text{eff}}={\left({\sum }_i{w}_i\right)}^2/{\sum }_i{w}_i^2$, substituted for $n$ in the standard formulas.

This is the best of the four and is often close. But it is a special case, not a general solution — Section 2.9 derives exactly which assumption it makes and when that assumption fails. It is scale-invariant, which is a genuine virtue, and it is what to reach for when linearisation is intractable (correlations, for instance). It is not equivalent to doing the calculation properly.

The actual solution

Compute the standard error using Taylor series linearisation (TSL), also called the delta method or the linearisation method. Every serious survey package implements it. Bock's article stops here, declining to give the mathematics on the grounds that it would take a page and not help anyone.

It takes about three pages, and it does help — because once you can write down the influence values, you can implement any weighted test yourself, audit what your software is doing, and know when the assumptions have broken. That derivation is Part 2.

Part 2 — Taylor series linearisation, derived

2.1 Notation and setup

• $s$ — the achieved sample, of size $n$.
• $w_i$ — the final weight for respondent $i$ (design weight × non-response adjustment × calibration factor). Interpreted as the number of population units respondent $i$ represents.
• $y_i$ — the outcome. For a proportion, $y_i\in \{0,1\}$.
• $\hat{N}={\sum }_{i\in s}{w}_i$ — the estimated population size.

Two estimator families matter:

Horvitz–Thompson (HT) estimator of a population total, where HT refers to the standard design-unbiased estimator for sampling with unequal probabilities:

$${\hat{T}}_y=\sum _{i\in s}{w}_i{y}_i$$

Hájek estimator of a mean or proportion — the ratio of two estimated totals:

$$\hat{R}=\frac{{\sum }_i{w}_i{y}_i}{{\sum }_i{w}_i}=\frac{{\hat{T}}_y}{\hat{N}}$$

The whole problem lives in the gap between these two. ${\hat{T}}_y$ is linear in the data, so its variance can be written down exactly. $\hat{R}$ is a ratio of two random quantities — both numerator and denominator vary from sample to sample — and has no exact closed-form variance. Linearisation bridges the gap.

2.2 The one thing you can compute exactly: variance of a total

For a single-stage design treated as with-replacement (the standard convention, and conservative when the sampling fraction is small):

$$\boxed{\hat{V}({\hat{T}}_y)=\frac{n}{n-1}\sum _{i\in s}{\left(w_i{y}_i-\overline{wy}\right)}^2},\overline{wy}=\frac{1}{n}\sum _{i\in s}{w}_i{y}_i$$

This is just the sample variance of the weighted contributions $w_i{y}_i$, scaled up. Everything that follows reduces to this formula.

For a stratified, multi-stage design, use the ultimate cluster form: aggregate weighted values to the primary sampling unit (PSU) — the first unit selected, e.g. a postcode sector or an output area — then take variances between PSUs within stratum:

$$\hat{V}({\hat{T}}_y)=\sum _{h=1}^H\frac{n_h}{n_h-1}\sum _{i=1}^{n_h}{\left(z_{hi}-{\bar{z}}_h\right)}^2,z_{hi}=\sum _{j\in \text{PSU }hi}{w}_{hij}{y}_{hij}$$

with $H$ strata and $n_h$ PSUs in stratum $h$. Later stages of selection are ignored — under with-replacement PSU sampling this is exactly right, and otherwise slightly conservative.

2.3 The Taylor expansion

Let $\hat{\theta }=g({\hat{T}}_1,{\hat{T}}_2,\dots ,{\hat{T}}_k)$ be a smooth function of estimated totals. Expand $g$ about the true totals $\mathbf{T}$ and discard terms of second order and above:

$$\hat{\theta }-\theta \approx \sum _k{\frac{\partial g}{\partial T_k}|}_{\mathbf{T}}\left({\hat{T}}_k-T_k\right)$$

The right-hand side is a linear combination of totals. Its variance is computable by Section 2.2. Define the linearised variable (equivalently the influence value, or score) for respondent $i$:

$$u_i=\sum _k\frac{\partial g}{\partial T_k}{y}_{ki}$$

so that $\hat{\theta }-\theta \approx {\sum }_i{w}_i{u}_i-\text{const}$, and therefore:

$$\boxed{\hat{V}(\hat{\theta })=\hat{V}\left(\sum _{i\in s}{w}_i{u}_i\right)}$$

This is the entire method. Work out the influence value $u_i$ for your statistic, then feed $w_i{u}_i$ into the total-variance formula. The approximation is first-order, so it is asymptotically correct — accurate at realistic survey sample sizes, less reliable for very small bases (roughly $n<30$ per cell) or extremely skewed weights.

2.4 Applying it to a proportion or mean

For the Hájek ratio, $g(T_y,T_N)=T_y/T_N$. The partial derivatives are:

$$\frac{\partial g}{\partial T_y}=\frac{1}{T_N},\frac{\partial g}{\partial T_N}=-\frac{T_y}{T_N^2}=-\frac{R}{T_N}$$

Each respondent contributes $y_i$ to the numerator total and $1$ to the denominator total, so:

$$u_i=\frac{1}{T_N}\cdot y_i+\left(-\frac{R}{T_N}\right)\cdot 1=\frac{y_i-R}{T_N}$$

Substituting sample estimates, the influence value is the residual from the weighted mean, divided by the estimated population size:

$$u_i=\frac{y_i-\hat{R}}{\hat{N}}$$

Now apply Section 2.2. A useful simplification: because ${\sum }_i{w}_i(y_i-\hat{R})=0$ exactly, by the definition of $\hat{R}$, the centring term vanishes:

$$\boxed{\hat{V}(\hat{R})=\frac{n}{n-1}\cdot \frac{{\sum }_{i\in s}{w}_i^2{\left(y_i-\hat{R}\right)}^2}{{\left({\sum }_{i\in s}{w}_i\right)}^2}}$$

For a proportion, with $y_i\in \{0,1\}$ and $\hat{R}=\hat{p}$:

$$\hat{V}(\hat{p})=\frac{n}{n-1}\cdot \frac{{\sum }_i{w}_i^2(y_i-\hat{p}{)}^2}{{\hat{N}}^2}$$

This is implementable in three lines of any language, and it is what survey::svymean computes internally.

2.5 Two sanity checks

Equal weights must recover the textbook formula. Set $w_i=c$ for all $i$. Then $\hat{N}=nc$ and:

$$\hat{V}(\hat{R})=\frac{n}{n-1}\cdot \frac{c^2{\sum }_i(y_i-\bar{y}{)}^2}{n^2{c}^2}=\frac{{\sum }_i(y_i-\bar{y}{)}^2}{n(n-1)}=\frac{s^2}{n}$$

and for a binary outcome, $\hat{p}(1-\hat{p})/(n-1)$. Correct.

Scale invariance. Multiplying every weight by a constant $c$ scales the numerator by $c^2$ and the denominator by $c^2$. The variance is unchanged.

This second property is the formal reason the rescaling hacks are unnecessary. Done properly, it makes no difference whether your weights average 1, sum to 1, or gross up to 67 million — the p-value is identical. Any method whose answer changes when you rescale the weights is, by that fact alone, wrong.

2.6 Subgroups: the crosstab case

Almost every significance test in applied social research is a subgroup (domain) comparison. Let ${\delta }_i=1$ if respondent $i$ is in domain $d$, else $0$. Then:

$${\hat{p}}_d=\frac{{\sum }_i{w}_i{\delta }_i{y}_i}{{\sum }_i{w}_i{\delta }_i},{\hat{N}}_d=\sum _i{w}_i{\delta }_i$$

This is again a ratio of two totals, so the same derivation gives:

$$u_i=\frac{{\delta }_i\left(y_i-{\hat{p}}_d\right)}{{\hat{N}}_d}$$

Respondents outside the domain get $u_i=0$ — they contribute nothing to the estimate, but they still occupy their place in the design.

Practical warning. Compute $u_i$ over the full sample, then apply the variance formula across the full design. Do not filter the dataset to the subgroup and re-run. Under an unclustered design the two agree; under stratification or clustering they do not, because filtering discards strata and PSUs that legitimately contribute to the degrees of freedom, and can produce strata with a single remaining PSU. This is why survey uses subset() on the design object rather than on the data frame, and why filtering a survey dataset before analysis is a recurring source of wrong standard errors in government reporting.

2.7 Testing the difference between two subgroups

The instinct is to compute $\hat{V}({\hat{p}}_A)$ and $\hat{V}({\hat{p}}_B)$ separately and add them. Linearise the contrast instead. Let $\theta =p_A-p_B$; by linearity of the expansion:

$$\boxed{u_i=\frac{{\delta }_i^A\left(y_i-{\hat{p}}_A\right)}{{\hat{N}}_A}-\frac{{\delta }_i^B\left(y_i-{\hat{p}}_B\right)}{{\hat{N}}_B}}$$

Then $\hat{V}(\hat{\theta })=\hat{V}\left({\sum }_i{w}_i{u}_i\right)$ and the test statistic is:

$$t=\frac{{\hat{p}}_A-{\hat{p}}_B}{\sqrt{\hat{V}(\hat{\theta })}}$$

referred to a $t$ distribution on the design degrees of freedom (Section 2.8).

Why linearise the contrast rather than add variances? Because the contrast form automatically captures the covariance. For disjoint domains under an unclustered design, every respondent has $u_i$ from only one of the two terms, the cross-product is zero, and the covariance genuinely vanishes — adding the variances is fine. But as soon as there is clustering, both domains appear inside the same PSU, the PSU-level totals are correlated, and $\text{Cov}({\hat{p}}_A,{\hat{p}}_B)\ne 0$. Adding variances then overstates the standard error, typically making you less likely to detect real differences.

The same construction generalises to any linear contrast — differences of differences, wave-on-wave change in a panel, net promoter–style composites. Write down the contrast, differentiate, sum the influence values.

2.8 Degrees of freedom

Use a $t$ reference distribution, not the normal, with:

$$\text{df}=(\text{number of PSUs})-(\text{number of strata})$$

For an unclustered, unstratified sample this is $n-1$. For a clustered design it can be brutally small: a survey of 5,000 people across 120 sampling points stratified into 30 strata has $\text{df}=90$, not 4,999. Software that ignores this and uses the normal distribution will give confidence intervals that are too narrow.

For domain estimates, the conventional advice is to retain the full-design degrees of freedom rather than recomputing from the subgroup, unless the domain is confined to a small number of PSUs.

2.9 Where Kish's effective sample size comes from — and when it breaks

Kish's formula is not a rival to linearisation; it is a special case of it. Start from Section 2.4 and assume the squared residuals $(y_i-\hat{p}{)}^2$ are statistically independent of the weights, so that ${\sum }_i{w}_i^2(y_i-\hat{p}{)}^2\approx \frac{1}{n}\left({\sum }_i{w}_i^2\right)\left({\sum }_i(y_i-\hat{p}{)}^2\right)$. Writing ${\tilde{s}}^2=\frac{1}{n}{\sum }_i(y_i-\hat{p}{)}^2$:

$$\hat{V}(\hat{p})\approx \frac{n}{n-1}\cdot \frac{{\tilde{s}}^2{\sum }_i{w}_i^2}{{\left({\sum }_i{w}_i\right)}^2}$$

Divide by the simple random sampling (SRS) variance ${\tilde{s}}^2/(n-1)$ to get the design effect (deff) — the ratio of the true variance to what it would have been under simple random sampling:

$${\text{deff}}_{\text{Kish}}=\frac{n{\sum }_i{w}_i^2}{{\left({\sum }_i{w}_i\right)}^2},n_{\text{eff}}=\frac{n}{\text{deff}}=\frac{{\left({\sum }_i{w}_i\right)}^2}{{\sum }_i{w}_i^2}$$

So Kish's approximation is exact when the weights carry no information about the outcome, and approximate otherwise. Its error grows with the correlation between weight and outcome, and can run in either direction:

(600 simulated samples of n = 400 per row; see Part 3 for the setup.)

The uncomfortable point for applied work: weights are correlated with outcomes precisely when weighting was most necessary. If your weights are doing real work — correcting a serious non-response bias on a variable related to your outcome — that is exactly the situation in which Kish's approximation is least trustworthy. In a small worked case, I have seen the Kish design effect come out at 1.85 when the realised design effect was 0.48.

Kish also assumes a single design effect applies across the whole table, which is false whenever subgroups have different weight distributions.

2.10 Any smooth statistic, including regression

The recipe generalises. For a weighted linear or generalised linear model, the influence values are:

$$u_i={\left(\sum _j{w}_j{x}_j{x}_j^{\top }\right)}^{-1}{x}_i\left(y_i-x_i^{\top }\hat{\beta }\right)$$

— that is, the inverse information matrix times the score contribution. Applying the total-variance formula to $w_i{u}_i$ produces the sandwich estimator:

$$\hat{V}(\hat{\beta })={\left(X^{\top }WX\right)}^{-1}\left[\frac{n}{n-1}\sum _i{\left(w_i{x}_i{\hat{e}}_i-\overline{wx\hat{e}}\right)}^{\otimes 2}\right]{\left(X^{\top }WX\right)}^{-1}$$

which is what svyglm returns and why its standard errors differ from glm(weights = w). Design-based survey standard errors and heteroskedasticity-robust standard errors are the same idea applied to different sources of variability.

The pattern for any statistic: express it as a function of totals, differentiate, evaluate at the estimates, apply Section 2.2.

2.11 Chi-square tests: Rao–Scott corrections

The Pearson chi-square statistic assumes multinomial sampling and is invalid under a complex design. The Rao–Scott corrections rescale it using design effects computed by linearisation:

• First-order correction: divide the Pearson statistic by the mean of the generalised design effects, referring the result to the usual chi-square distribution. Simple, and slightly anti-conservative.
• Second-order (Satterthwaite) correction: additionally adjusts the degrees of freedom using the variability of the design effects. More accurate, and the default in Displayr, Q, and survey::svychisq.

This is the correct replacement for a chi-square test on a weighted crosstab. Note the corollary: there is no valid chi-square test you can perform on a weighted table alone. The correction needs the microdata, because design effects cannot be recovered from the cell counts. If a subcontractor sends you weighted tables and asks you to test them, the honest answer is that it cannot be done without the respondent-level file.

2.12 Two refinements that matter in government work

Calibrated weights are not fixed. If weights were calibrated (raked, post-stratified) to known population totals, the weights are themselves estimated, and treating them as fixed constants overstates the variance — often substantially, since calibration removes real sampling variability. The correct treatment replaces the residual $y_i-\hat{R}$ with the residual from a weighted regression of $y$ on the calibration variables. In R, use calibrate() or postStratify() on the design object so svyrecvar accounts for it; do not simply pass calibrated weights to svydesign and stop. Doing so is conservative rather than wrong, but you are discarding precision you paid for.

When linearisation is intractable, replicate. Some statistics have no tractable influence function — correlations, quantiles, medians, and most machine-learning outputs. Displayr's own documentation concedes it falls back to effective sample size in these cases. The better fallback is replication:

• Delete-one jackknife (JKn) — stratified, drops one PSU at a time and reweights.
• Balanced repeated replication (BRR) — for designs with exactly two PSUs per stratum; Fay's variant is more stable.
• Rao–Wu rescaled bootstrap — the most general, and the one to reach for by default.

Replication has a decisive practical advantage: the analyst re-estimates the statistic on each replicate and takes the variance across replicates, so any estimator can be handled, including a whole modelling pipeline. Many national statistical agencies now distribute replicate weight sets alongside microdata for precisely this reason. If your data provider offers them, use them — and re-run the calibration inside each replicate.

Part 3 — A worked example

A simulated survey of $n=400$ in two reporting subgroups (180 and 220), with lognormal calibration weights (coefficient of variation, the standard deviation divided by the mean, $=0.91$; maximum weight 8.0) and a modest correlation between weight and outcome — mimicking a design where hard-to-reach respondents both receive larger weights and answer differently.

True subgroup difference: small. Weighted estimates: 34.4% vs 43.8%, a gap of 9.4 percentage points.

Three things to take from this.

First, the two most common approaches disagree with the truth in opposite ways, and both would have been publishable. The unweighted test flags a significant difference. The naive weighted test on mean-1 weights sits on the wrong side of 0.05. The correct answer is a clear null.

Second, rescaling weights does not save you. Rows two and three use identical data and identical weights, differing only by a multiplicative constant. Their p-values differ by five orders of magnitude. The linearised answer would be unchanged by the same rescaling.

Third, the linearised standard error is right and the naive one is badly wrong. Repeating the whole exercise over 4,000 independent replications:

A standard error 26% too small inflates the type I error rate — the rate of false positives — from 5% to roughly 13%. On a crosstab with 200 significance tests, that is about 16 spurious "significant" findings instead of 10, every one of which is liable to end up in a summary of findings.

Reproducible code for the example (Python 3, numpy and scipy):

import numpy as np
from scipy import stats

def hajek(y, w):
    """Weighted (Hajek) proportion or mean."""
    return np.sum(w * y) / np.sum(w)

def lin_var_ratio(y, w):
    """Taylor-linearised variance of a weighted mean/proportion.
    Single-stage, with-replacement, unstratified design."""
    n = len(y)
    p = hajek(y, w)
    u = (y - p) / np.sum(w)        # influence value, Section 2.4
    z = w * u                      # mean(z) == 0 by construction
    return (n / (n - 1)) * np.sum(z**2)

def lin_var_contrast(y, w, dom_a, dom_b):
    """Variance of the difference between two domain proportions.
    dom_a, dom_b are boolean arrays over the FULL sample."""
    n = len(y)
    pa, pb = hajek(y[dom_a], w[dom_a]), hajek(y[dom_b], w[dom_b])
    Na, Nb = w[dom_a].sum(), w[dom_b].sum()
    u = dom_a * (y - pa) / Na - dom_b * (y - pb) / Nb   # Section 2.7
    z = w * u
    return (n / (n - 1)) * np.sum((z - z.mean())**2), pa - pb

# --- test ---
V, diff = lin_var_contrast(y, w, g == "A", g == "B")
se = np.sqrt(V)
df = n - 1                          # PSUs minus strata; Section 2.8
t  = diff / se
p  = 2 * stats.t.sf(abs(t), df)

Part 4 — Implementation

R, survey package (Lumley) — the reference implementation:

library(survey)

# Unclustered, weights only
des <- svydesign(ids = ~1, weights = ~w, data = dat)

# Realistic government design: stratified, clustered, calibrated
des <- svydesign(ids = ~psu, strata = ~stratum, weights = ~w,
                 data = dat, nest = TRUE)
des <- calibrate(des, ~ age_band + region + sex, pop.totals)  # see 2.12

svymean(~outcome, des)                      # linearised SE
svyby(~outcome, ~subgroup, des, svymean)    # subgroup estimates
svychisq(~outcome + subgroup, des,
         statistic = "F")                   # second-order Rao-Scott
svyglm(outcome ~ subgroup + age, des,
       family = quasibinomial())            # correct SEs; NOT glm(weights=)

# Direct contrast with covariance, per Section 2.7
svycontrast(svyby(~outcome, ~subgroup, des, svymean),
            quote(A - B))

# Replicate weights (Section 2.12)
rep_des <- as.svrepdesign(des, type = "bootstrap", replicates = 1000)

Use quasibinomial() rather than binomial() to suppress a spurious non-integer-successes warning; the point estimates are identical.

Stata — declare the design once with svyset, then prefix every command:

svyset psu [pweight = w], strata(stratum)
svy: proportion outcome, over(subgroup)
svy: tabulate subgroup outcome, row pearson
svy: logit outcome i.subgroup age

Note pweight (sampling weight), not fweight (frequency weight) or aweight (analytic weight). Choosing the wrong one silently reproduces the error in Part 1.

Python — statsmodels offers survey functionality but coverage is thinner than R's; the samplics package is more complete. For proportions and contrasts, the twenty lines in Part 3 are self-contained and auditable, which for regulated government outputs is arguably preferable to a dependency.

SPSS — use the Complex Samples module (CSPLAN then CSTABULATE / CSLOGISTIC). The base WEIGHT BY command does not do this.

Displayr / Q — applied automatically when a weight is set. No user action needed.

Part 5 — Practical checklist

For weighted survey outputs where significance is reported:

1. Never use WEIGHT BY in SPSS, or glm(weights=) in R, for anything with a p-value attached. These give correct percentages and wrong tests.
2. Declare the full design, not just the weight. Strata and PSUs change both the standard errors and the degrees of freedom, usually by more than the weights do.
3. Subset the design object, not the data frame. Filtering to a subgroup before analysis is a common and invisible error.
4. Check scale invariance as a smoke test. Multiply every weight by 1,000 and re-run. If any p-value moves, the method is wrong. This is a two-minute check that catches most misconfigurations, and it is worth building into an automated quality-assurance harness.
5. Report the design effect and effective sample size alongside the base. An unweighted base of 400 with a design effect of 2.0 is a base of 200 for inferential purposes, and readers should be told.
6. Watch small weighted subgroups. Linearisation is a first-order approximation and degrades at small bases; below roughly 30–50 respondents per cell, prefer descriptive reporting to significance flags regardless of method.
7. Specify the method in the technical annex. "Significance tests use Taylor series linearisation accounting for the stratified clustered design and calibration weights, with degrees of freedom equal to primary sampling units minus strata" is a sentence that survives peer review. "Tests at the 95% level" is not.
8. If you receive weighted tables without microdata, say you cannot test them. This is a defensible position, not an inconvenience — the design effects are not recoverable from the cells.

Sources and further reading

• Bock, T. The Correct Treatment of Sampling Weights in Statistical Tests, Displayr — the article condensed in Part 1.
• Lohr, S. (2019) Sampling: Design and Analysis, 3rd edn, Chapter 9 — the standard reference for linearisation, and the one Bock points to.
• Lumley, T. (2010) Complex Surveys: A Guide to Analysis Using R — the survey package explained by its author; the most directly useful text for practitioners.
• Wolter, K. (2007) Introduction to Variance Estimation, 2nd edn — comprehensive on both linearisation and replication.
• Rao, J.N.K. & Scott, A.J. (1984) "On chi-squared tests for multiway contingency tables with cell proportions estimated from survey data", Annals of Statistics 12(1) — the original chi-square corrections.
• Kish, L. (1965) Survey Sampling — the origin of the effective sample size approximation, including Kish's own caveats about it.