Selective inference across NHANES cycles: lmFScreen example
Olivia McGough
2025-07-10
lmFScreen.RmdThis analysis is meant to mimic the analysis in the paper “Decreasing Trend of Bone Mineral Density in US Multiethnic Population: Analysis of Continuous NHANES 2005 - 2014”. The results of this analysis can be found in the 2025 paper “Valid F-screening in linear regression” by McGough, Witten, and Kessler (arxiv preprint: https://arxiv.org/abs/2505.23113).
Introduction
In this analysis we use the lmFScreen package to compare
femur neck bone mineral density (BMD) across National Health and
Nutrition Examination Survey (NHANES) cycles for men over age 30. NHANES
is a biennial survey that reports a series of health metrics in a large
sample of people in the United States. To conduct this analysis, we
first perform an F-test of the “overall hypothesis of no change in BMD
across cycles. This overall null hypothesis is rejected, so we
subsequently conduct post hoc tests for the individual coefficients that
correspond to each year/cycle in the data. We adjust for the fact that
we would not have conducted the post-hoc tests if we had not rejected
the overall null hypothesis.
Load and Merge NHANES Data
We load both BMD and demographic data from four NHANES cycles and
merge by participant ID (SEQN).
read_nhanes_cycle <- function(year) {
BMD <- read_xpt(paste0("DXXFEM_", year, ".xpt"))
DEMO <- read_xpt(paste0("DEMO_", year, ".xpt"))
merge(BMD, DEMO, by = "SEQN")
}
merged_2005_2006 <- read_nhanes_cycle("2005_2006")
merged_2007_2008 <- read_nhanes_cycle("2007_2008")
merged_2009_2010 <- read_nhanes_cycle("2009_2010")
merged_2013_2014 <- read_nhanes_cycle("2013_2014")Clean and Prepare Data
We retain only age, gender, and femur neck BMD, rename the columns, remove missing values, and restrict to men older than 30.
prep_bmd_data <- function(df) {
df <- df[, c("RIDAGEYR", "DXXNKBMD", "RIAGENDR")]
colnames(df) <- c("Age", "Femur_Neck_BMD", "Gender")
df <- na.omit(df)
df %>% filter(Age > 30)
}
bmd_men <- bind_rows(
prep_bmd_data(merged_2005_2006) %>% filter(Gender == 1) %>% mutate(Year = "2005-2006"),
prep_bmd_data(merged_2007_2008) %>% filter(Gender == 1) %>% mutate(Year = "2007-2008"),
prep_bmd_data(merged_2009_2010) %>% filter(Gender == 1) %>% mutate(Year = "2009-2010"),
prep_bmd_data(merged_2013_2014) %>% filter(Gender == 1) %>% mutate(Year = "2013-2014")
)Construct Design Matrix
We create indicator variables for each NHANES cycle.
bmd_men <- bmd_men %>% mutate(
Year = factor(Year, levels = c("2005-2006", "2007-2008", "2009-2010", "2013-2014")),
Year_2005_2006 = as.integer(Year == "2005-2006"),
Year_2007_2008 = as.integer(Year == "2007-2008"),
Year_2009_2010 = as.integer(Year == "2009-2010"),
Year_2013_2014 = as.integer(Year == "2013-2014")
)Conduct overall F-test
We run an overall F-test to check if there are significant differences in BMD across the years.
##
## Call:
## lm(formula = Femur_Neck_BMD ~ Year + Age, data = bmd_men)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37743 -0.09450 -0.00985 0.07919 0.70088
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.039434 0.007076 146.891 < 2e-16 ***
## Year2007-2008 -0.003687 0.004599 -0.802 0.423
## Year2009-2010 -0.003280 0.004567 -0.718 0.473
## Year2013-2014 -0.019821 0.004900 -4.045 5.29e-05 ***
## Age -0.003448 0.000112 -30.801 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1349 on 7130 degrees of freedom
## Multiple R-squared: 0.1235, Adjusted R-squared: 0.123
## F-statistic: 251.1 on 4 and 7130 DF, p-value: < 2.2e-16The overall F-statistic is significant, so we proceed with inference on the individual coefficients for each year/cycle.
Prepare for Pairwise Year Comparisons
We now define a function that:
Constructs a contrast
Projects out age and the intercept
Runs
lmFScreen.fiton the resulting matrix, which conducts inference on the individual coefficients for each year in the model while accounting for the overall F-test.
do_year_comparison <- function(X, test_col, seed_val) {
set.seed(seed_val)
n <- nrow(X)
Y <- bmd_men$Femur_Neck_BMD
project_out <- cbind(bmd_men$Age)
svdP <- svd(project_out, nu = nrow(project_out))
U_perp <- svdP$u[, (sum(svdP$d > .Machine$double.eps) + 1):ncol(svdP$u)]
X <- t(U_perp) %*% X
Y <- t(U_perp) %*% Y
result <- lmFScreen.fit(X, Y, test_cols = test_col, alpha = 0.05, alpha_ov = 0.05)
summary(result)
}
do_year_comparison <- function(X, test_col, seed_val) {
set.seed(seed_val)
n <- nrow(X)
Y <- bmd_men$Femur_Neck_BMD
project_out <- cbind(bmd_men$Age, 1)
svdP <- svd(project_out, nu = nrow(project_out))
U_perp <- svdP$u[, (sum(svdP$d > .Machine$double.eps) + 1):ncol(svdP$u)]
X <- t(U_perp) %*% X
Y <- t(U_perp) %*% Y
result <- lmFScreen(Y~X+0, test_cols = test_col, alpha = 0.05, alpha_ov = 0.05)
summary(result)
}Run Pairwise Year Comparisons
# Compare 2005-2006 vs 2007-2008
do_year_comparison(X=cbind(bmd_men$Year_2005_2006 + bmd_men$Year_2007_2008,
bmd_men$Year_2007_2008, bmd_men$Year_2009_2010),
test_col = 2, seed_val = 5678)## lmFScreen Model Summary
## --------------------------------------
## Overall F-statistic: 6.8067
## --------------------------------------
##
## Number of post hoc tests: 1
## --------------------------------------
##
## Selective Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X2 -0.003686 -0.0127 0.0053 0.4225
##
## Standard Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X2 -0.003686 -0.0127 0.0053 0.4228
##
##
## Significance levels: * < 0.05 ** < 0.01 *** < 0.001
# Compare 2005-2006 vs 2009-2010
do_year_comparison(X=cbind(bmd_men$Year_2005_2006 + bmd_men$Year_2009_2010,
bmd_men$Year_2007_2008, bmd_men$Year_2009_2010),
test_col = 3, seed_val = 56910)## lmFScreen Model Summary
## --------------------------------------
## Overall F-statistic: 6.8067
## --------------------------------------
##
## Number of post hoc tests: 1
## --------------------------------------
##
## Selective Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.003280 -0.0123 0.0057 0.4731
##
## Standard Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.003280 -0.0122 0.0057 0.4727
##
##
## Significance levels: * < 0.05 ** < 0.01 *** < 0.001
# Compare 2005-2006 vs 2013-2014
do_year_comparison(X=cbind(bmd_men$Year_2005_2006 + bmd_men$Year_2013_2014,
bmd_men$Year_2009_2010,
bmd_men$Year_2013_2014),
test_col = 3, seed_val = 561314)## lmFScreen Model Summary
## --------------------------------------
## Overall F-statistic: 6.8067
## --------------------------------------
##
## Number of post hoc tests: 1
## --------------------------------------
##
## Selective Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.018638 -0.0294 -0.0105 0.0010**
##
## Standard Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.019821 -0.0294 -0.0102 0.0001***
##
##
## Significance levels: * < 0.05 ** < 0.01 *** < 0.001
# Compare 2007-2008 vs 2009-2010
do_year_comparison(X=cbind(bmd_men$Year_2005_2006,
bmd_men$Year_2007_2008+bmd_men$Year_2009_2010,
bmd_men$Year_2009_2010),
test_col = 3, seed_val = 78910)## lmFScreen Model Summary
## --------------------------------------
## Overall F-statistic: 6.8067
## --------------------------------------
##
## Number of post hoc tests: 1
## --------------------------------------
##
## Selective Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 0.000407 -0.0079 0.0087 0.9227
##
## Standard Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 0.000407 -0.0079 0.0087 0.9233
##
##
## Significance levels: * < 0.05 ** < 0.01 *** < 0.001
# Compare 2007-2008 vs 2013-2014
do_year_comparison(X=cbind(bmd_men$Year_2009_2010,
bmd_men$Year_2007_2008+bmd_men$Year_2013_2014,
bmd_men$Year_2013_2014),
test_col = 3, seed_val = 781314)## lmFScreen Model Summary
## --------------------------------------
## Overall F-statistic: 6.8067
## --------------------------------------
##
## Number of post hoc tests: 1
## --------------------------------------
##
## Selective Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.016130 -0.0251 -0.0071 0.0004***
##
## Standard Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.016134 -0.0251 -0.0072 0.0004***
##
##
## Significance levels: * < 0.05 ** < 0.01 *** < 0.001
# Compare 2009-2010 vs 2013-2014
do_year_comparison(X=cbind(bmd_men$Year_2005_2006,
bmd_men$Year_2009_2010+bmd_men$Year_2013_2014,
bmd_men$Year_2013_2014),
test_col = 3, seed_val = 9101314)## lmFScreen Model Summary
## --------------------------------------
## Overall F-statistic: 6.8067
## --------------------------------------
##
## Number of post hoc tests: 1
## --------------------------------------
##
## Selective Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.016300 -0.0255 -0.0085 0.0006***
##
## Standard Estimates:
## Predictor Estimate Lower.CI Upper.CI P-value
## -------------------------------------------------------------
## X3 -0.016541 -0.0255 -0.0076 0.0003***
##
##
## Significance levels: * < 0.05 ** < 0.01 *** < 0.001Interpretation
Our selective p-values, point estimates, and intervals are valid conditional on rejecting the overall null hypothesis, in contrast to standard (naive) linear regression output. However, in this particular case, the selective p-values, estimates, and confidence intervals tend to coincide with their standard counterparts, indicating that selection (i.e. the face that the post hoc tests were conducted conditional on rejecing the overall null hypothesos) did not have a large impact on the analysis.