cassionData Analysis

Back to the lessonLesson 5 of 8Estimating properly

Declare the design once

The same deck as the downloads, rendered as a page. Start the slideshow to present it full screen — arrow keys or a click advance one slide, Escape leaves.

Slides · PDFSlides · PowerPoint

  1. Slide 1 / 21

    What this lesson covers

    • Stop carrying the design by hand
    • The design object in R
    • The same thing in Python
    • The ultimate cluster approximation
    • The four ways to get a plain average by accident
    • Estimates other than proportions
    • The sub-sample design
    • Write the design down beside the estimates
    • What comes next
    Speaker notes
    One design object, and every estimate after it inherits the weights, the strata and the clusters. The R survey package, the equivalent arithmetic in Python, and the four mistakes that silently produce a plain average.

  2. Slide 2 / 21

    Stop carrying the design by hand

    • Everything so far has been computed explicitly, which was the point — you should know what the arithmetic is.
    Speaker notes
    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.

  3. Slide 3 / 21

    The design object in R — 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

  4. Slide 4 / 21

    The design object in R

    • ids is the primary sampling unit — the enumeration area, not the household. Passing ~household_id here is the…
    • 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…
    Speaker notes
    Four arguments and every one is load-bearing. Then every estimate comes off the design.

  5. Slide 5 / 21

    The design object in R — In R

    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")
    Speaker notes
    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.

  6. Slide 6 / 21

    The same thing in Python — In Python (cont.)

    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)
    Speaker notes
    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.

  7. Slide 7 / 21

    The same thing in Python — In Python (cont.)

            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")

  8. Slide 8 / 21

    The same thing in Python — In Python (cont.)

    print(design.mean(lambda d: d["food_insecure"] == "true"))

  9. Slide 9 / 21

    The same thing in Python — In R

    # The R equivalent of the whole class above:
svymean(~I(food_insecure == "true"), design, deff = TRUE)
    Speaker notes
    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.

  10. Slide 10 / 21

    The ultimate cluster approximation

    • Both implementations use the same shortcut, and it is worth naming because it looks like a simplification and is not.
    Speaker notes
    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.

  11. Slide 11 / 21

    The four ways to get a plain average by accident

    • Passing the household as the cluster — ids = ~household_id says each household is its own PSU, so there is no…
    • Forgetting the weights — svydesign(..., weights = NULL) gives every household equal weight, which is the 32.8% from…
    • Filtering the data frame instead of subsetting the design — This one is subtle and common:
    Speaker notes
    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:

  12. Slide 12 / 21

    The four ways to get a plain average by accident — In R

    # 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")

  13. Slide 13 / 21

    The four ways to get a plain average by accident — In Python

    # 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")

  14. Slide 14 / 21

    The four ways to get a plain average by accident

    • Computing on a merged frame — Join the survey to anything one-to-many and the weights are now wrong, because a…
    Speaker notes
    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.

  15. Slide 15 / 21

    Estimates other than proportions — In R

    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)

  16. Slide 16 / 21

    Estimates other than proportions — In Python

    print(design.mean(lambda d: d["household_size"] > 7))
    Speaker notes
    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.

  17. Slide 17 / 21

    The sub-sample design — In R

    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")
    Speaker notes
    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.

  18. Slide 18 / 21

    The sub-sample design — In Python

    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))
    Speaker notes
    The clusters are the same areas, so the clustering still applies. Only the weight and the population change.

  19. Slide 19 / 21

    Write the design down beside the estimates — Example

    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
    Speaker notes
    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.

  20. Slide 20 / 21

    What comes next

    • The design object assumes everyone you selected was interviewed and every area you drew was visited.
    Speaker notes
    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.

  21. Slide 21 / 21

    Where this goes next

    Read the full lesson, with runnable code Back to the lesson