A public-health town: younger and older adults

Imagine a town of 1,000 adults. We deliberately sample older adults at a much higher rate so that we have enough of them to study.

In R: the toy weight

town <- data.frame(
  age_group = c("Younger", "Older"),
  in_town = c(900, 100),
  sampled = c(100, 100)
)

town$chance <- town$sampled / town$in_town
town$weight <- 1 / town$chance
town
#>   age_group in_town sampled    chance weight
#> 1   Younger     900     100 0.1111111      9
#> 2     Older     100     100 1.0000000      1

Hypertension prevalence: sample versus town

Suppose hypertension prevalence is 15% in younger adults and 40% in older adults.

• Sample: \((15+40)/2 = 27.5\%\).
• Town: \([900(0.15)+100(0.40)]/1000 = 17.5\%\).

In R: build the town sample

library(survey)

s <- data.frame(
  age_group = rep(c("Younger", "Older"), each = 100),
  weight = rep(c(9, 1), each = 100)
)

s$htn <- c(
  rep(1, 15), rep(0, 85),
  rep(1, 40), rep(0, 60)
)

mean(s$htn)
#> [1] 0.275

The ordinary mean gives the prevalence among the 200 sampled people.

In R: let the weights reconstruct the town

des_town <- svydesign(
  ids = ~1,
  weights = ~weight,
  data = s
)

svymean(~htn, des_town)
#>      mean     SE
#> htn 0.175 0.0327

Age can matter for two different reasons

Back to NHANES: the final weight is more than 1 / probability

The toy isolated the base selection weight. NHANES final weights also adjust for two additional features.

Which NHANES weight?

Different variables come from different components of NHANES.

MASLD: follow the measurement path

Combining ordinary two-year cycles

For 2007-2018, we combine six ordinary two-year cycles.

Each two-year weight represents a full population for its cycle, so the pooled weight is

\[ w_{2007-2018}=\frac{1}{6}w_{2yr}. \]

Dividing every weight by the same constant:

• changes weighted population totals;
• does not change weighted percentages, means, or most regression coefficients;
• does not change standard errors that are invariant to a common scaling of weights.

In R: pooling 1999-2004

nhanes$wt6 <- NA_real_

# 1999-2002: start with the NCHS 4-year weight
nhanes$wt6[nhanes$SDDSRVYR %in% c(1, 2)] <-
  (2 / 3) * nhanes$WTMEC4YR[nhanes$SDDSRVYR %in% c(1, 2)]

# 2003-2004: ordinary 2-year weight
nhanes$wt6[nhanes$SDDSRVYR == 3] <-
  (1 / 3) * nhanes$WTMEC2YR[nhanes$SDDSRVYR == 3]

Second special case: 2017-March 2020 is a 3.2-year file

The public-use 2017-March 2020 pre-pandemic file represents 3.2 years, not 2 years.

Combine 2015-2016 (2 years) with 2017-March 2020 (3.2 years): total period = 5.2 years.

\[ w^*_{2015-16}=\frac{2}{5.2}w_{2015-16} \]

\[ w^*_{2017-Mar2020}=\frac{3.2}{5.2}w_{2017-Mar2020}. \]

COVID-era example in code

nhanes$wt52 <- NA_real_

nhanes$wt52[nhanes$SDDSRVYR == 9] <-
  (2 / 5.2) * nhanes$WTMEC2YR[nhanes$SDDSRVYR == 9]

nhanes$wt52[nhanes$SDDSRVYR == 66] <-
  (3.2 / 5.2) * nhanes$WTMECPRP[nhanes$SDDSRVYR == 66]

For three components - 2013-14, 2015-16, and 2017-March 2020 - the total is 7.2 years:

• each ordinary 2-year cycle gets 2 / 7.2;
• the 3.2-year pre-pandemic file gets 3.2 / 7.2.

Three pooling situations to remember

Three design ingredients, three jobs

Why strata and PSUs matter

People in the same sampled geographic unit often resemble each other. Treating them as completely independent can make uncertainty too optimistic.

Stratification does something different: sampling PSUs separately within planned strata can increase precision.

Design-based and model-based inference: the useful distinction

Model-based analysis asks what follows under a statistical model for how outcomes are generated.

Design-based survey analysis also uses the known sampling design to make inference about a finite target population represented by the survey.

Survey weights are not treatment/IP weights

Build the survey design before defining the analysis subgroup

# 1. People eligible for the fasting-subsample weight
fast <- subset(nhanes, WTSAF2YR > 0)
fast$wt <- fast$WTSAF2YR / 6

# 2. Tell R how NHANES sampled them
library(survey)
des_full <- svydesign(
  ids = ~SDMVPSU,
  strata = ~SDMVSTRA,
  weights = ~wt,
  nest = TRUE,
  data = fast
)

A domain is just the subgroup you want to analyse

Examples:

• adults with MASLD;
• women;
• Group I;
• adults aged 60+;
• people with complete data for a particular model.

Wrong and right ways to analyse a subgroup

masld <- subset(fast, in_rows)

des_wrong <- svydesign(
  ids = ~SDMVPSU,
  strata = ~SDMVSTRA,
  weights = ~wt,
  nest = TRUE,
  data = masld
)

des_full <- svydesign(
  ids = ~SDMVPSU,
  strata = ~SDMVSTRA,
  weights = ~wt,
  nest = TRUE,
  data = fast
)

des <- subset(des_full, in_rows)

A small domain shows why the order matters

# Group I women: obese, lower waist fat
gw <- subset(fast, in_rows & group == "I" & female == 1)

c(people = nrow(gw), strata = length(unique(gw$SDMVSTRA)))
#> people strata
#>     57     44

tryCatch(svymean(~age, svydesign(
  ids = ~SDMVPSU, strata = ~SDMVSTRA,
  weights = ~wt, nest = TRUE, data = gw
)), error = function(e) cat("Error:", conditionMessage(e)))
#> Error: Stratum (59) has only one PSU at stage 1

svymean(~age, subset(des_full, in_rows & group == "I" & female == 1))
#>       mean     SE
#> age 41.711 2.3874

Complete cases are another domain - but that does not fix missing-data bias

fast$complete_model <- complete.cases(
  fast[, c("time_yr", "dead", "obese", "central",
           "age", "female", "race", "smk_current", "sed")]
)

# survey design already exists
des_full <- update(des_full, complete_model = fast$complete_model)
analysis_des <- subset(
  des_full,
  MASLD == 1 & complete_model
)

A weighted estimate can still be based on little information

A subgroup may represent millions of people but contain only a small number of respondents.

For descriptive survey estimates, inspect:

• unweighted n;
• the design effect and effective sample size;
• design degrees of freedom;
• confidence-interval width.

Table 1 for a complex survey

A useful reporting convention is:

• report unweighted n;
• report weighted percentages or weighted means;
• report design-based uncertainty when appropriate.

Example:

Fit the model on the design object

fm <- Surv(time_yr, dead) ~ obese * central +
  age + female + race + smk_current + sed

fit_w <- svycoxph(fm, design = des)

The design object supplies:

• the analysis rows in the domain;
• their survey weights;
• strata and PSU information for variance estimation.

Tests use the survey design too

regTermTest(
  fit_w,
  ~obese:central,
  df = degf(des)
)
#> Wald test for obese:central
#>  in svycoxph(formula = fm, design = des)
#> F =  0.3185855  on  1  and  94  df: p= 0.5738

For our pooled design:

\[ \text{design df} = \text{PSUs} - \text{strata} = 184 - 90 = 94. \]

Weighting changes the real MASLD population mix only modestly

What can we say about the published Cox analysis?

The paper describes NHANES as a stratified, multistage survey designed for national inference.

Its reported statistical methods do not appear to specify survey weights, strata, or PSUs, and an ordinary Cox model closely reproduces the published unadjusted Group IV vs I hazard ratio.

Start with one comparison

High vs lower waist fat among the non-obese (Group IV vs III), adjusted for Week 4’s L:

Adjustment and survey design are different choices

• Left -> right within weighting status: confounder adjustment.
• Unweighted -> weighted within adjustment status: survey design / represented population.

Week 5’s effect-modification result, redone with the survey design

The survey-design analysis gives less evidence and much more uncertainty. That is not evidence of no heterogeneity - the interval is very wide.

In R: the Week 5 interaction model, survey-aware

fit_u <- coxph(fm, data = d)          # ordinary model
fit_w <- svycoxph(fm, design = des)   # NHANES design

waist_hr <- function(fit, obese) {
  b <- coef(fit)
  exp(b[["central"]] + obese * b[["obese:central"]])
}

round(c(
  unweighted_not_obese = waist_hr(fit_u, 0),
  weighted_not_obese   = waist_hr(fit_w, 0),
  unweighted_obese     = waist_hr(fit_u, 1),
  weighted_obese       = waist_hr(fit_w, 1)
), 2)
#> unweighted_not_obese   weighted_not_obese
#>                 1.09                 1.43
#>     unweighted_obese       weighted_obese
#>                 2.99                 2.21

The correct weight is chosen by the measurement process

Same 5,911 rows; waist fat adjusted for L.

In R: sensitivity to the exam weight

mec <- subset(nhanes, WTMEC2YR > 0)
mec$wt <- mec$WTMEC2YR / 6

des_mec <- subset(
  svydesign(
    ids = ~SDMVPSU, strata = ~SDMVSTRA,
    weights = ~wt, nest = TRUE, data = mec
  ),
  in_rows
)

fit_m <- svycoxph(fm, design = des_mec)

round(rbind(
  fasting = c(not_obese = waist_hr(fit_w, 0), obese = waist_hr(fit_w, 1)),
  exam    = c(not_obese = waist_hr(fit_m, 0), obese = waist_hr(fit_m, 1))
), 2)
#>         not_obese obese
#> fasting      1.43  2.21
#> exam         1.46  2.10

Why can weighted and unweighted regression differ?

Three broad reasons are worth considering:

1. Sampling variation: two noisy estimators can differ by chance.
2. Effect heterogeneity: the effect differs across age/race/income/etc., and weighting changes their population mix.
3. Sampling-model dependence: variables related to sampling or participation are not adequately represented in the outcome model.

Why weighting can change an overall RR

Hypothetical 5-year hypertension risks:

Week 5 + Week 6: effect modification by age + unequal sampling of age groups. High-sodium diet is 50% in both age groups here, so age does not confound this toy comparison.

In R: build the high-sodium example

# 100 younger + 100 older; half exposed in each age group
s2 <- data.frame(
  age_group = rep(c("Younger", "Older"), each = 100),
  high_sodium = rep(rep(c(0, 1), each = 50), times = 2),
  weight = rep(c(9, 1), each = 100)
)

cases <- function(k, n = 50) c(rep(1, k), rep(0, n - k))
s2$incident_htn <- c(
  cases(5), cases(10),    # younger: 10% -> 20%
  cases(10), cases(30)    # older:   20% -> 60%
)

des2 <- svydesign(ids = ~1, weights = ~weight, data = s2)

In R: sample RR versus town RR

fit_sample <- glm(
  incident_htn ~ high_sodium,
  family = poisson(link = "log"), data = s2
)

fit_town <- svyglm(
  incident_htn ~ high_sodium,
  family = quasipoisson(link = "log"), design = des2
)

round(exp(c(
  sample = coef(fit_sample)[["high_sodium"]],
  town = coef(fit_town)[["high_sodium"]]
)), 2)
#> sample   town
#>   2.67   2.18

What survey design cannot fix

Survey weighting/design does not repair:

• confounding of the exposure-outcome comparison;
• selection created by defining the MASLD cohort using variables related to the exposures;
• bias from complete-case analysis due to item missingness;
• poorly defined exposures or causal estimands;
• sparse cells - four deaths remain four deaths after weighting.

Report the same decisions you made

The algorithm to take into your own paper

Before trusting a result, ask:

• Is the weight valid for every variable?
• Did cycle pooling follow the correct time period?
• Was the full survey design declared before subgroup restriction?
• Is there enough empirical information in the subgroup?
• Does the population claim match the analysis?

Key takeaways

1. NHANES is not a simple random sample; weights, strata and PSUs each have a job.
2. Choose the weight from the smallest component needed by your variables, not from the answer it produces.
3. Pool weights according to what each supplied weight represents: ordinary 2-year cycles use / number of cycles; analyses including 1999-2000 must start with the NCHS 1999-2002 4-year weight; the 2017-March 2020 3.2-year file uses time-proportional scaling.
4. Declare the design first, then subset the domain. Complete cases are another domain, but that does not solve missing-data bias.
5. Weighted regression and unweighted regression may differ for several reasons; the gap is a question to investigate, not a diagnosis.
6. Survey design changes who is represented and how uncertainty is calculated. It does not make a causal comparison valid.

Sources for the core workflow

Selected references - survey foundations

Selected references - regression and reporting

Appendix: design effect in the MASLD model

Waist fat among the not obese, adjusted for L.

Appendix: variance estimation methods

Two broad families are common in complex surveys:

Taylor linearization

• approximates the variance of a complex estimator with a linearized form;
• uses design strata and PSUs;
• is the standard approach used in the main NHANES examples in this lecture.

Replication methods

• repeatedly recompute the statistic under replicate weights;
• examples: jackknife, BRR, Fay’s BRR, survey bootstrap.

Appendix: design-adjusted tests

Complex survey analysis changes more than standard errors.

Examples include:

• Rao-Scott chi-square instead of an ordinary Pearson chi-square for contingency tables;
• design-based Wald/F tests for regression terms;
• survey-adapted goodness-of-fit procedures such as Archer-Lemeshow for logistic models;
• survey-aware versions of model summaries such as pseudo-R2 when appropriate.

Appendix: survey-weighted survival curves

Group IV has the lowest survival in both ordinary and survey-weighted curves. The weighted curves shift because the represented population is younger and differently distributed across sampled groups.

Appendix: R code for survey survival curves

km <- svykm(Surv(time_yr, dead) ~ group, design = des)

died6 <- function(k) {
  100 * (1 - k$surv[max(which(k$time <= 6))])
}

round(sapply(km, died6), 1)
#>    I   II  III   IV
#>  1.2  5.1  4.5 11.3

lr <- svylogrank(Surv(time_yr, dead) ~ group, design = des)
round(lr[[2]], 4)
#>   Chisq       p
#> 41.5624  0.0000

Appendix: Week 5’s sex-modification analysis, design-aware

fs <- Surv(time_yr, dead) ~
  IV * female + age + race + smk_current + sed

d34 <- subset(des, group %in% c("III", "IV"))
fit_s <- svycoxph(fs, design = d34)

round(exp(c(
  ratio = coef(fit_s)[["IV:female"]],
  confint(fit_s)["IV:female", ]
)), 2)
#>  ratio  2.5 % 97.5 %
#>   2.60   0.89   7.58

regTermTest(fit_s, ~IV:female, df = degf(des))
#> Wald test for IV:female
#>  in svycoxph(formula = fs, design = d34)
#> F =  3.067102  on  1  and  94  df: p= 0.083154

Appendix: the debate about regression weights

There is broad agreement that survey weights are needed for population descriptive estimates.

Regression is more nuanced because the coefficient can depend on:

• the target estimand;
• whether sampling is informative after conditioning on modeled covariates;
• model misspecification;
• effect heterogeneity;
• efficiency costs from highly variable weights.

Appendix: extreme weights and trimming

Highly variable survey weights can make a small number of respondents extremely influential and increase variance.

Possible responses include:

• inspect whether the large weights are expected from the documented survey design;
• verify that the correct component weight was chosen;
• report the weight distribution and design effect;
• consider trimming only as a sensitivity analysis when scientifically justified.

Appendix: sample and population causal effects - SATE and PATE

When the estimand is genuinely causal, survey generalization is sometimes described using:

• SATE: sample average treatment effect - average causal effect in the sampled people;
• PATE: population average treatment effect - average causal effect in the target population.

Appendix: sample target vs population target

The old shorthand “weighted = population, unweighted = sample” is most literal for descriptive statistics.

For regression, distinguish the target explicitly:

• Sample-average / model-based association: what does the fitted model say for the observed analytic sample under its assumptions?
• Population-representative survey analysis: what coefficient/contrast is obtained when the survey design is used to represent the target population?

Appendix: extra COVID-era pooling example

Combine 2013-14, 2015-16, and 2017-March 2020.

Total represented period:

\[ 2 + 2 + 3.2 = 7.2\text{ years}. \]

So:

\[ w^*_{2013-14}=\frac{2}{7.2}w_{2yr}, \qquad w^*_{2015-16}=\frac{2}{7.2}w_{2yr}, \]

\[ w^*_{2017-Mar2020}=\frac{3.2}{7.2}w_{prepandemic}. \]

Appendix: COVID-era pooling in code

nhanes$wt72 <- NA_real_

nhanes$wt72[nhanes$SDDSRVYR %in% c(8, 9)] <-
  (2 / 7.2) * nhanes$WTMEC2YR[nhanes$SDDSRVYR %in% c(8, 9)]

nhanes$wt72[nhanes$SDDSRVYR == 66] <-
  (3.2 / 7.2) * nhanes$WTMECPRP[nhanes$SDDSRVYR == 66]

Check the relevant data-file documentation before assuming the same variable definition or weight name exists in every cycle.

Appendix: reliability standards - the formal idea

For descriptive proportions, NCHS reliability guidance considers multiple features, not only raw n.

Common diagnostics include:

• unweighted sample size;
• effective sample size;
• design degrees of freedom;
• confidence-interval width / relative precision;
• whether the estimate is close to the boundary 0 or 1.

Appendix: what if a lonely PSU truly remains?

Sometimes a stratum has only one PSU even after the design is declared correctly. Then the variance for that stratum cannot be estimated in the usual way.

Possible strategies depend on the survey documentation and software:

• treat a certainty PSU appropriately when that is how the design was constructed;
• combine or adjust singleton strata using a documented rule;
• use an approved replication method when replicate weights/design information are available.

Appendix: other complex surveys

The same principles recur in many public-health datasets:

• Canadian Community Health Survey (CCHS)
• Behavioral Risk Factor Surveillance System (BRFSS)
• Korean NHANES (KNHANES)
• Health and Retirement Study (HRS)
• National Comorbidity Survey Replication (NCS-R)
• European Social Survey (ESS)

Appendix: software

Common software has dedicated complex-survey procedures:

• R: survey package (svydesign, svymean, svyglm, svycoxph, svykm)
• Stata: svyset, svy: estimators
• SAS: SURVEY procedures
• SUDAAN: survey-specific procedures
• SPSS: Complex Samples procedures

Appendix: future extensions

Survey principles also matter for modern prediction workflows.

Examples include:

• sampling weights in machine-learning models;
• survey-aware cross-validation;
• calibration and performance measures for complex samples.

These require more than simply passing a weight column to an algorithm and are beyond the core Week 6 goals.