Function on Scalar Regression in Complex Survey Settings

Lily Koffman

UVA School of Data Science

Motivation: functional data in large epidemiological studies

National Health and Nutrition Examination Survey (NHANES)

NHANES

Accelerometer-measured physical activity

Accelerometer-measured physical activity

Accelerometer-measured physical activity

  • Does activity depend on covariates?

  • Are associations different at different times of the day?

Accelerometer-measured physical activity

Function on scalar regression (FoSR)

Model: \[\mathbb{E}[\mathrm{PA}_i(s)] = \beta_0(s) + \beta_1(s)\mathbb{I}\{\mathrm{sex}_i=\mathrm{male}\} + \sum_{k=2}^{7}\beta_k(s)\mathbb{I}\{\mathrm{agecat}_i=k\} \]

\(i\): participant;

\(s \in \{1, \dots, 1440\}\): each minute of the day

FoSR: model fitting

  • Fast univariate inference (FUI) (Cui et al. 2021)

  • Fit separate (univariate) GLM at each point \(s\), smooth the resulting point estimates

  • Bootstrap to get confidence bands

  • Other methods for FoSR exist: benefits of FUI in this application will become clear shortly

NHANES is not a random sample!

Individuals sampled in geographic clusters, minority groups oversampled

Are estimates from FoSR valid for population-level inference?

NHANES Sample Design Tutorial

NHANES is not a random sample!

  • Ignoring the sampling structure can lead to:

    • Biased estimation (need to account for weights to make survey look like population)
    • Incorrect inference (individuals in same stratum/PSU have correlated data; failure to account can under/overestimate variance)

  • No method exists for function on scalar regression under complex surveys

Survey terminology

  • Survey weight: number of people in the population represented by the sampled person

    • \(\text{weight} \propto \frac{1}{\text{probability of selection}}\)
    • In complex surveys like NHANES, weights also account for complex, multistage design, nonresponse, and post-stratification to match US population1

  • Stratum: larger geographic location

  • Cluster or Primary Sampling Unit (PSU): smaller area within stratum, such as zip code or city block

Proposed solution: let’s combine some existing methods

Survey function on scalar regression

  • GLMs: methods and software exist to take into account weights (estimation) and correlation between clusters for inference (\(\texttt{svyglm}\), \(\texttt{svycoxph}\)) (Lumley 2010)

  • FUI is built on separate GLMs

  • Idea: incorporate survey weights into the GLMs (estimation) and use survey-aware replication/bootstrap methods for inference

Implementation

Fit survey-weighted GLM at each point along the functional domain

Fit survey-weighted GLM at each point along the functional domain

Point estimates

Smoothed point estimates

What about inference?

  • Standard errors from weighted GLMs don’t account for correlation between individuals in same stratum/PSU
  • Survey-aware replication methods generate \(B\) sets of survey weights (analogous to bootstrap replicates), accounting for strata/PSU structure
  • Fit a model on each replicate and estimate standard errors from set of replicates
  • Several methods exist, we focus on balanced repeated replication (BRR) (Rao and Shao 1999)

Inference: BRR

Choose one PSU from each stratum, double weights for each individual in selected PSU, fit model

Inference: BRR

\(B\) coefficient function estimates

Inference: BRR

From \(B\) BRR replicates, get pointwise variance

Inference: BRR

Use pointwise variance to get pointwise confidence intervals

Joint inference

  • Pointwise intervals ask “are females more active than males at 8AM?”

  • Might also want a confidence band for the true function across entire functional domain

  • We have a multiplicity problem if we use the pointwise intervals: we’re doing \(1440\) tests

Joint inference: Bonferroni correction

Bonferroni: \(\alpha = \frac{0.05}{1440}; z_{1-\alpha/2}=3.98\)

Joint inference: Bonferroni correction

  • Bonferroni gives valid joint coverage, but it is conservative because tests are not independent
    • \(\hat{\beta}(8:00\text{am})\) and \(\hat{\beta}(8:01\text{am})\) are highly correlated
  • How can we adjust our CIs to take into account correlation \(and\) multiplicity (CMA)? (Crainiceanu et al. 2024)

Joint inference: intuition

  • Want a \(95\)% confidence band that covers the true \(\beta(s)\) everywhere simultaneously

  • We could think of a z-score at each time point: \(\frac{\hat{\beta}(s) - \beta(s)}{SE_{\beta(s)}}\)

  • We want to find a critical value \(q\) where probability that all z-scores are less than \(q\) in absolute value is \(1-\alpha\)

  • Equivalent to the max absolute z-score being less than \(q\)

  • Value of \(q\) depends on correlation between coefficients across domain

Joint inference: CMA

Step 1: recognize that coefficient estimates are multivariate Normal \[\hat{\boldsymbol{\beta}} = \big(\hat{\beta}(s_1), \hat{\beta}(s_2), \dots, \hat{\beta}(s_L)\big)^T \sim \mathcal{N}(\boldsymbol{\beta}, \mathbf{V}_\beta)\]

Joint inference: CMA

Step 1: recognize that coefficient estimates are multivariate Normal \[\hat{\boldsymbol{\beta}} = \big(\hat{\beta}(s_1), \hat{\beta}(s_2), \dots, \hat{\beta}(s_L)\big)^T \sim \mathcal{N}(\boldsymbol{\beta}, \mathbf{V}_\beta)\]

Step 2: Standardize to get z-scores

\[\mathbf{Z} = \frac{(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})}{\mathbf{D}_\beta} \sim \mathcal{N}(\mathbf{0}, \mathbf{C}_\beta)\] where \(\mathbf{D}_\beta = \sqrt{\text{diag}(\mathbf{V}_\beta)}\) and \(\mathbf{C}_\beta\) is the correlation matrix, in practice estimated via BRR replicates

Joint inference: CMA

Step 1: recognize that coefficient estimates are multivariate Normal \[\hat{\boldsymbol{\beta}} = \big(\hat{\beta}(s_1), \hat{\beta}(s_2), \dots, \hat{\beta}(s_L)\big)^T \sim \mathcal{N}(\boldsymbol{\beta}, \mathbf{V}_\beta)\]

Step 2: Standardize to get z-scores

\[\mathbf{Z} = \frac{(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})}{\mathbf{D}_\beta} \sim \mathcal{N}(\mathbf{0}, \mathbf{C}_\beta)\] where \(\mathbf{D}_\beta = \sqrt{\text{diag}(\mathbf{V}_\beta)}\) and \(\mathbf{C}_\beta\) is the correlation matrix, in practice estimated via BRR replicates

Step 3: find critical value \(q(\mathbf{C}_\beta, 1-\alpha)\) such that

\[P\big(\max_s|Z_s|\leq q\big) = 1-\alpha \text{, where } \mathbf{Z} \sim \mathcal{N}(\mathbf{0}, \mathbf{C}_\beta)\]

Joint inference: CMA

Bonferroni critical value = \(3.98\); CMA critical value = \(3.03\); pointwise critical value = \(1.96\)

Final result

Easy enough to implement – but how do we know if it works?

Compare survey-aware vs. unweighted approach with simulation

  • Generate a superpopulation of individuals who belong to a stratum and cluster
  • Each individual has functional outcome and covariate1
  • [Optionally] create correlation in functional outcomes within strata and clusters
  • [Optionally] create heterogeneity in covariate effects by strata/cluster
  • Sample from the superpopulation in each iteration of the simulation
    • [Optionally] allow for informative sampling: probability of being included in sample related to functional outcome
  • Evaluate coverage and estimation error (integrated squared error; ISE = \(\int_s \{\hat{\beta}(s) - \beta(s)\}^2ds\))

Simulated survey data

Survey sampling: sample individuals from superpopulation

Informative sampling: allow selection probability to depend on functional outcome: \(\pi_i \propto \frac{1}{L}\sum_{s=1}^LY_i(s)\)

Selected Simulation Results

Estimation error: intercept

Estimation error: coefficient

Coverage: intercept

Coverage: coefficient

Simulation summary

  • Unweighted only achieves adequate coverage/comparable estimation error w/ uniform sampling and no stratum/PSU random effects
  • Our method achieves adequate coverage and lower bias than unweighted methods under informative sampling and/or stratum/PSU effects
  • Results hold for Gaussian, Poisson, and binomial data
  • Empirical simulation resampling NHANES curves validated results

Assumptions

Survey design

  • Survey weights are bounded
  • Inclusion probabilities are positive and correctly specified
  • For BRR: two PSUs per stratum

Assumptions

Survey design

  • Survey weights are bounded
  • Inclusion probabilities are positive and correctly specified
  • For BRR: two PSUs per stratum

Model

  • \(\beta(s)\) is smooth
  • The mean is correctly specified: \(E(Y_i(s)|X_i) = g^{-1}(X_i\beta(s))\)
  • Functional outcomes \(Y_i(s)\) are independent conditional on covariates and survey design

Asymptotics

Consistency

  • Survey-weighted pointwise estimator \(\hat{\beta}(s)\) is consistent for the population parameter \(\beta(s)\)
  • Smoothed estimator \(\tilde{\beta}(s)\) is consistent for \(\beta(s)\) provided that smoothing parameter \(\lambda\) approaches zero at the appropriate rate

Asymptotics

Consistency

  • Survey-weighted pointwise estimator \(\hat{\beta}(s)\) is consistent for the population parameter \(\beta(s)\)
  • Smoothed estimator \(\tilde{\beta}(s)\) is consistent for \(\beta(s)\) provided that smoothing parameter \(\lambda\) approaches zero at the appropriate rate

Coverage

  • Pointwise CIs achieve nominal rate due to consistency of BRR variance estimators, provided smoothing bias of smaller order than standard error
  • Joint CMA intervals achieve nominal coverage under assumption that \(\hat{\beta}(s)\) is multivariate normal across the functional domain

Did we need to do all that survey stuff anyway?

Software and Computation

R Package

model_fit = svyfosr::svyfui(steps_mat ~ age + sex,
                            weights = survey_weight,
                            data = steps_df,
                            family = gaussian(),
                            boot_type = "BRR",
                            num_boots = 500,
                            parallel = TRUE,
                            nknots_min = 20,
                            seed = 2213)
Koffman and Muschelli (2025)

Computation

  • NHANES model (\(B = 500\), \(n = 12500\), \(L = 1440\)) takes under 4 minutes to run on standard computer
  • Speedups:
    • If predictors \(X\) are same across functional domain, instead of fitting a separate model at each time point, we can create \(n\times L\) matrix \(Y\) by stacking our outcomes and fit \(B = (X'WX)^{-1}X'WY\)
    • Use QR decomposition instead of calculating \((X'WX)^{-1}\) directly
    • For GLMs, apply above at each step of IRLS or Newton-Raphson
    • Parallelize bootstrapping

Summary, Limitations, Future Directions

Summary

  • Developed a method for function on scalar regression under complex survey designs, motivated by physical activity in NHANES
  • Demonstrated importance of approach in presence of informative sampling, within-PSU/strata covariate effects, and/or within-PSU/strata correlation
  • Computationally efficient software for implementation, support for Gaussian, Poisson, and binomial data
  • Paper on arXiv

Limitations and Future Directions

  • Multilevel data (e.g. multiple days)
  • Full asymptotic treatment
  • Smoothing order sensitivity

Appendix

Simulation data generation

Simulation data generation: add stratum-specific covariate effects

\(\beta_{1(h)}(s) = \gamma_h\cdot\beta_1(s)\)

Where \(\gamma_h \sim N(1, 0.125^2)\)

Simulation data generation: induce strata and cluster-level correlation with functional random effects

Generate strata-level (\(b_{h(i)}(s)\)) and cluster-level (\(b_{c(h, i)}(s)\)) random effects using B-spline basis

Individual-level linear predictor: \(\eta_i(s) = \beta_0(s) + X_i\beta_{1(h)}(s) + b_{h(i)}(s) + b_{c(h, i)}(s)\)

Simulation data generation: induce strata and cluster-level correlation with functional random effects

Generate strata-level (\(b_{h(i)}(s)\)) and cluster-level (\(b_{c(h, i)}(s)\)) random effects using B-spline basis

Individual-level linear predictor: \(\eta_i(s) = \beta_0(s) + X_i\beta_{1(h)}(s) + b_{h(i)}(s) + b_{c(h, i)}(s)\)

Simulation data generation: induce strata and cluster-level correlation with functional random effects

Generate strata-level (\(b_{h(i)}(s)\)) and cluster-level (\(b_{c(h, i)}(s)\)) random effects using B-spline basis

Individual-level linear predictor: \(\eta_i(s) = \beta_0(s) + X_i\beta_{1(h)}(s) + b_{h(i)}(s) + b_{c(h, i)}(s)\)

Simulation data generation: get functional outcomes from linear predictors

\(Y_i(s)\) sampled from \(\mathcal{N}(\eta_i(s), \sigma^2)\)

Smoothing order

  • Our approach is to fit pointwise models, then smooth coefficients
  • An alternative approach would be to smooth outcomes first, then fit pointwise models to the smoothed data
  • With Gaussian data and same linear smoother, approaches are equivalent
  • Difference is in how the smoothing parameter \(\lambda\) is selected
  • When smoothing applied to coefficients, whose pointwise variance is \(\sigma^2/n\), GCV selects smaller smoothing parameter as \(n\) grows, ensuring the estimator is consistent
  • When smoothing is applied to outcomes, pointwise variance is \(\sigma^2\), \(O(1)\) bias is introduced that doesn’t change with sample size
  • If the smooth applied to the outcomes is calibrated to have the same EDF as the smooth applied to the coefficients, comparable performance can be achieved between the two approaches
  • The next few slides illustrate this phenomenon

Smoothing order

Mean smoother EDF across 200 iterations
Varying sample size and smooth order
Sample size Mean EDF, smooth first Mean EDF, smooth second
500 5.14 11.28
1000 5.12 12.07
5000 5.13 14.22

Smoothing order

Smoothing order

Smoothing order

Smoothing order

References

Crainiceanu, C. M., J. Goldsmith, A. Leroux, and E. Cui. 2024. Functional Data Analysis with R. Springer New York, NY, USA.
Cui, Erjia, Andrew Leroux, Ekaterina Smirnova, and Ciprian M. Crainiceanu. 2021. “Fast Univariate Inference for Longitudinal Functional Models.” Journal of Computational and Graphical Statistics 31 (1): 219–30. https://doi.org/10.1080/10618600.2021.1950006.
Koffman, Lily, and John Muschelli. 2025. svyfosr: Survey-Weighted Function on Scalar Regression. https://github.com/jhuwit/svyfosr.
Lumley, Thomas. 2010. Complex Surveys: A Guide to Analysis Using R: A Guide to Analysis Using R. John Wiley; Sons.
Rao, J. N. K., and J. Shao. 1999. “Modified Balanced Repeated Replication for Complex Survey Data.” Biometrika 86 (2): 403–15. http://www.jstor.org/stable/2673521.