Add Bootstrap Wald Tests for 3-Step LCA Distal Outcomes#74
Open
CravingCrates wants to merge 2 commits into
Open
Add Bootstrap Wald Tests for 3-Step LCA Distal Outcomes#74CravingCrates wants to merge 2 commits into
CravingCrates wants to merge 2 commits into
Conversation
Introduce WaldTest3Step to perform non-parametric bootstrap inference for 3-step LCA structural models: computes bootstrap estimates, SEs, percentile CIs, pairwise Wald tests (with Bonferroni / BH corrections), and omnibus Wald tests. Add helper routines (_fdr_bh, _bonferroni, _build_contrast_matrix, _stars) and internal logic to build result tables and pretty-print summaries. Add a StepMix.wald_test convenience method to run the procedure from a fitted model. Include a comprehensive test suite (test/test_bootstrap_inference.py) covering helper functions, construction checks, continuous/binary/soft-assignment integrations, multiple-comparison corrections, and statistical sanity checks.
There was a problem hiding this comment.
Pull request overview
This PR adds a new bootstrap-based inference utility (WaldTest3Step) for 3-step StepMix distal outcome models, plus a StepMix.wald_test() convenience wrapper and a comprehensive test suite to validate the resulting estimates, confidence intervals, and Wald test outputs.
Changes:
- Added
WaldTest3Steptostepmix.bootstrapto compute bootstrapped SE/CI plus pairwise and omnibus Wald tests (with optional multiple-testing corrections). - Added
StepMix.wald_test()as a convenience API to run the bootstrap + return a fittedWaldTest3Step. - Added a new
test_bootstrap_inference.pymodule with unit tests for helpers and integration tests for continuous, binary, and soft-assignment workflows.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 9 comments.
| File | Description |
|---|---|
stepmix/bootstrap.py |
Adds Wald-test inference class + multiple-testing correction helpers built on top of existing bootstrap() output. |
stepmix/stepmix.py |
Adds wald_test() convenience method that instantiates and fits WaldTest3Step. |
test/test_bootstrap_inference.py |
Introduces a new (large) integration test suite for bootstrap inference + helper corrections. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| TestBinaryOutcome, | ||
| TestSoftAssignment, | ||
| TestCorrections, | ||
| TestConvenienceFunction, |
Comment on lines
+361
to
+362
| adjusted = np.full(n, np.nan) | ||
| adjusted[finite] = np.minimum(1.0, p_values[finite] * n) |
Comment on lines
+478
to
+479
| ``None`` → no correction; ``"Bonferroni"`` → Bonferroni-Holm; | ||
| ``"BH"`` → Benjamini-Hochberg FDR. |
Comment on lines
+328
to
+355
| def _fdr_bh(p_values: np.ndarray) -> np.ndarray: | ||
| """Benjamini-Hochberg FDR correction (two-stage). | ||
|
|
||
| Parameters | ||
| ---------- | ||
| p_values : 1-D array of raw p-values (NaN is preserved). | ||
|
|
||
| Returns | ||
| ------- | ||
| adjusted : 1-D array of BH-corrected p-values. | ||
| """ | ||
| n = len(p_values) | ||
| finite = np.isfinite(p_values) | ||
| adjusted = np.full(n, np.nan) | ||
| if finite.sum() == 0: | ||
| return adjusted | ||
|
|
||
| idx = np.where(finite)[0] | ||
| p = p_values[idx] | ||
| order = np.argsort(p) | ||
| ranks = np.empty_like(order) | ||
| ranks[order] = np.arange(1, len(p) + 1) | ||
| adj = np.minimum(1.0, p * len(p) / ranks) | ||
| # Enforce monotonicity (step-down) | ||
| for i in range(len(adj) - 2, -1, -1): | ||
| adj[order[i]] = min(adj[order[i]], adj[order[i + 1]]) | ||
| adjusted[idx] = adj | ||
| return adjusted |
Comment on lines
+578
to
+598
| def _get_boot_matrix(self, model_name: str, param: str, variable: str) -> np.ndarray: | ||
| """Return a (n_bootstrap × K) matrix of draws for a given variable. | ||
|
|
||
| Rows = bootstrap replications, columns = latent classes (0 … K-1). | ||
| """ | ||
| sm_df = self._get_structural_boot() | ||
| K = self.model.n_components | ||
| n_reps = sm_df["rep"].nunique() | ||
|
|
||
| mat = np.full((n_reps, K), np.nan) | ||
| for k in range(K): | ||
| mask = ( | ||
| (sm_df.index.get_level_values("model_name") == model_name) | ||
| & (sm_df.index.get_level_values("param") == param) | ||
| & (sm_df.index.get_level_values("class_no") == k) | ||
| & (sm_df.index.get_level_values("variable") == variable) | ||
| ) | ||
| sub = sm_df.loc[mask].sort_values("rep") | ||
| if len(sub) > 0: | ||
| mat[:, k] = sub["value"].values[:n_reps] | ||
| return mat # (B, K) |
Comment on lines
+1155
to
+1185
| def wald_test(self, X, Y=None, n_repetitions=500, ci_level=0.95, correction=None, progress_bar=True, random_state=None): | ||
| """Runs non-parametric bootstrap over the full 3-step procedure and returns Wald test results. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| X : array-like of shape (n_samples, n_features) | ||
| Measurement data (same as used to fit the model). | ||
| Y : array-like of shape (n_samples, n_outcome_features) | ||
| Outcome data (same as used to fit the structural model). | ||
| n_repetitions : int, default=500 | ||
| Number of bootstrap replications. >= 500 recommended. | ||
| ci_level : float, default=0.95 | ||
| Confidence level for bootstrap percentile CIs. E.g. 0.95 -> 95 % CI. | ||
| correction : {None, "Bonferroni", "BH"}, default=None | ||
| Multiple-comparison correction applied to pairwise p-values. | ||
| progress_bar : bool, default=True | ||
| Show a tqdm progress bar. | ||
| random_state : int, default=None | ||
| Random seed for the bootstrap. | ||
|
|
||
| Returns | ||
| ------- | ||
| wt : WaldTest3Step | ||
| A fitted object providing `.summary()`, `.pairwise_`, and `.omnibus_`. | ||
| """ | ||
| from stepmix.bootstrap import WaldTest3Step | ||
| wt = WaldTest3Step(self, ci_level=ci_level) | ||
| wt.fit_bootstrap( | ||
| X, Y, n_repetitions=n_repetitions, correction=correction, | ||
| progress_bar=progress_bar, random_state=random_state | ||
| ) |
Comment on lines
+1155
to
+1164
| def wald_test(self, X, Y=None, n_repetitions=500, ci_level=0.95, correction=None, progress_bar=True, random_state=None): | ||
| """Runs non-parametric bootstrap over the full 3-step procedure and returns Wald test results. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| X : array-like of shape (n_samples, n_features) | ||
| Measurement data (same as used to fit the model). | ||
| Y : array-like of shape (n_samples, n_outcome_features) | ||
| Outcome data (same as used to fit the structural model). | ||
| n_repetitions : int, default=500 |
Comment on lines
+1155
to
+1185
| def wald_test(self, X, Y=None, n_repetitions=500, ci_level=0.95, correction=None, progress_bar=True, random_state=None): | ||
| """Runs non-parametric bootstrap over the full 3-step procedure and returns Wald test results. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| X : array-like of shape (n_samples, n_features) | ||
| Measurement data (same as used to fit the model). | ||
| Y : array-like of shape (n_samples, n_outcome_features) | ||
| Outcome data (same as used to fit the structural model). | ||
| n_repetitions : int, default=500 | ||
| Number of bootstrap replications. >= 500 recommended. | ||
| ci_level : float, default=0.95 | ||
| Confidence level for bootstrap percentile CIs. E.g. 0.95 -> 95 % CI. | ||
| correction : {None, "Bonferroni", "BH"}, default=None | ||
| Multiple-comparison correction applied to pairwise p-values. | ||
| progress_bar : bool, default=True | ||
| Show a tqdm progress bar. | ||
| random_state : int, default=None | ||
| Random seed for the bootstrap. | ||
|
|
||
| Returns | ||
| ------- | ||
| wt : WaldTest3Step | ||
| A fitted object providing `.summary()`, `.pairwise_`, and `.omnibus_`. | ||
| """ | ||
| from stepmix.bootstrap import WaldTest3Step | ||
| wt = WaldTest3Step(self, ci_level=ci_level) | ||
| wt.fit_bootstrap( | ||
| X, Y, n_repetitions=n_repetitions, correction=correction, | ||
| progress_bar=progress_bar, random_state=random_state | ||
| ) |
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
Author
|
@copilot apply changes based on the comments in this thread |
Collaborator
|
Thanks! We'll merge if/when we find someone to review this thoroughly. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What this new code does
This PR introduces the
WaldTest3Stepclass to bootstrap.py, which provides a non-parametric bootstrap procedure to compute standard errors, confidence intervals, and Wald tests for 3-step Latent Class Analysis (LCA) distal outcome models.I love this package and use it a lot to do analysis, but one thing that bugged me is that the hypothesis testing is still left to the reader and many commercial statistic applications solve this issue, so I wanted to try and make a small contribution.
Feel free to deny it, if it's just too much work to review or doesn't fit into the scope of the project right now. No hard feelings!
Key features include:
summary()method to quickly view estimates, omnibus stats, and pairwise test results with significance stars.Example usage:
The output could look like this: