library(tidyverse)
subgroups <- expand.grid(
age = c("18-35", "36-50", "51-65", "66-80", "80+"),
sex = c("Male", "Female"),
hand = c("Left", "Right"),
coffee = c("0", "1", "2+"),
alcohol = c("0", "1", "2", "3+"),
physical_activity = c("0", "1", "2", "3+")
)
set.seed(21)
n <- 100 # Sample size
# Simulate recruitment (y: outcome)
y <- rnorm(n)
# Simulate recruitment (assume each new person has random subgroup)
x_idx <- sample(1:nrow(subgroups), n, replace = TRUE)
X <- subgroups[x_idx,]
simdata <- tibble(y = y, X)
# Create model matrix w/all interactions
X2 <- model.matrix(y ~ . * ., data = simdata)[,-1]
# combine into data set for model fitting
simdata2 <- data.frame(y=y, X2)
NoteReviewed by Ryan Peterson on 2026-08-13
In a previous post, Dr. P. Hackman demonstrated the remarkable effectiveness of unadjusted post-selection inference (UPSI) for achieving publishable, significant results. By searching for interesting “subgroups” through a large collection of main effects and interactions and then applying ordinary inferences as though no selection had occurred, Hackman virtually guaranteed statistically significant findings even when none existed.
The statistical problem I want to focus on here is not about causation: it’s about what happens when we use the data to select which subgroups to test, and then act as if those subgroups had been specified in advance. This issue applies equally to randomized trials and observational studies. Randomization protects against confounding, but not against the false positives that come from data-driven model selection. The analysis in the original post is simply performing model selection and then pretending the selected model had been specified in advance.
While this strategy may be excellent for a fictional statistician seeking ill-gotten fame and glory, sane analysts are interested in a slightly different goal: obtaining inference that is actually valid. This raises a natural question: if I select a model, and UPSIs are invalid, what should I use instead? In this post, I’ll revisit Hackman’s simulation setup and illustrate a more rigorous approach to post-selection inference: selective inference.
TipThe
selectInferToolkit R package
Selective inference, forward stepwise selection, and many other selection and inference methods are available in the selectInferToolkit R package, available on GitHub.
UPSI Alternatives
One proposed remedy, discussed in Maggie Qian’s follow-up post, The Antidote to Guaranteed Significance, is multiplicity adjustment. Qian shows that a Bonferroni or similar correction adjusts p-values for the number of tests performed, and therefore substantially “brings the false positive rate back to earth.” However, even Bonferroni, conservative as it is, doesn’t bring it down to the nominal \(\alpha=0.05\) level; the null still gets rejected nearly three times the nominal rate. A key reason is that multiplicity adjustment requires knowing how many tests you’re adjusting for. In the model selection setting, that number is genuinely ambiguous. Is it the \(2^p\) possible models? Something proportional to \(p\)? The number of steps taken during the forward search? The answer lies somewhere in between, and no single correction cleanly resolves the problem.
More importantly, the subgroup discovery framing is a special case of a more general problem. In many analyses, there is no pre-specified treatment and no enumerable list of candidate subgroups: we simply have a large collection of predictors and interactions, and we use the data to select a model. In the general setting, model selection (usually, but not always!) presumes no natural pre-specified hierarchy and no obvious grouping structure. The post-selection inference problem remains, and we need methods that address it directly.
Selective inference is designed for exactly this situation. Rather than counting comparisons after the fact, it directly accounts for the model selection procedure by asking: given that this particular model was selected from data, what can we still infer? In effect, it conditions on the selection event and recalibrates the usual p-values and confidence intervals accordingly. The price of this intellectual humility can be quite steep: the resulting intervals are often wider - sometimes infinite in length. Worse, the impressively small p-values driven by Hackman’s underhanded UPSIs will, more often than not, lose their sheen of significance.
TipHow does selective inference actually work though?
For readers interested in more theoretical properties of selective inference can refer to slides here.
Example: a study for treatment of chronic pain
Hackman’s UPSI approach
Let’s revisit the (rather extreme) simulation in the original post which uses interactions between 6 covariates to define subgroups:
As we saw in Hackman’s original post, when using forward stepwise selection with BIC, every single selected subgroup was declared statistically significant as reproduced below. Importantly, Hackman’s UPSI approach finds highly significant results despite the relationship between treatment and response being completely random!
library(selectInferToolkit) #github.com/petersonR/selectInferToolkit
# Fit the stepwise BIC model
fit_bic <- select_stepwise_ic(y ~ ., data = simdata2,
direction = "forward",
penalty = "BIC")
# Perform unadjusted post-selection inference (UPSI)
infer_upsi(fit_bic, data = simdata2) %>%
tidy() %>%
filter(coef != 0) | term | coef | ci_low | ci_high | p_value |
|---|---|---|---|---|
| (Intercept) | 0.07 | −0.09 | 0.24 | 0.393 |
| age36.50 | 0.27 | 0.10 | 0.45 | 0.003 |
| age66.80 | −0.38 | −0.57 | −0.18 | <0.001 |
| age66.80.coffee1 | 0.23 | 0.04 | 0.42 | 0.019 |
| age51.65.coffee2. | 0.29 | 0.11 | 0.46 | 0.002 |
| age51.65.physical_activity1 | −0.24 | −0.42 | −0.07 | 0.007 |
| sexFemale.physical_activity2 | −0.32 | −0.49 | −0.15 | <0.001 |
Presenting UPSI-based results like this can be deeply misleading. Imagine a research team applies a similar approach to an observational study of pain management, discovers that patients aged 36-50 show a “significantly” stronger association with better outcomes, and publishes accordingly. A clinical lab takes the finding at face value and designs a confirmatory trial targeting that age group. The trial fails - not because the treatment doesn’t work, but because the original association was an artifact of searching through dozens of subgroups without accounting for the search. Time, money, and patient participation were spent chasing a signal that was never there.
CautionUPSIs are “oopsies”
The fact that a variable was selected and appears statistically significant in the unadjusted post-selection inferences of the final model (e.g. a significant UPSI-based p-value) does not imply that the effect is real. Much of the apparent certainty can disappear once we account for the model selection process itself.
Selective inference to the rescue
When performing inference after selection, we know that since we’ve used the data both to select the model and now are using the same data to perform inference, these results are overly-optimistic - they are UPSIs. So next let’s use selective inference, which explicitly accounts for the fact that the model itself was selected using the data.
# infer_selective can perform selective inference
infer_selective(fit_bic, data = simdata2) %>%
tidy() %>%
filter(coef != 0) | term | coef | ci_low | ci_high | p_value |
|---|---|---|---|---|
| (Intercept) | 0.07 | NA | NA | NA |
| age36.50 | 0.27 | −9.26 | 0.33 | 0.879 |
| age66.80 | −0.38 | −Inf | 0.72 | 0.080 |
| age66.80.coffee1 | 0.23 | −1.30 | 1.63 | 0.381 |
| age51.65.coffee2. | 0.29 | −5.85 | 1.21 | 0.759 |
| age51.65.physical_activity1 | −0.24 | −1.57 | 1.06 | 0.333 |
| sexFemale.physical_activity2 | −0.32 | −6.74 | 0.61 | 0.099 |
Here we see that after accounting for the selection process, none of the selected effects remain statistically significant! Let’s take a moment to appreciate this victory over Hackman’s UPSI approach.
Of course, as mentioned previously, one caveat is confidence intervals become substantially wider. In this example, the lower CI bound for the 66-80 age subgroup extends to \(-\infty\). This isn’t an error or numerical instability. Rather, it reflects the uncertainty introduced by model selection: when an effect is identified only after searching through many possible models, accounting for that selection process can result in substantially wider confidence intervals.
A victory? Or a fluke?
I can almost hear Hackman coughing back: this was only a single simulation, I might have lost that battle but surely not the war. Even Bonferroni and his ridiculous adjustment failed to fully rein in my UPSIs.
Selective inference did successfully wash the false sheen of significance off of Hackman’s UPSI-based p-values once. Taking a page from Hackman himself, let’s repeat this 1000 times, as though we run the study in 1000 parallel universes to confirm without any doubt that selective inference controls the false-positive results. I will also add the simple antidote with Bonferroni multiplicity correction (with the number of tests set to \(p\) as in Qian’s post) to see how it performs.
As noted earlier, selective inference controls type I error per selected variable conditional on selection, but not across all selected variables within a single analysis/dataset. When BIC selects multiple variables and each is tested at \(\alpha=0.05\), the chance of at least one false positive per simulation compounds. To address this, I also apply a within-simulation Bonferroni correction on top of selective inference, adjusting each SI p-value by the number of selected variables in that simulation. This adds a second layer of multiplicity control: SI accounts for the search process, while the Bonferroni step accounts for testing multiple selected variables simultaneously.
# create list to store different results
sim_results_UPSI <- list()
sim_results_SI <- list() # selective inference only
sim_results_SI_BF <- list() # selective inference with bonferroni correction
sim_results_BONF <- list()
set.seed(1234)
## The below takes about 10 minutes to run
for(s in 1:1000) {
simdata2$y <- rnorm(n)
fit_bic <- select_stepwise_ic(
y ~ ., data = simdata2,
direction = "forward",
penalty = "BIC"
)
sim_results_UPSI[[s]] <-
infer_upsi(fit_bic, data = simdata2) %>%
tidy() %>%
filter(coef != 0)
sim_results_BONF[[s]] <-
infer_upsi(fit_bic, data = simdata2) %>%
tidy() %>%
filter(coef != 0) %>%
mutate(
n_tests = ncol(simdata2) - 1, # number of predictors, p,
p_value = pmin(p_value * n_tests, 1) # Bonferroni adjusted
)
sim_results_SI[[s]] <-
infer_selective(fit_bic, data = simdata2) %>%
tidy() %>%
filter(coef != 0)
sim_results_SI_BF[[s]] <-
infer_selective(fit_bic, data = simdata2) %>%
tidy() %>%
filter(coef != 0) %>%
mutate(
n_selected = sum(term != "(Intercept)"), # number of selected non-intercept terms
p_value = pmin(p_value * n_selected, 1)
)
}
upsi_results <- bind_rows(sim_results_UPSI, .id = "sim")
si_results <- bind_rows(sim_results_SI, .id = "sim")
si_bonf_results <- bind_rows(sim_results_SI_BF, .id = "sim")
bonf_results <- bind_rows(sim_results_BONF, .id = "sim")I can now summarize the results as shown below.
Code
summarise_results <- function(results) {
results <- results %>% filter(term != "(Intercept)")
per_sim <- results %>%
group_by(sim) %>% # simulation level summary
summarise(
n_disc = sum(!is.na(term)),
n_sig = sum(p_value < 0.05, na.rm = TRUE),
any_sig = n_sig > 0,
.groups = "drop"
)
per_sim %>%
summarise(
n_sims = n(), # sims with any selections
avg_disc = mean(n_disc), # avg selections per sim ( should be same)
avg_n_sig = mean(n_sig), # avg significant finding per simulation
pct_any_sim = mean(any_sig) * 100, # % sims with >=1 false positive
pct_sig_var = sum(n_sig)/sum(n_disc) * 100 # % of selected vars signif.
)
}
upsi_summary <-summarise_results(upsi_results)
si_summary <- summarise_results(si_results)
si_bonf_summary <- summarise_results(si_bonf_results)
bonf_summary <- summarise_results(bonf_results)Code
tibble(
Method = c("UPSI", "Bonferroni (BF)",
"Selective Inference (SI)", "SI + BF correction"),
`Avg selected` = c(upsi_summary$avg_disc, bonf_summary$avg_disc,
si_summary$avg_disc, si_bonf_summary$avg_disc),
`% vars sig` = c(upsi_summary$pct_sig_var, bonf_summary$pct_sig_var,
si_summary$pct_sig_var, si_bonf_summary$pct_sig_var),
`% sims ≥1 FP` = c(upsi_summary$pct_any_sim, bonf_summary$pct_any_sim,
si_summary$pct_any_sim, si_bonf_summary$pct_any_sim),
`Avg FP/sim` = c(upsi_summary$avg_n_sig, bonf_summary$avg_n_sig,
si_summary$avg_n_sig, si_bonf_summary$avg_n_sig)
)| Method | Avg selected | % vars sig | % sims ≥1 FP | Avg FP/sim |
|---|---|---|---|---|
| UPSI | 3.12 | 98.70 | 100.00 | 3.08 |
| Bonferroni (BF) | 3.12 | 14.31 | 28.10 | 0.45 |
| Selective Inference (SI) | 3.12 | 4.71 | 13.72 | 0.15 |
| SI + BF correction | 3.12 | 1.12 | 3.40 | 0.04 |
Across 1000 simulated null datasets, BIC selected on average of about 3 variables per simulation. The table above summarizes false positive rates across four approaches.
UPSI declared 99% of selected variables significant and produced at least one false positive in every simulation, averaging 3 spurious significant discoveries per simulation. In other words, by combining model selection with ordinary inference, we (like Hackman) were able to manufacture statistically significant findings 100% of the time.
A Bonferroni correction applied to the UPSI p-values, adjusting for all 92 candidate predictors, reduced the variable-level false positive rate to 14% and the simulation-level rate to about 28%. This is a meaningful improvement, but still far above the nominal 5% level. Hackman can still conjure false significance all too easily. The core difficulty with a simple multiplicity adjustment like Bonferroni is that it requires knowing how many comparisons were made, and that number is genuinely ambiguous after model selection.
Selective inference brought the variable-level rate down to the 5% nominal rate. Yet still, 13.7% of simulations yielded a falsely significant finding! What gives? It turns out that SI controls type I error rate per selected variable conditional on selection, but when BIC selects \(\sim\) 3 variables per simulation, each is tested at 5%, and thus the probability of at least one false positive per simulation roughly becomes \(1 - (0.95)^3 \sim 14\%\), even with perfect variable-level type I error control. Adding a Bonferroni correction within each simulation, that is, adjusting for the number of selected variables tested within each simulation, brings the simulation-level family-wise error rate down at last to ~5%.
Final thoughts/Conclusion
The original post framed the analysis as subgroup discovery, with a pre-specified treatment as the primary hypothesis. But the post-selection inference problem arises in any setting where the data are used to select a model, whether that’s using lasso to build a prediction model from observational data, screening hundreds of genes to identify relevant biomarkers, or an analyst cycling through several covariate sets before settling on a final specification. In none of these cases is there a pre-specified treatment hypothesis driving the analysis, and in all of them ordinary inference will be overly optimistic after selection.
Here I’ve only evaluated selective inference in terms of false positive control. While selective inference can address the false positive rate, it comes at the cost of very wide CIs (sometimes infinite!). In a future post, we’ll explore how powerful selective inference is when true effects actually exist, as well as how precise are CI obtained from selective inference.
Key Takeaways
Performing inference after variable selection using the UPSI approach nearly guarantees false positives.
Multiplicity adjustment like Bonferroni can’t fix the post-selection inference problem alone.
Selective inference can control false-positive at nominal alpha level. Our new R package, the
selectInferToolkitpackage leverages theselectiveInferencepackage to render adjusting inferences for selection highly practical.
Future Threads
What are alternative post-selection inference methods?
How do other post-selection inference methods (i.e., bootstrap, PIPE) implemented in the
selectInferToolkitpackage compare in terms of false discoveries vs power?Can we still make some inference about non-selections?
Appendix
R Session Info
sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.2 (2025-10-31)
os macOS Tahoe 26.6.2
system aarch64, darwin20
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/Chicago
date 2026-09-08
pandoc 3.9.0.2 @ /opt/homebrew/bin/ (via rmarkdown)
quarto 1.9.36 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
adaptMCMC 1.5 2024-01-29 [1] CRAN (R 4.5.0)
backports 1.5.1 2026-04-03 [1] CRAN (R 4.5.2)
broom * 1.0.13 2026-05-14 [1] CRAN (R 4.5.2)
class 7.3-23 2025-01-01 [2] CRAN (R 4.5.2)
cli 3.6.6 2026-04-09 [1] CRAN (R 4.5.2)
coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.5.0)
codetools 0.2-20 2024-03-31 [2] CRAN (R 4.5.2)
data.table 1.18.4 2026-05-06 [1] CRAN (R 4.5.2)
digest 0.6.39 2025-11-19 [1] CRAN (R 4.5.2)
dplyr * 1.2.1 2026-04-03 [1] CRAN (R 4.5.2)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.5.0)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.5.0)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.5.0)
forcats * 1.0.1 2025-09-25 [1] CRAN (R 4.5.0)
foreach 1.5.2 2022-02-02 [1] CRAN (R 4.5.0)
fs 2.1.0 2026-04-18 [1] CRAN (R 4.5.2)
future 1.75.0 2026-07-20 [1] CRAN (R 4.5.2)
future.apply 1.20.2 2026-02-20 [1] CRAN (R 4.5.2)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.5.0)
ggplot2 * 4.0.3 2026-04-22 [1] CRAN (R 4.5.2)
glmnet 5.0 2026-05-04 [1] CRAN (R 4.5.2)
globals 0.19.1 2026-03-13 [1] CRAN (R 4.5.2)
glue 1.8.1 2026-04-17 [1] CRAN (R 4.5.2)
gower 1.0.2 2024-12-17 [1] CRAN (R 4.5.0)
gt * 1.3.0 2026-01-22 [1] CRAN (R 4.5.2)
gtable 0.3.6 2024-10-25 [1] CRAN (R 4.5.0)
hardhat 1.4.3 2026-04-04 [1] CRAN (R 4.5.2)
here * 1.0.2 2025-09-15 [1] CRAN (R 4.5.0)
hms 1.1.4 2025-10-17 [1] CRAN (R 4.5.0)
htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.5.2)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.5.0)
intervals 0.15.5 2024-08-23 [1] CRAN (R 4.5.0)
ipred 0.9-15 2024-07-18 [1] CRAN (R 4.5.0)
iterators 1.0.14 2022-02-05 [1] CRAN (R 4.5.0)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.5.0)
knitr 1.51 2025-12-20 [1] CRAN (R 4.5.2)
lattice 0.22-9 2026-02-09 [1] CRAN (R 4.5.2)
lava 1.9.2 2026-06-30 [1] CRAN (R 4.5.2)
lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.5.2)
listenv 1.0.0 2026-06-22 [1] CRAN (R 4.5.2)
lubridate * 1.9.5 2026-02-04 [1] CRAN (R 4.5.2)
magrittr 2.0.5 2026-04-04 [1] CRAN (R 4.5.2)
MASS 7.3-65 2025-02-28 [2] CRAN (R 4.5.2)
Matrix 1.7-5 2026-03-21 [1] CRAN (R 4.5.2)
ncvreg 3.16.0 2025-10-09 [1] Github (pbreheny/ncvreg@5fecc8c)
nnet 7.3-20 2025-01-01 [1] CRAN (R 4.5.0)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.5.0)
parallelly 1.48.0 2026-06-29 [1] CRAN (R 4.5.2)
pbapply 1.7-4 2025-07-20 [1] CRAN (R 4.5.0)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.5.0)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.5.0)
prodlim 2026.03.11 2026-03-11 [1] CRAN (R 4.5.2)
purrr * 1.2.2 2026-04-10 [1] CRAN (R 4.5.2)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.5.0)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.5.0)
Rcpp 1.1.2 2026-07-05 [1] CRAN (R 4.5.2)
readr * 2.2.0 2026-02-19 [1] CRAN (R 4.5.2)
recipes 1.3.3 2026-05-30 [1] CRAN (R 4.5.2)
rlang 1.3.0 2026-07-05 [1] CRAN (R 4.5.2)
rmarkdown 2.31 2026-03-26 [1] CRAN (R 4.5.2)
rpart 4.1.27 2026-03-27 [1] CRAN (R 4.5.2)
rprojroot 2.1.1 2025-08-26 [1] CRAN (R 4.5.0)
rstudioapi 0.18.0 2026-01-16 [1] CRAN (R 4.5.2)
S7 0.2.2 2026-04-22 [1] CRAN (R 4.5.2)
sass 0.4.10 2025-04-11 [1] CRAN (R 4.5.0)
scales 1.4.0 2025-04-24 [1] CRAN (R 4.5.0)
selectInferToolkit * 0.4.4 2026-08-13 [1] Github (petersonR/selectInferToolkit@d41796d)
selectiveInference 1.2.5 2019-09-07 [1] CRAN (R 4.5.0)
sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.5.0)
shape 1.4.6.1 2024-02-23 [1] CRAN (R 4.5.0)
sparsevctrs 0.3.6 2026-01-27 [1] CRAN (R 4.5.2)
stringi 1.8.9 2026-08-04 [1] CRAN (R 4.5.2)
stringr * 1.6.0 2025-11-04 [1] CRAN (R 4.5.0)
survival 3.8-6 2026-01-16 [1] CRAN (R 4.5.2)
tibble * 3.3.1 2026-01-11 [1] CRAN (R 4.5.2)
tidyr * 1.3.2 2025-12-19 [1] CRAN (R 4.5.2)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.5.0)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.5.0)
timechange 0.4.0 2026-01-29 [1] CRAN (R 4.5.2)
timeDate 4052.112 2026-01-28 [1] CRAN (R 4.5.2)
tzdb 0.5.0 2025-03-15 [1] CRAN (R 4.5.0)
vctrs 0.7.3 2026-04-11 [1] CRAN (R 4.5.2)
withr 3.0.3 2026-06-19 [1] CRAN (R 4.5.2)
xfun 0.57 2026-03-20 [1] CRAN (R 4.5.2)
xml2 1.5.2 2026-01-17 [1] CRAN (R 4.5.2)
yaml 2.3.12 2025-12-10 [1] CRAN (R 4.5.2)
[1] /Users/rpterson/Library/R/arm64/4.5/library
[2] /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────Citation
For attribution, please cite this work as:
Shah, Jinal. 2026. “I Selected a Model. Can I Trust My p-Values
Yet?” Data Diction (blog). September 8, 2026. https://doi.org/10.59350/pcr99-e8q06.
