Survey-Weighted Survival Analysis: Fusing svycoxph and rms

Why not just use cph()?

cph() uses the inverse observed information matrix for variance estimation, which assumes independent observations from a simple random sample. Applied to NHANES, this underestimates standard errors because it ignores correlation among participants sampled from the same PSU.

Why not just use svycoxph()?

svycoxph() produces correct point estimates and a sandwich variance-covariance matrix that accounts for the survey design. However, its output class does not support the rms generics: anova.rms(), Predict(), nomogram(), and survplot() do not dispatch on svycoxph objects. Testing nonlinearity of an RCS term requires manual use of survey::regTermTest(), and plotting smooth effects requires manual construction.

The fusion approach

The key insight is that anova.rms() and Predict() depend only on:

• The coefficient vector $coefficients
• The variance-covariance matrix $var
• The $Design object that rms attaches during cph() fitting, which maps coefficients to predictors and identifies spline structure

These components can be sourced from two separate fits and combined into a single object:

• Coefficients and vcov from svycoxph() — survey-correct
• Design and structural metadata from cph() — provides rms dispatch

Survival curve generation additionally requires a baseline hazard estimate. For this, survival::basehaz() applied to a cph() object uses an unweighted Breslow estimator. We implement a survey-weighted alternative following Lin (2000).

Key references

• Binder, D.A. (1992). Fitting Cox’s proportional hazards models from survey data. Biometrika, 79(1), 139–147.
• Lin, D.Y. (2000). On fitting Cox’s proportional hazards models to survey data. Biometrika, 87(1), 37–47.
• Lumley, T. (2004). Analysis of complex survey samples. Journal of Statistical Software, 9(1), 1–19.
• Harrell, F.E. (2015). Regression Modeling Strategies, 2nd ed. Springer. [rms package]

Comparing weighted and unweighted estimates

Both point estimates and standard errors differ between the two fits, and both differences matter. The coefficient differences (up to ~0.18 log-HR units for spline terms) mean that the unweighted cph() estimates represent the sample, not the US population. The SE ratios (1.15–1.61) confirm that ignoring the cluster design produces standard errors that are too small, with the largest design effect on race/ethnicity (RIDRETH1 ratio = 1.61), reflecting its strong geographic clustering.

tbl_coef <- data.frame(
  Term    = names(coef(fit_cph)),
  cph     = round(coef(fit_cph), 4),
  svy     = round(coef(fit_svy), 4),
  diff    = round(coef(fit_svy) - coef(fit_cph), 4),
  SE_cph  = round(sqrt(diag(vcov(fit_cph))), 4),
  SE_svy  = round(sqrt(diag(vcov(fit_svy))), 4),
  SE_ratio = round(sqrt(diag(vcov(fit_svy))) / sqrt(diag(vcov(fit_cph))), 3),
  row.names = NULL
)

flextable(tbl_coef) |>
  set_header_labels(
    Term     = "Term",
    cph      = "β (cph)",
    svy      = "β (svycoxph)",
    diff     = "Δβ",
    SE_cph   = "SE (cph)",
    SE_svy   = "SE (svycoxph)",
    SE_ratio = "SE ratio"
  ) |>
  colformat_double(digits = 4) |>
  bold(j = "SE_ratio", bold = TRUE) |>
  add_footer_lines("SE ratio > 1 indicates design effect from cluster sampling. Both β and SE differ materially, requiring substitution of both from svycoxph.") |>
  autofit()

Comparing anova.rms() output: survey-correct vs naive

The survey-correct Wald statistics are uniformly smaller because the sandwich vcov is larger. The magnitude of the difference reflects the design effect for each predictor. Race/ethnicity (RIDRETH1) shows the largest discrepancy (chi-square 8.60 vs 34.59, ratio ~4×), consistent with its strong geographic clustering within NHANES PSUs.

anova(fit_fused)   # survey-correct
anova(fit_cph)     # naive — overstates significance

# Predict() works: survey-correct CIs on the GGT smooth effect
p <- Predict(fit_fused, GGT = seq(5, 150, by = 5), fun = exp)
plot(p, ylab = "Hazard Ratio (vs median GGT)",
     xlab = "GGT (U/L)")

What works

The fusion approach is validated end to end on NHANES 1999–2018 data (N = 34,456 non-statin adults, 3,975 deaths). All three target outputs function correctly on the fused object:

Empirical demonstration of design effects

The NHANES cluster design inflates all Wald statistics in cph(). The magnitude of over-statement differs by predictor, reflecting how strongly each variable is geographically clustered within PSUs:

Race/ethnicity and BMI show the largest design effects (~4×), consistent with their strong geographic clustering. Sex shows the smallest (1.16×), as it is nearly uniformly distributed across PSUs. Ignoring survey design would lead to substantial over-confidence for these predictors.

Both GGT and albumin exhibit statistically robust nonlinearity under the survey-correct tests (GGT nonlinear χ² = 58.7, albumin nonlinear χ² = 25.6, both p < 0.0001), confirming that linear or categorical treatment of these biomarkers would misrepresent their mortality relationships.

When survey weighting is and is not necessary

The empirical design effects in the table above point toward a general principle with practical consequences for NHANES analysts.

The intraclass correlation coefficient (ICC) measures the proportion of a variable’s total variance that lies between PSUs rather than within them. It is the structural property that determines how much the cluster design inflates standard errors. The design effect is approximately \(\text{DEFF} \approx 1 + (m - 1) \cdot \text{ICC}\), where \(m\) is mean PSU size, and is directly estimable as the square of the observed SE ratio.

Race/ethnicity has a high ICC because people of the same race predominantly live in the same neighbourhoods and therefore the same PSUs. When race is the predictor of interest — as in studies of outcome disparities by race — the cluster design creates substantial extra-binomial variance in the score contributions, and unweighted standard errors are materially too small. Survey weighting is not optional in that setting.

Biochemical analytes measured in blood or urine occupy a different position. Their values are determined by individual physiology, not geography. The NHANES sampling probabilities are set by demographic strata, not by lab values, so conditional on the demographic adjusters included in the model, the correlation between survey weight and analyte value is near zero. This is the condition Harrell (2015) identifies for model-based (unweighted) estimates to be consistent: the sampling mechanism must be non-informative with respect to the variable of interest, given the covariates in the model. When that condition holds, the weighted and unweighted coefficient estimates converge and the design correction provides little benefit for coefficient inference.

The Lin (2000) design variance for the baseline hazard illustrates this directly: it is negligibly small precisely because the model has absorbed the major sources of geographic clustering through its covariates. A large Lin variance would indicate residual geographic heterogeneity not explained by the model — a signal of model misspecification, not just a reason to weight.

Efficient screening. Whether a given analyte warrants the full fusion machinery can be assessed cheaply before fitting the Cox model:

library(lme4)
library(performance)

# ICC for each analyte across PSUs — low ICC suggests non-informative sampling
analytes <- c("GGT", "LBXSAL", "TC", "BMI", "RIDAGEYR")
icc_tbl  <- lapply(analytes, function(v) {
  m   <- lmer(as.formula(paste(v, "~ 1 + (1|SDMVPSU)")), data = dat2, REML = TRUE)
  icc <- performance::icc(m)$ICC_adjusted
  data.frame(analyte = v, ICC = round(icc, 4))
})
do.call(rbind, icc_tbl)

# DEFF from design: (design SE / naive SE)^2 for each analyte mean
deff_tbl <- lapply(analytes, function(v) {
  se_design <- SE(svymean(reformulate(v), sub_design))
  se_naive  <- sd(dat2[[v]], na.rm = TRUE) / sqrt(sum(!is.na(dat2[[v]])))
  data.frame(analyte = v, DEFF = round((se_design / se_naive)^2, 3))
})
do.call(rbind, deff_tbl)

Analytes with ICC < 0.02 or DEFF < 1.1 are unlikely to require survey weighting for coefficient inference. Those with ICC > 0.05 or DEFF > 1.3 should use the full fusion approach. Race/ethnicity and other demographic variables used in the NHANES sampling design will always fall in the latter category and should never be analysed with unweighted standard errors when disparities are the focus.

This screening step — applied across all analytes of interest before modelling — is itself a useful contribution to NHANES analysis practice, as it allows analysts to apply the full svycph_fuse() workflow selectively where it matters rather than uniformly across all predictors.

Implications for prior literature. The ICC framework offers a more precise account of where prior NHANES survival analyses are methodologically vulnerable than a blanket criticism of survey weight omission. Biochemical analytes measured in blood or urine — including lipids (TC, LDL-C, HDL-C, triglycerides), liver enzymes (GGT, ALT, AST, alkaline phosphatase), renal markers (creatinine, albumin), and nutritional biomarkers — are determined by individual physiology rather than geography. Their ICC values are expected to be low, making non-survey-weighted coefficient estimates defensible for these predictors when demographic adjusters are included in the model. Prior analyses that omitted survey weighting while studying these analytes may therefore have produced valid conditional associations despite the omission.

The more consequential methodological limitations in that literature are elsewhere: the use of logistic regression on a binary mortality outcome when follow-up time and censoring information were available, discarding the entire time dimension of the data; and the discretization of continuous biomarkers into quantile or tertile categories. The latter is not merely an aesthetic limitation. For biomarkers with U-shaped mortality relationships — as observed here for GGT and albumin, and widely reported for TC — quantile categorization systematically biases toward the null: observations on the ascending and descending arms of the relationship fall within the same category and their opposing risk contributions cancel, attenuating the estimated effect and obscuring the shape entirely. This bias operates whether or not survey weights are applied, making the choice of continuous spline modelling the more consequential methodological decision. These criticisms are unconditional — they apply regardless of analyte type, sampling design, or weighting strategy, and they bear directly on the biological interpretability of reported findings.

Current limitations

Confidence band scale. weighted_basehaz() provides two variance options. The Lin (2000) design variance (se_type = "lin") is theoretically correct for population inference but produces nearly invisible bands for NHANES-scale data (~\(10^{-6}\) on the log scale) because PSU-selection uncertainty is negligible for rare events. The Greenwood-weighted option (se_type = "greenwood") gives interpretable band widths (~\(2 \times 10^{-4}\)) at the cost of reflecting population-scale rather than sample-scale precision. Neither exactly matches the sample-scale $std.err from cph() (~\(5 \times 10^{-3}\)), which treats the data as a simple random sample. The three-regime structure — design variance, population-scale statistical variance, and sample-scale statistical variance — reflects a genuine open question in the survival analysis literature: no unified variance estimator simultaneously captures PSU-selection uncertainty and finite-sample event-count uncertainty for weighted Cox models. Deriving such an estimator, likely via a joint influence function that propagates both sources of randomness, represents a natural extension of Lin (2000) and a direction for future methodological work.

Bootstrap validation and calibration. validate() and calibrate() from rms internally refit the model using cph(), not svycoxph(), so the survey correction is lost silently during resampling. These functions should not be used on fused objects without a custom resampling wrapper.

Singleton PSU strata. Subsetting NHANES data before creating the svydesign() object can produce strata with a single PSU, causing svycoxph() to fail. The correct approach is to create the design on the full dataset and subset the design object, as demonstrated in Step 1.

The nhanesR package

The functions described in this vignette — svycph_fuse(), weighted_basehaz(), svycph_set_basehaz(), and the ICC and DEFF screening utilities — are exported components of the nhanesR package. Together with the data infrastructure functions (nhanes_download_analyte(), nhanes_harmonize(), nhanes_mortality_link(), nhanes_survival_prep()), they constitute a complete pipeline from raw NHANES data download through survey-correct survival analysis within a single R package.

The design intention is that an analyst wishing to study the mortality relationship of any biochemical analyte measured across NHANES cycles can: (1) retrieve and harmonize the analyte across cycles using nhanes_download_analyte() with automatic CDC-catalog resolution of per-cycle filename changes; (2) link to NCHS public-use mortality data via nhanes_mortality_link(); (3) screen the analyte’s ICC to determine whether the full survey-correction machinery is warranted; and (4) fit a svycph_fuse() object that brings both the inferential correctness of svycoxph() and the display richness of the rms ecosystem to bear on the analysis, including spline nonlinearity tests, smooth effect plots, and survey-weighted survival curves.

The complete workflow demonstrated in this vignette — from NHANES download to survplot() with confidence bands — is fully reproducible using only functions in nhanesR and its declared dependencies.

A broader methodological agenda

The current implementation is presented as a practically useful approximation: it provides design-correct point estimates and coefficient inference via the fused object, and interpretable survival curves via the Greenwood-weighted baseline hazard, but it does not yet resolve the unified variance problem identified above. The JSS paper accompanying nhanesR documents this approximation honestly and provides the ICC/DEFF screening framework to guide analysts in deciding when the approximation is adequate.

A companion paper targeting Biostatistics would address the deeper unsolved problems. The central contribution would be a design-aware REML (dREML) criterion for penalized Cox regression: modifying the mgcv::gam() smoothing parameter selection with family = cox.ph() to use the survey-corrected sandwich variance rather than the observed information matrix. This is the analogue of what @lumley2004 did for svyglm() — adapting the estimating equations to the cluster design — but extended to the penalized regression setting where the smoothing parameter governs the shape of the estimated effect. A dREML-selected 2D tensor product smooth via te() would produce joint hazard ratio surfaces whose confidence regions properly account for the NHANES cluster design rather than treating PSU-clustered observations as independent; for race-stratified surface analyses, this distinction is material. The theory would be validated against the three-regime SE structure identified here, with a simulation study varying ICC from 0 to 0.15 to establish the empirical threshold below which the design correction leaves spline shape estimates unchanged.

The primary NHANES application — TC × HDL and GGT × albumin joint mortality surfaces — would showcase a visualization approach not yet standard in the biomarker-mortality literature. Each figure would combine three layers: (1) a smooth 2D hazard ratio surface from the design-correct GAM-Cox; (2) an empirical support envelope from Hmisc::perimeter(), which computes the convex-hull boundary of the joint observed predictor distribution and, when passed to bplot(perim = ...), restricts the displayed surface to the region where the model is not extrapolating; and (3) the observed mortality events plotted as points within that envelope. This three-layer construction — smooth surface, bounded by empirical support, annotated with actual deaths — immediately answers the question of whether the high-risk regions of the estimated surface correspond to where deaths actually occurred, and whether those deaths fall within clinically recognizable covariate ranges. The mgcv analogue is vis.gam(..., too.far = 0.1), which achieves the same masking automatically from the prediction grid. The event overlay is the novel element: to our knowledge no published 2D risk surface in the biomarker literature has overlaid the event locations within the perimeter-bounded surface, though the information is trivially available from the model’s response variable. It is a presentation innovation made straightforwardly available by the perimeter() function already present in Hmisc.