Declare the design once — Survey Analysis, Sampling and Weighting — Data Analysis

Stop carrying the design by hand

Everything so far has been computed explicitly, which was the point — you should know what the arithmetic is. From here it should be declared once and inherited, because the alternative is remembering to apply weights, strata and clusters to every estimate in a fifty-line script, and the estimate you forget is the one that goes in the summary.

The design object in R

library(survey)

design <- svydesign(
  ids     = ~ea_id,      # the primary sampling unit
  strata  = ~stratum,    # the stratification
  weights = ~weight,     # the weight built in lesson 2
  data    = survey,
  nest    = TRUE         # cluster ids are nested within strata
)

design

Four arguments and every one is load-bearing.

ids is the primary sampling unit — the enumeration area, not the household. Passing ~household_id here is the commonest error in this course and it silently produces a simple-random-sample variance.
strata removes the between-stratum variance, which is the gain the previous lesson measured.
weights makes the point estimate unbiased.
nest = TRUE tells the package that area identifiers are unique only within a stratum. Omit it where they are reused across strata and the areas get merged.

Then every estimate comes off the design.

svymean(~I(food_insecure == "true"), design, deff = TRUE)
svytotal(~I(food_insecure == "true"), design)
svyby(~I(food_insecure == "true"), ~stratum, design, svymean, vartype = c("se", "ci"))
svyciprop(~I(food_insecure == "true"), design, method = "logit")

svyby is the workhorse: any estimate, by any grouping, with the design carried through. And svyciprop is what you want for a proportion — the next lesson explains why the interval from svymean is the wrong shape near 0 and 1.

The same thing in Python

Python has no equivalent of survey that this sector has settled on, so the examples compute the estimator directly. That is a cost and also a benefit: the arithmetic is visible.

import numpy as np
import pandas as pd


class SurveyDesign:
    """Stratified two-stage design, ultimate-cluster variance."""

    def __init__(self, data, ids, strata, weights):
        self.data, self.ids, self.strata, self.weights = data, ids, strata, weights

    def mean(self, indicator):
        d = self.data
        y = indicator(d).astype(float)
        w = d[self.weights]
        p = np.average(y, weights=w)

        variance, total_w = 0.0, w.sum()
        for _, stratum in d.groupby(self.strata):
            residual = stratum[self.weights] * (indicator(stratum).astype(float) - p)
            totals = residual.groupby(stratum[self.ids]).sum()
            n = len(totals)
            if n < 2:
                continue
            variance += n / (n - 1) * ((totals - totals.mean()) ** 2).sum()
        se = np.sqrt(variance) / total_w

        srs = np.sqrt(p * (1 - p) / len(d))
        return {"estimate": p, "se": se, "deff": (se / srs) ** 2,
                "ci_low": p - 1.96 * se, "ci_high": p + 1.96 * se}


design = SurveyDesign(survey, ids="ea_id", strata="stratum", weights="weight")
print(design.mean(lambda d: d["food_insecure"] == "true"))

# The R equivalent of the whole class above:
svymean(~I(food_insecure == "true"), design, deff = TRUE)

The Python class is thirty lines and the R call is one. That asymmetry is real and it is why survey analysis in this sector is usually done in R. Use R for survey estimation where you have the choice; the Python path is for teams that do not, and for understanding what the R call is doing.

The ultimate cluster approximation

Both implementations use the same shortcut, and it is worth naming because it looks like a simplification and is not.

The variance is computed from the totals of the primary sampling units, ignoring the second stage entirely. That is exact when the first stage is sampled with replacement, and slightly conservative otherwise — it very slightly overstates the variance because it ignores the finite population correction at stage two.

Slightly conservative is the right direction to be wrong in, which is why every major package does this by default. survey calls it the ultimate cluster approximation and so does every DHS report.

The four ways to get a plain average by accident

Each of these produces a number that looks fine and is a simple random sample estimate.

Passing the household as the cluster. ids = ~household_id says each household is its own PSU, so there is no clustering to account for and the standard error collapses.

Forgetting the weights. svydesign(..., weights = NULL) gives every household equal weight, which is the 32.8% from lesson 1.

Filtering the data frame instead of subsetting the design. This one is subtle and common:

# WRONG: the design object no longer knows about the dropped strata
rural <- svydesign(ids = ~ea_id, strata = ~stratum, weights = ~weight,
                   data = filter(survey, stratum != "urban"), nest = TRUE)

# RIGHT
rural <- subset(design, stratum != "urban")

# Same rule: subset the design, keeping the full strata structure in view
rural = SurveyDesign(survey[survey["stratum"] != "urban"],
                     ids="ea_id", strata="stratum", weights="weight")

Subsetting a design keeps the strata and cluster structure intact for variance estimation. Rebuilding it from a filtered frame throws away areas with no remaining households, which changes the degrees of freedom and can silently produce a stratum with one cluster in it — where the variance contribution is undefined and quietly dropped.

Computing on a merged frame. Join the survey to anything one-to-many and the weights are now wrong, because a household appears several times. Aggregate to the household before estimating, which is the grain lesson from the joining course arriving in a new setting.

Estimates other than proportions

svymean(~household_size, design, na.rm = TRUE)
svytotal(~I(food_insecure == "true"), design)
svyquantile(~minutes_to_water, design, quantiles = c(0.5, 0.9), na.rm = TRUE)
svyratio(~I(food_insecure == "true"), ~I(children_under5 > 0), design)

print(design.mean(lambda d: d["household_size"] > 7))

Two notes. A weighted total is the estimate people most often want and most often cannot get, because it needs population-scaled weights — the normalised weights from lesson 2 will give a total of 996. And na.rm = TRUE on a survey mean is a decision, not a convenience: it assumes the missing households resemble the observed ones within their stratum, which is the same assumption the non-response adjustment makes and should be stated once for both.

The sub-sample design

The measured children are a sub-sample with their own weight, and they need their own design object — not a filter on the household one.

children <- subset(design, !is.na(child_muac_mm))
children <- update(children, child_weight = weight * children_under5)

child_design <- svydesign(ids = ~ea_id, strata = ~stratum,
                          weights = ~child_weight,
                          data = subset(survey, !is.na(child_muac_mm)), nest = TRUE)

svyciprop(~I(child_muac_mm < 125), child_design, method = "logit")

measured = survey[survey["child_muac_mm"].notna()].copy()
measured["child_weight"] = measured["weight"] * measured["children_under5"]

child_design = SurveyDesign(measured, ids="ea_id", strata="stratum",
                            weights="child_weight")
print(child_design.mean(lambda d: d["child_muac_mm"] < 125))

The clusters are the same areas, so the clustering still applies. Only the weight and the population change.

Write the design down beside the estimates

Design           Stratified two-stage cluster
Strata           urban, rural-accessible, rural-remote
PSU              enumeration area (75 sampled, 25 per stratum)
Stage 1          PPS, size measure = households at frame listing
Stage 2          SRS, 14 households per area
Weights          base = M_h / (25 * 14), non-response adjusted within area
Variance         ultimate cluster, Taylor linearisation
Software         R survey 4.x / equivalent direct computation in Python

Seven lines. Any reader can now reproduce your standard errors, and any reviewer can tell in ten seconds whether you did it properly — which is more than can be said for most survey annexes.

What comes next

The design object assumes everyone you selected was interviewed and every area you drew was visited. Neither was true here. The next lesson handles the 54 households that were not interviewed and the two areas that had to be replaced.