Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #133 +/- ##
==========================================
- Coverage 73.79% 73.61% -0.18%
==========================================
Files 26 20 -6
Lines 1408 1190 -218
==========================================
- Hits 1039 876 -163
+ Misses 369 314 -55
🚀 New features to boost your workflow:
|
|
Looks good so far - let me know if you need help with the regression for site correction in the prediction part |
for more information, see https://pre-commit.ci
|
For the site correction for the prediction metrics, I had to change the way we were doing the cross-validation. Instead of using the function from scikit-learn Let me know what you think of this approach. |
htwangtw
left a comment
htwangtw
left a comment
There was a problem hiding this comment.
A few things -
- The function forgot to fit each site separately.
- It would be clearer if we could create a customised function to fit into the sklearn pipeline function for cross-validation. There's
FunctionTransformerto create a sklearn-compatible object for the pipeline. This way you can also drop the duplication inregress_site
| beta, *_ = np.linalg.lstsq( | ||
| design_train, | ||
| X_train, | ||
| rcond=None, | ||
| ) | ||
|
|
||
| beta_site = beta[1:, :] |
There was a problem hiding this comment.
You need to fit the linear regression on the training and testing sets separately
beta_train, *_ = np.linalg.lstsq(design_train, X_train, rcond=None,)
beta_test, *_ = np.linalg.lstsq(design_test, X_test, rcond=None,)
X_train_corr = X_train - design_train[:, 1:] @ beta_train[1:, :]
X_test_corr = X_test - design_test[:, 1:] @ beta_test[1:, :]| X_train_corr = X_train - design_train[:, 1:] @ beta_site | ||
| X_test_corr = X_test - design_test[:, 1:] @ beta_site |
|
Thanks for the feedback! I added the fitting step for each train/test set separately. I agree that creating a scikit-learn–compatible pipeline object would make the implementation much cleaner. I'll refactor the code as soon as the current implementation is approved. |
|
@claraElk I think we can merge this so we have a record of a version that we know works. If the refactor doesn't pan out, we have something to fall back on. Do you feel confident implementing with a sklearn function transformer? |
|
I am working on it today. I'll let you know how it goes. So far I tried with FunctionTransformer, but I don't think I might work here (as we might need to extract the beta only from our training set, and use these values for the test set; see my question below). I am trying to create a class instead (based on @pbergeret12 code for another project). If I run into issues I'll let you know and we can merge this version instead. Example of the class: However, I have a question regarding the current implementation. At the moment, we are performing two separate linear regressions (one on the training set and one on the testing set). I was discussing this with @pbergeret12, and it seems that we should estimate the beta values only on the training set and then apply these same coefficients to the testing set to avoid information leakage. Would we need to modify the code before merging, for example: What do you think @htwangtw ? |
|
@pbergeret12, thank you for raising this point, and it got me to rethink my suggestion. Let's go with your suggestion now. Below is my train of thought: You are correct about data leakage. If we treat the site as a property the model should learn and generalise to correct, then it makes sense to only fit on the training set. (Side note - for this to work, we need proper stratification.) However, we are working with an issue that doesn't have a real solution. The subject-per-site distribution can be way too imbalanced in neuroimaging, so we cannot really properly stratify, and leaving one site out just doesn't work. This leaves us with one option that is not correct and will potentially introduce a positive bias due to data leakage: regress out the site as a nuisance variable on the whole dataset before CV. This is the same logic as why I suggested doing the nuisance regression on both - I was treating it like a timeseries-level confound, but this doesn't really work. If we implement @pbergeret12's solution, we might get user complaints that things are not working, etc., depending on the nature of the dataset. Implementing the second one will work more universally. It will be defensible but incorrect. And we will need to add a huge disclaimer to the documentation that no one would read. Might be worth bringing this up to Lune and seeing what she thinks. |
| def regress_site( | ||
| X_train: NDArray[np.float32], | ||
| X_test: NDArray[np.float32], | ||
| site_train: NDArray[np.str_], | ||
| site_test: NDArray[np.str_], | ||
| ) -> tuple[NDArray[np.float32], NDArray[np.float32]]: | ||
| """Regress out site effects from the training and test data.""" | ||
|
|
||
| train_site = pd.get_dummies(site_train, drop_first=True, dtype=float) | ||
| test_site = pd.get_dummies(site_test, drop_first=True, dtype=float) | ||
|
|
||
| test_site = test_site.reindex( | ||
| columns=train_site.columns, | ||
| fill_value=0.0, | ||
| ) | ||
|
|
||
| design_train = np.column_stack([np.ones(len(train_site)), train_site.to_numpy()]) | ||
|
|
||
| design_test = np.column_stack([np.ones(len(test_site)), test_site.to_numpy()]) | ||
|
|
||
| beta_train, *_ = np.linalg.lstsq( | ||
| design_train, | ||
| X_train, | ||
| rcond=None, | ||
| ) | ||
| beta_test, *_ = np.linalg.lstsq( | ||
| design_test, | ||
| X_test, | ||
| rcond=None, | ||
| ) | ||
| X_train_corr = X_train - design_train[:, 1:] @ beta_train[1:, :] | ||
| X_test_corr = X_test - design_test[:, 1:] @ beta_test[1:, :] | ||
|
|
||
| return X_train_corr, X_test_corr |
There was a problem hiding this comment.
Instead of this function I propose using this class to perform the fitting of the linear regression on the training set, and use these beta values for the test set as well:
class SiteRegressor(BaseEstimator, TransformerMixin):
"""
Regress out site effects from connectivity features.
"""
def __init__(self, n_connectivity_features: int):
self.n_connectivity_features = n_connectivity_features
def fit(self, X, y=None):
X_feat = X[:, :self.n_connectivity_features]
site = X[:, self.n_connectivity_features:]
# Add intercept
design = np.column_stack([np.ones(len(site)), site])
# Estimate coefficients on the training fold only
self.beta_ = np.linalg.lstsq(design, X_feat, rcond=None)[0]
return self
def transform(self, X):
X_feat = X[:, :self.n_connectivity_features]
site = X[:, self.n_connectivity_features:]
# Add intercept
design = np.column_stack([np.ones(len(site)), site])
# Remove only the site contribution (keep the intercept)
return X_feat - design[:, 1:] @ self.beta_[1:]
| with parallel_backend("threading", n_jobs=n_jobs): | ||
| cv_results = cross_validate( | ||
| pipe, | ||
| connectivity_data, | ||
| y_train, | ||
| cv=cv_strategy, | ||
| scoring=scoring_metrics, | ||
| n_jobs=n_jobs, | ||
| ) | ||
|
|
||
| scores_df = pd.DataFrame({k.replace("test_", ""): v for k, v in cv_results.items() if k.startswith("test_")}) |
There was a problem hiding this comment.
With the class function proposed above, we could go back to the previous workflow, and adding site correction directly within the scikit-learn pipeline:
pipe = Pipeline(
[
("imputer", SimpleImputer(strategy="median")),
("site_regression", SiteRegressor(connectivity_data.shape[1])),
("scaler", StandardScaler()),
("pca", PCA(n_components=n_pca, svd_solver="randomized", random_state=random_state,)),
("estimator", estimator),
]
)
with parallel_backend("threading", n_jobs=n_jobs):
cv_results = cross_validate(
pipe,
X,
y_train,
cv=cv_strategy,
scoring=scoring_metrics,
n_jobs=n_jobs,
)
|
@htwangtw I implemented the code with @pbergeret12 suggestion for now. But I agress, we should discuss this in more details |
htwangtw
left a comment
htwangtw
left a comment
There was a problem hiding this comment.
We are missing tests for the new features. I would prefer to have it added at some point, either through a test dataset with the site as an argument or a unit test. @mathdugre @pbergeret12 are you guys okay with this direction? We can open an issue to log it.
|
Hi everyone, |
There was a problem hiding this comment.
Pull request overview
Adds optional site-effect correction across the QC-FC and prediction portions of the wonkyconn workflow (Issue #132), plumbing a new CLI/config flag through the pipeline and extending tests to cover both corrected/uncorrected runs.
Changes:
- Add
--site-correctionflag and propagate it throughworkflow()/config, enforcing presence of asitephenotype column when enabled. - Include site as a covariate in QC-FC and add site-effect regression into the age/sex prediction cross-validation pipeline.
- Expand/adjust test coverage (including new prediction tests) and update typing/tooling/CI coverage configuration.
Reviewed changes
Copilot reviewed 22 out of 22 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
| wonkyconn/workflow.py | Wires site-correction flag into QC-FC + prediction; adds phenotype validation for site. |
| wonkyconn/run.py | Adds --site-correction CLI flag. |
| wonkyconn/config.py | Persists site_correction in config and namespace conversion. |
| wonkyconn/features/quality_control_connectivity.py | Adds optional site covariate to QC-FC partial-correlation model. |
| wonkyconn/features/age_sex_prediction.py | Introduces site-effect regression in the CV pipeline and plumbs sites through prediction scoring. |
| wonkyconn/features/tests/test_prediction.py | New tests validating that site-confounded signals drop to chance when correction is enabled. |
| wonkyconn/tests/test_cli.py | Runs smoketests with/without --site-correction. |
| wonkyconn/textual_app.py | Adjusts Textual imports/type-checker directives. |
| wonkyconn/tests/test_gradients_correlation.py | Minor typing-related adjustments to nibabel loading. |
| wonkyconn/tests/test_correlation.py | Adds pyright suppression for a call-site typing mismatch. |
| wonkyconn/tests/test_atlas.py | Updates atlas tests (signatures and typing suppression). |
| wonkyconn/features/tests/test_network.py | Adds pyright suppressions for datalad calls in tests. |
| wonkyconn/features/network.py | Renames accumulators + suppresses runtime warnings around NaN statistics/corrcoef. |
| wonkyconn/features/distance_dependence.py | Normalizes spearmanr return handling and scalar conversion. |
| wonkyconn/features/calculate_gradients_correlation.py | Tightens typing (Nifti1Image), adjusts nibabel load usage, and return typing. |
| wonkyconn/features/calculate_degrees_of_freedom.py | Minor typing annotation adjustment. |
| wonkyconn/atlas.py | Removes import-not-found ignores; tweaks masker usage; adds pyright suppressions. |
| pyproject.toml | Expands mypy ignore-missing-imports module list (brainspace/joblib/nilearn/sklearn). |
| codecov.yml | Adjusts Codecov comment settings and ignore globs. |
| .pre-commit-config.yaml | Adds additional mypy env deps (e.g., pytest/textual). |
| .github/workflows/test.yml | Appends coverage in smoketest light job. |
Suppressed comments (3)
wonkyconn/workflow.py:307
- This check also defaults
site_correctiontoTrue, which can raise on phenotype files that legitimately don't have asitecolumn when the flag wasn't provided. The fallback should match the CLI/config default (False).
if getattr(args, "site_correction", True):
wonkyconn/features/age_sex_prediction.py:48
- The
fit()docstring refers to aconnectivity_data_siteargument and suggests site dummy variables are already part ofX, butSiteRegressorderives dummies from the separatesitesarray. This mismatch makes the API hard to understand.
Args:
connectivity_data_site: A 2D array where the first n_connectivity_features columns
are connectivity features and the remaining columns are site dummy variables.
y: Ignored. This parameter exists for compatibility with the scikit-learn API.
wonkyconn/features/age_sex_prediction.py:69
- Similarly,
transform()'s docstring describesXas including site dummy variables, but the transformer computes dummies internally fromsitesandX.index. Align the docstring with the actual behavior.
Transform the connectivity data by regressing out site effects.
Args:
connectivity_data_site: A 2D array where the first n_connectivity_features columns
are connectivity features and the remaining columns are site dummy variables.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Also simplify SiteRegressor by using pandas DataFrames and add tests to ensure that site effects are removed and do not confound prediction after regression
- The stochastic gradient descent (saga) solver for logistic regression fails to converge on small datasets, which is a core use case for wonkyconn. We switch to lbfgs for greater performance on small data and better convergence - Replace fit_transform call
Also add type annotations to tests
Argparse already makes sure that every argument exists as an attribute, set to either the user-specified value, the default value or None if no default was given. This means that the default values defined in the getattr call will be ignored in most cases, so keeping them may lead to misunderstandings down the line.
325d709 to
a2b562b
Compare
|
Hi @HippocampusGirl ! Looks good to me! I am currently testing the new implementation on additional datasets, and so far it seems to be working correctly. Many thanks! |
…ts in their subgroups from the training pipeline
328982f to
e221347
Compare
Smoke test fails for ds000030
Previously, each group would be assigned the same file name, leading to each subsequent group overwriting the previous group's results. Now, each group will have its own file name.
Related to #132
Few changes made so far:
--site_correction, default: False)."age + C(gender) + C(site)". We could also use dummy variables if that seems like a cleaner approach—what do you think?We might need to update our test dataset to have the scanner column named "site" as required.
Todo:
--site_correction