The Confidence Interval in Your Tables Is Probably Wrong

Every survey table with a small base carries a confidence interval, and in most software that interval is computed by the formula everyone learns first:

$$\hat{p}\pm z_{0.975}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$$

This is the Wald interval. It is the one in the textbooks, the one in the significance testing appendix of most tabulation packages, and the one behind the "±3%" that gets quoted in the press.

It is also, for small bases and for proportions near 0 or 1, badly and systematically wrong — not slightly conservative, not approximately right, but capable of delivering 33% coverage while claiming 95%. This has been documented since at least 1998 and the fix predates it by seventy years. It is still the default nearly everywhere.

Part 1 — Exact coverage, not a simulation

Coverage of a binomial interval can be computed exactly, without simulation. For a given $n$ and true $p$, enumerate every possible outcome $x=0,1,\dots ,n$, construct the interval for each, and sum the binomial probabilities of the outcomes whose interval contains $p$:

$$C(p,n)=\sum _{x=0}^n\left(\genfrac{}{}{0ex}{}{n}{x}\right)p^x(1-p{)}^{n-x}\cdot 1\left[p\in \left(L(x,n),U(x,n)\right)\right]$$

No Monte Carlo error, no random seed. The numbers below are the actual coverage probabilities.

Four methods, nominal 95%:

Three things follow.

The failure is severe, not marginal. At $n=20$, $p=0.02$ the Wald interval covers the truth a third of the time. An analyst reporting "2%, ±95% CI" on that base is making a statement that is false two times in three.

The failure persists at bases nobody thinks of as small. At $n=100$ and a 5% incidence — an entirely routine cell in a survey table — Wald coverage is 87.8%. At $n=400$, still 92.7%. The interval is too narrow across the whole range of low-incidence reporting.

It is not only about small $n$. At $n=20$, $p=0.5$, Wald achieves 95.9% — essentially correct. The problem is not the base alone but the combination of base and distance from 0.5. Low-incidence items on moderate bases are the danger zone, and low-incidence items are exactly what surveys are often commissioned to measure.

The degenerate case

There is a failure mode worse than under-coverage. When $x=0$, the Wald interval is:

$$0\pm z\sqrt{\frac{0\cdot 1}{n}}=[0,0]$$

A zero-width interval. The software reports, with 95% confidence, that the population proportion is exactly zero — on the basis of having seen nobody. The same happens at $x=n$, where it asserts exactly 100%.

At $n=30$ and a true rate of 2%, $x=0$ occurs about 55% of the time. So the single most likely output of that cell is an interval claiming certainty about a quantity the data cannot pin down at all.

And when $x$ is small but non-zero, the interval routinely runs outside $[0,1]$. Probability that a Wald interval crosses a boundary:

Software usually hides this by clipping at 0 and 1. Clipping produces a plausible-looking interval and disguises the fact that the method has failed.

Part 2 — Why it fails, and what to use instead

The diagnosis

Wald inverts the wrong test. It builds the interval as "estimate ± error", using the standard error evaluated at the estimate $\hat{p}$ rather than at the hypothesised value $p$. When $\hat{p}$ is near a boundary, $\hat{p}(1-\hat{p})$ underestimates $p(1-p)$, the standard error comes out too small, and the interval is too narrow — in precisely the cases where the estimate is least reliable.

Wilson (1927)

Invert the test properly. Instead of solving for an interval around $\hat{p}$, find the set of $p$ that the data would not reject:

$$\left|\frac{\hat{p}-p}{\sqrt{p(1-p)/n}}\right|<z$$

Squaring and solving the resulting quadratic in $p$ gives a closed form:

$$\frac{\hat{p}+\frac{z^2}{2n}\pm z\sqrt{\frac{\hat{p}(1-\hat{p})}{n}+\frac{z^2}{4n^2}}}{1+\frac{z^2}{n}}$$

Uglier, entirely elementary, and it fixes the problem. It never leaves $[0,1]$, never produces a zero-width interval, and its coverage across the table above stays close to nominal everywhere.

Note what the algebra does: the centre is pulled from $\hat{p}$ towards $1/2$, by an amount that matters when $n$ is small and vanishes as $n$ grows. The interval is no longer symmetric about the point estimate — which is correct, because the sampling distribution of a proportion near a boundary is not symmetric either.

Agresti–Coull (1998)

An approximation to Wilson that can be explained in one sentence: add two successes and two failures, then run Wald.

$$\tilde{p}=\frac{x+2}{n+4},\tilde{p}\pm z\sqrt{\frac{\tilde{p}(1-\tilde{p})}{n+4}}$$

(Strictly it is $z^2/2\approx 1.92$ rather than 2, but 2 is the version worth remembering.) Coverage is good — sometimes conservative, as the table shows — and it has the practical advantage that you can implement it in a spreadsheet without explaining a quadratic to anyone.

Clopper–Pearson (1934)

The "exact" interval, obtained by inverting the binomial test directly. It guarantees at least nominal coverage for every $p$ — but that guarantee is bought with conservatism, often running at 98–99% when you asked for 95%, which means intervals wider than the data require.

Use it when under-coverage is unacceptable at any cost — regulatory submissions, safety reporting — and accept that you are giving away precision.

Which one

Wilson. It is the best-behaved across the realistic range, it is a closed form, and it is what the statistical literature has recommended for twenty-five years. Agresti–Coull if you need to explain the method to a non-technical audience. Clopper–Pearson if you are in a setting where guaranteed coverage is a requirement rather than a nicety.

Never Wald, at any base you would actually report on.

Part 3 — The survey complication

Everything above assumes simple random sampling. Real survey data is weighted, and the two problems compound.

The design-based standard error for a weighted proportion — derived in the companion article — is:

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

Substituting that into the Wald formula gives a design-correct standard error inside a method that is still wrong near the boundaries. You have fixed one problem and left the other.

The standard fix is a design-effect-adjusted Wilson interval: compute the design effect, convert to an effective sample size, and use $n_{\text{eff}}$ in the Wilson formula in place of $n$.

$$n_{\text{eff}}=\frac{n}{\text{deff}},\text{deff}=\frac{n\hat{V}(\hat{p})}{\hat{p}(1-\hat{p})}$$

This is what survey::svyciprop(..., method = "logit") and its relatives are doing, and it is why their intervals are asymmetric. Two caveats worth keeping in view:

• The effective sample size is itself estimated, and at small bases it is estimated badly. An $n_{\text{eff}}$ of 22 computed from a base of 40 carries real uncertainty that the interval does not show.
• The logit method — transform to the log-odds scale, build a symmetric interval there, transform back — is generally the better default in survey, and handles proportions near the boundary gracefully for the same structural reason Wilson does.

Part 4 — Implementation

R — base R already has Wilson, though the name gives no hint:

# prop.test uses the Wilson interval
prop.test(x = 3, n = 40, correct = FALSE)$conf.int
#> 0.02624 0.19965

# Wald, for comparison - note it dips below zero before clipping
p <- 3/40
p + c(-1, 1) * qnorm(0.975) * sqrt(p * (1 - p) / 40)
#> -0.00668  0.15668

# binom package for the full set
binom::binom.confint(3, 40,
  methods = c("asymptotic", "wilson", "agresti-coull", "exact"))

For weighted survey data:

library(survey)
des <- svydesign(ids = ~psu, strata = ~stratum, weights = ~w, data = dat)
svyciprop(~outcome, des, method = "logit")     # recommended default
svyciprop(~outcome, subset(des, region == "North"), method = "logit")

Python — statsmodels:

from statsmodels.stats.proportion import proportion_confint

proportion_confint(3, 40, alpha=0.05, method="wilson")
# (0.02624, 0.19965)
proportion_confint(3, 40, alpha=0.05, method="agresti_coull")
proportion_confint(3, 40, alpha=0.05, method="beta")     # Clopper-Pearson

Wilson from scratch, for a spreadsheet or a language without a library:

from math import sqrt

def wilson(x, n, z=1.959963985):
    """Wilson score interval for a binomial proportion."""
    p = x / n
    denom = 1 + z**2 / n
    centre = (p + z**2 / (2 * n)) / denom
    half = (z / denom) * sqrt(p * (1 - p) / n + z**2 / (4 * n**2))
    return centre - half, centre + half

SPSS — the base module does not offer Wilson for a single proportion; CTABLES confidence intervals are Wald. Compute Wilson in a COMPUTE block or move the calculation elsewhere.

Part 5 — Checklist

1. Default to Wilson. There is no case in survey reporting where Wald is preferable, and several where it is indefensible.
2. Never report a Wald interval on a base under 100, and treat 100–400 with low incidence as equally suspect.
3. Watch for zero-width intervals. If a table shows a CI of [0.0%, 0.0%] or [100%, 100%], the method has failed, not the population.
4. Watch for clipped intervals. A lower bound sitting exactly on 0.0% is usually a clipped negative number, not a real bound.
5. Combine the design correction with the boundary correction. Use design-effect-adjusted Wilson or the logit method; a design-correct Wald is still wrong.
6. Set a minimum reporting base and hold it. Below roughly 30–50 respondents, report the count rather than a percentage with an interval. "4 of 38 respondents" is honest; "10.5% ±9.8%" implies a precision the data does not have.
7. Quote the interval, not the flag. For low-incidence items the interval is the entire story, and it is usually much wider than readers expect.

Reproducing the coverage table

import numpy as np
from scipy import stats

z = stats.norm.ppf(0.975)

def wald(x, n):
    p = x / n
    h = z * np.sqrt(p * (1 - p) / n)
    return p - h, p + h

def wilson(x, n):
    p = x / n
    d = 1 + z**2 / n
    c = (p + z**2 / (2 * n)) / d
    h = (z / d) * np.sqrt(p * (1 - p) / n + z**2 / (4 * n**2))
    return c - h, c + h

def coverage(method, n, p):
    """Exact coverage: enumerate every outcome, weight by its probability."""
    xs = np.arange(n + 1)
    pmf = stats.binom.pmf(xs, n, p)
    return sum(pmf[x] for x in xs
               if method(x, n)[0] <= p <= method(x, n)[1])

print(coverage(wald, 20, 0.02))     # 0.3318
print(coverage(wilson, 20, 0.02))   # 0.9399

Sources and further reading

• Wilson, E.B. (1927) "Probable inference, the law of succession, and statistical inference", Journal of the American Statistical Association 22(158) — the original interval, and still the right default.
• Brown, L.D., Cai, T.T. & DasGupta, A. (2001) "Interval estimation for a binomial proportion", Statistical Science 16(2) — the paper that ought to have ended the use of Wald. Comprehensive, and the source of the "chaotic" description of Wald coverage.
• Agresti, A. & Coull, B.A. (1998) "Approximate is better than 'exact' for interval estimation of binomial proportions", The American Statistician 52(2) — the add-two-successes method, and the argument against reflexive use of exact intervals.
• Newcombe, R.G. (1998) "Two-sided confidence intervals for the single proportion: comparison of seven methods", Statistics in Medicine 17(8) — the systematic comparison, including intervals for differences of proportions.