--- title: "Survey-Weighted Survival Analysis: Fusing svycoxph and rms" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Survey-Weighted Survival Analysis: Fusing svycoxph and rms} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE ) ``` ## Background and Motivation NHANES is a complex probability sample of the US civilian non-institutionalized population. Valid inference from NHANES survival data requires three things that standard Cox model implementations do not provide jointly: 1. **Survey-weighted partial likelihood** — accounts for unequal selection probabilities, stratification, and clustering within primary sampling units (PSUs). Ignoring these produces standard errors that are too small and confidence intervals that are too narrow. 2. **Flexible nonlinear covariate effects** — biomarkers such as total cholesterol, GGT, and albumin have U-shaped or threshold relationships with mortality. Restricted cubic splines (RCS) capture these shapes without requiring a priori categorization. 3. **A rich output environment** — the `rms`/`Hmisc` ecosystem provides nonlinearity tests (`anova.rms()`), effect displays (`Predict()`, `plot.Predict()`), and survival curve plots (`survplot()`) that go far beyond what base `survival` offers. The `survey` package provides `svycoxph()` for (1); the `rms` package provides `cph()` for (2) and (3). Neither alone provides all three. The goal of this vignette is to describe and implement a fusion approach that combines the survey-correct inference of `svycoxph()` with the output machinery of the `rms` ecosystem. ### 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] --- ## Step 1: Examine object structures Before implementing the fusion, we need to understand exactly which slots in a `cph()` object are used by each `rms` generic, and which corresponding components are available in a `svycoxph()` object. We fit both models on the nhanesR analytic dataset using the same formula and compare structures. ```{r load-data} library(nhanesR) library(rms) library(survey) library(survival) library(flextable) dat <- readRDS("~/Documents/R.code/nhanesR/analytic_survival.rds") # Analysis population: non-statin users, adults >= 20, landmark > 2yr, # complete GGT / albumin / TC / BMI / PIR dat2 <- subset(dat, ELIGSTAT == 1 & !is.na(time) & time > 2 & statin == FALSE & !is.na(GGT) & !is.na(LBXSAL) & !is.na(TC) & !is.na(BMI) & !is.na(INDFMPIR) & RIDAGEYR >= 20 ) # N = 34,456 events = 3,975 ``` ```{r survey-design} # NHANES design: create on the FULL dataset, then subset the design object. # Creating the design on the already-subsetted data can leave some strata with # a single PSU, causing svycoxph() to fail at variance estimation. full_design <- svydesign( ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, nest = TRUE, data = dat ) sub_design <- subset(full_design, ELIGSTAT == 1 & !is.na(time) & time > 2 & statin == FALSE & !is.na(GGT) & !is.na(LBXSAL) & !is.na(TC) & !is.na(BMI) & !is.na(INDFMPIR) & RIDAGEYR >= 20 ) ``` ```{r fit-models} # Shared formula: RCS(4 knots) on GGT and albumin; linear adjusters f <- Surv(time, event) ~ rcs(GGT, 4) + rcs(LBXSAL, 4) + RIDAGEYR + RIAGENDR + RIDRETH1 + log(BMI) # cph() fit — x=TRUE, y=TRUE, surv=TRUE needed for Predict() and survplot() # Inference from this fit is NOT survey-correct; used only for $Design structure dd <- datadist(dat2) options(datadist = "dd") fit_cph <- cph(f, data = dat2, x = TRUE, y = TRUE, surv = TRUE) # svycoxph() fit — survey-correct coefficients and sandwich vcov fit_svy <- svycoxph(f, design = sub_design) ``` ```{r examine-cph-structure} names(fit_cph) ``` ```{r examine-svy-structure} names(fit_svy) ``` ### 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. ```{r compare-coefs-table} 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() ``` --- ## Step 2: Identify Design slot dependencies The `rms` generics dispatch based on the `Design` attribute attached to `cph` objects during fitting. Inspection reveals that `cph` and `svycoxph` share a common core of slots but diverge in the metadata needed by each ecosystem. Slots present in `cph` but **not** `svycoxph` — the structural components we must preserve from `cph`: | Slot | Purpose | |---|---| | `$Design` | Maps coefficients to terms; identifies spline nonlinear components | | `$surv` | Baseline survival S_0(t) as numeric vector (not a list) | | `$time` | Time points corresponding to `$surv` | | `$std.err` | SEs of baseline survival (used by `survplot()` for confidence bands) | | `$maxtime` | Maximum observed time | | `$time.inc` | Time axis increment for `survplot()` | | `$center`, `$scale.pred` | Centering constants for `Predict()` | | `$x`, `$y` | Design matrix and survival outcome (needed for `Predict()`) | Slots present in `svycoxph` but **not** `cph` — the survey-correct inference components we substitute in: | Slot | Purpose | |---|---| | `$var` | Sandwich vcov (replaces naive information-matrix vcov) | | `$naive.var` | Naive vcov, stored separately | | `$degf.resid` | Survey df = n_PSU − n_strata (138 for this analysis population) | | `$survey.design` | The `svydesign` object used for fitting | | `$weights` | Participant sampling weights | **Coefficient naming** differs between the two fits. `rms` shortens `rcs(GGT, 4)GGT` to `GGT` and `log(BMI)` to `BMI`; `svycoxph` preserves the full formula terms. The positional order is identical (same formula), so `svycph_fuse()` copies values by position and applies `cph` names. ```{r design-structure} str(fit_cph$Design, max.level = 2) ``` ```{r slot-names} # Confirm positional correspondence; names will differ length(coef(fit_cph)) == length(coef(fit_svy)) # TRUE names(coef(fit_cph)) # rms short names names(coef(fit_svy)) # full formula names ``` --- ## Step 3: Implement the fusion The `svycph_fuse()` function takes a fitted `cph` object and a fitted `svycoxph` object with the same formula, and returns a modified `cph`-class object with survey-correct coefficients and vcov. ```{r svycph-fuse-source} # Source the implementation (see R/svycph_fuse.R) # devtools::load_all("~/Documents/R.code/nhanesR") ``` ```{r apply-fusion} fit_fused <- svycph_fuse(fit_cph, fit_svy) ``` A note on coefficient naming: `rms` shortens `rcs(GGT, 4)GGT` to `GGT` and `log(BMI)` to `BMI`, while `svycoxph` preserves the full term names. `svycph_fuse()` copies values by position (same formula guarantees same order) and applies `cph` names so that downstream `rms` generics resolve terms correctly. ### 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. ```{r test-anova} anova(fit_fused) # survey-correct anova(fit_cph) # naive — overstates significance ``` ```{r test-predict} # 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)") ``` --- ## Step 4: Survey-weighted baseline hazard `survplot()` requires a baseline hazard estimate. The standard `survival::basehaz()` uses the unweighted Breslow estimator; we implement the survey-weighted version following Lin (2000). The weighted cumulative baseline hazard at event time $t_i$ is: $$\hat{H}_0^w(t) = \sum_{t_i \leq t} \frac{w_i}{\sum_{j \in \mathcal{R}(t_i)} w_j \exp(\mathbf{X}_j^\top \hat{\boldsymbol{\beta}})}$$ where $w_i$ is the survey weight and $\hat{\boldsymbol{\beta}}$ comes from `svycoxph()`. The `weighted_basehaz()` function computes both the weighted Breslow point estimate and its standard error. The `se_type` argument controls the variance estimator: **`se_type = "lin"` (default)** implements Lin (2000) eq. 2.4 — the linearization variance via PSU-level totals of the influence function $\Phi_i(t) = \sum_{t_k \leq t} \phi_i(t_k)$. This is the design-based variance: it measures how much $\hat{H}_0^w(t)$ would change if different PSUs had been selected. For NHANES-scale populations with rare events, this is very small (~$10^{-6}$ on the log scale) because event rates are similar across PSUs. **`se_type = "greenwood"`** uses the survey-weighted Greenwood formula: $\sum_{t_k \leq t} n^w(t_k) / [Y^w(t_k)]^2$. This measures statistical precision from the weighted event count — the population-scale analog of the Nelson-Aalen variance. For `survplot()` confidence bands this gives widths proportional to population-scale uncertainty (~$2 \times 10^{-4}$ on the log scale). The unweighted `cph()` `$std.err` (~$5 \times 10^{-3}$) reflects sample-scale precision from $n \approx 34{,}000$ observations and $\sim 4{,}000$ events, treated as if they were a simple random sample. For visualization, `se_type = "greenwood"` produces the most interpretable confidence bands; `se_type = "lin"` is appropriate for formal design-based inference. ```{r weighted-basehaz} # Lin design variance (default) — correct for population inference h0_lin <- weighted_basehaz(fit_svy, design = sub_design, se_type = "lin") # Greenwood-weighted — interpretable survplot() confidence bands h0_gw <- weighted_basehaz(fit_svy, design = sub_design, se_type = "greenwood") head(h0_gw) ``` ```{r compare-se-scale} # SE scale comparison (log H0 scale, late follow-up): # cph unweighted std.err ~ 0.005 (sample-scale statistical precision) # Greenwood-weighted ~ 0.0002 (population-scale statistical precision) # Lin design ~ 1e-6 (PSU-selection uncertainty) data.frame( method = c("cph unweighted", "Greenwood-weighted", "Lin design"), std.err = c( mean(tail(fit_cph$std.err, 5), na.rm = TRUE), mean(tail(h0_gw$std.err, 5)), mean(tail(h0_lin$std.err, 5)) ) ) ``` ```{r compare-basehaz} h0_naive <- basehaz(fit_cph, centered = TRUE) plot(h0_naive$time, h0_naive$hazard, type = "s", xlab = "Time (years)", ylab = "Cumulative baseline hazard", main = "Weighted vs. unweighted baseline hazard") lines(h0_gw$time, h0_gw$hazard, type = "s", col = "steelblue") legend("topleft", c("Unweighted (cph)", "Weighted (svycoxph)"), col = c("black", "steelblue"), lty = 1) ``` ```{r substitute-basehaz} # Substitute Greenwood-weighted hazard for survplot() with visible bands fit_fused <- svycph_set_basehaz(fit_fused, h0_gw) ``` ```{r test-survplot} survplot(fit_fused, GGT = c(20, 50, 100), conf = "bands", xlab = "Follow-up (years)", ylab = "Survival", label.curves = list(keys = "lines")) ``` --- ## Step 5: Degrees of freedom correction `anova.rms()` uses the rank of the contrast matrix for degrees of freedom. For proper survey inference, F-test df should be based on the number of PSUs minus the number of strata. This section documents the correction. ```{r survey-df} # svycoxph stores degf.resid = n_PSU - n_strata directly; no manual computation needed fit_svy$degf.resid # e.g. 138 for the NHANES 1999-2018 analysis population fit_fused$svycph_vcov_df # same value, copied into fused object by svycph_fuse() ``` `anova.rms()` reports chi-square statistics (df = rank of contrast matrix), not F-statistics. The survey df (138) therefore affects interpretation rather than the test statistic itself: very large chi-squares remain informative, but for borderline results `regTermTest()` should be used as a check, as it denominates using survey df and returns an F-statistic with correct finite- population correction. ```{r regTermTest-comparison} # regTermTest() as a check on borderline spline nonlinearity results regTermTest(fit_svy, ~ rcs(GGT, 4)) # overall GGT association regTermTest(fit_svy, ~ rcs(LBXSAL, 4)) # overall albumin association ``` --- ## Conclusions and current limitations ### 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: | Function | Status | Notes | |---|---|---| | `anova.rms()` | **Working** | Survey-correct Wald tests with nonlinearity decomposition | | `Predict()` | **Working** | Survey-correct effect displays with CIs | | `plot.Predict()` / `ggplot.Predict()` | **Working** | Inherits from `Predict()` | | `summary.rms()` | **Working** | Inherits from correct `$var` and `$Design` | | `survplot()` | **Working** | Requires `svycph_set_basehaz()` first | ### 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: | Predictor | χ² (fused) | χ² (cph) | Ratio | |---|---|---|---| | GGT (overall) | 144.2 | 227.9 | 1.58 | | GGT (nonlinear) | 58.7 | 90.9 | 1.55 | | Albumin (overall) | 162.3 | 251.2 | 1.55 | | Albumin (nonlinear) | 25.6 | 36.8 | 1.44 | | Age | 3407.7 | 5609.4 | 1.65 | | Sex | 106.0 | 123.1 | 1.16 | | Race/ethnicity | **8.6** | **34.6** | **4.02** | | BMI | 6.6 | 25.8 | 3.90 | 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: ```{r icc-screen} 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) ``` ```{r deff-screen} # 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`. --- ## Session information ```{r session-info} sessionInfo() ```