Generate metrics from external regressors using F stats

[handwerkerd]

Closes #1009. This is an alternative approach to #1008 and #1021.

If a user provides external regressors, this will calculate a fit to those regressors to include as metrics in the component table and for use in decision trees.

Changes proposed in this pull request:

• Create a new module, metrics.external, for calculating metrics from external regressors
  • An F statistic is used with a Legendre polynomial detrending baseline and p, F, and R2 values are saved as metrics in the component table
  • These stats are saved for the Full F model using all regressors as well as partial F models. This makes it possible to both examine when the full model fits, but to also log if it's the motion vs ROI-based vs cardiac/respiration-based regressors that fit to a specific component time series.
  • There is a special class of regressors called task_keep If regressors under this label are included, these will be excluded from the full F model and a separate full F model will be calculated using just these regressors. This was a suggestion from @dowdlelt so that it would be possible to identify and conservatively retain components that fit to the overall task design.
• Add a new parameter to the tedana CLI, --external, to pass in a TSV with external regressors.
• Demo decision trees that incorporate external regressor correlation metrics. There are called demo_ because they demonstrate how to use the new functionality. They might be similar to trees we'll want to recommend, but validation on actual data (and tweaking) will be easier when this is merged and more people are using it.
  • In the decision tree there is a new field external_regressor_config which is a dictionary with the following keys:
    • info A description of what the config does that's saved
    • report A description to add to report.txt
    • calc_stats Currently the only option is "F" but this leaves open possibilities for additional options.
    • detrend: Will automatically calculate the number of polynomial detrending regressors if true, but can also be a number for the users to specify. If this is false, then will be set to 0 (mean removal) and log a warning.
    • f_states_partial_models [optional] A list of titles for the partial models (i.e. ["Motion", "CSF"])
    • If f_states_partial_models has model names, each model name is its own key and contains either a list of column labels or a regular express wildcard (i.e. "Motion": ["^mot_.*$"] means the Motion partial model will include any external regressor column label that begins with mot_ and "Motion": ["mot_x", "mot_y", "mot_z"] specifies 3 specific column label names
    • task_keep [optional] Contents are a regex wildcard or specific names to define task regressors that will not be included in the full model
  • demo_minimal_external_regressors_motion_task_models.json uses all the above options, uses the partial models to add a classification_tag, but not change results and retains components that correlate to the task (R2>0.5), have kappa>elbow irregardless of what rho is.
  • demo_minimal_external_regressors_single_model.json uses the minimum number of options to run with external regressors.
• Includes tests that should cover all new code.
  • In adding testing, I started to use the matching parameter to confirm the correct errors are triggered and caplog.text to confirm expected warnings or info statements are logged.
  • Also starting to add tests that check for specific numerical outputs.

To do:

• [ ] Update documentation to explain all this new functionality and outputs (waiting on this for enough of us to agree that the approach is generally good and we're not making more substantive changes to the structure/parameters/etc

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 88.05970% with 32 lines in your changes are missing coverage. Please review.

Project coverage is 89.61%. Comparing base (12aea7f) to head (4223bf9). Report is 1 commits behind head on main.

Files	Patch %	Lines
tedana/metrics/external.py	91.36%	7 Missing and 5 partials :warning:
tedana/selection/component_selector.py	83.01%	4 Missing and 5 partials :warning:
tedana/selection/selection_utils.py	85.71%	2 Missing and 1 partial :warning:
tedana/workflows/tedana.py	70.00%	2 Missing and 1 partial :warning:
tedana/metrics/collect.py	86.66%	1 Missing and 1 partial :warning:
tedana/stats.py	80.00%	1 Missing and 1 partial :warning:
tedana/metrics/_utils.py	94.44%	0 Missing and 1 partial :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1064      +/-   ##
==========================================
- Coverage   89.91%   89.61%   -0.31%     
==========================================
  Files          26       27       +1     
  Lines        3621     3853     +232     
  Branches      629      686      +57     
==========================================
+ Hits         3256     3453     +197     
- Misses        214      233      +19     
- Partials      151      167      +16     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[handwerkerd]

I'd appreciate some reviewers for this: @tsalo @eurunuela @goodalse2019 @n-reddy @dowdlelt . From my end, the only thing left to do is update the ReadTheDocs part of the documentation. I'm holding off on that until a couple of you agree that the overall structure makes sense.

[handwerkerd]

@tsalo I just added type hints to external.py I'm still a bit of a novice with type hints and it looks like formatting is non-trivially changing across python versions. I tried to follow guidance that's valid for python 3.8 since that's our oldest permitted version. I'll keep working in other files and on your remaining review comments, but let me know if I'm appropriately adding type hints so far.

~~Update as I was writing this comment: Lots of tests just failed so the answer to that question is"no". It looks numpy v1.24.3 that's currently in my python 3.8 environment includes numpy.typing and not numpy._typing while numpy v1.26.3 in my python 3.11 environment includes numpy._typing and not numpy.typing. I suspect there's a way to set up backwards compatibility for this, but I wanted to share what I'm seeing before I switch to childcare stuff.~~ I figured out the issue and all tests are passing now. It seems like using import numpy.typing as npt works, but directly referencing np.typing does not work across all versions

[handwerkerd]

@tsalo Can you explain more what you mean by "It would be great if strings and lines in docstrings were broken on punctuation. This will result in cleaner diffs in the future."

[tsalo]

Absolutely. Take a paragraph in the text, like the one below:

This is sentence one. This is
sentence two. This is sentence
three- but it's still going.

Since sentences are a unit within the paragraph, we're more likely to change them than random words. Plus, since almost everything we write is in markdown or rst, the newlines don't impact the rendered text- just the code.

-This is sentence one. This is
+This is sentence one. This is a 
-sentence two. This is sentence
+new sentence two with other
-three- but it's still going.
+changes. This is sentence three-
+but it's still going.

The diff on that, where the only consideration is line length, is much more extensive (and harder to interpret for a reviewer) than the diff if we broke up the text on punctuation:

 This is sentence one.
-This is sentence two.
+This is a new sentence 
+two with other changes.
 This is sentence three-
 but it's still going.

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 92.30769% with 21 lines in your changes missing coverage. Please review.

Project coverage is 89.90%. Comparing base (7ca9c27) to head (8d09599).

:exclamation: Current head 8d09599 differs from pull request most recent head 2d45383

Please upload reports for the commit 2d45383 to get more accurate results.

Files	Patch %	Lines
tedana/selection/component_selector.py	83.33%	4 Missing and 5 partials :warning:
tedana/selection/selection_utils.py	85.71%	2 Missing and 1 partial :warning:
tedana/workflows/tedana.py	70.00%	2 Missing and 1 partial :warning:
tedana/metrics/collect.py	86.66%	1 Missing and 1 partial :warning:
tedana/stats.py	80.00%	1 Missing and 1 partial :warning:
tedana/metrics/_utils.py	94.44%	0 Missing and 1 partial :warning:
tedana/metrics/external.py	99.24%	0 Missing and 1 partial :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1064      +/-   ##
==========================================
- Coverage   89.91%   89.90%   -0.02%     
==========================================
  Files          26       27       +1     
  Lines        3620     3854     +234     
  Branches      629      685      +56     
==========================================
+ Hits         3255     3465     +210     
- Misses        214      226      +12     
- Partials      151      163      +12     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[handwerkerd]

I think I've now addressed all feedback regarding code and I've started to add documentation to building_decision_trees.rst I still need to make flow charts for the new trees and add them to included_decision_trees.rst but otherwise this PR is (again) ready for review & for people to try it out.

[handwerkerd]

I think/hope this is fully done and ready for review. It's important to have people look at the code, but I'd also like to have people read the documentation (tedana --help https://tedana--1064.org.readthedocs.build/en/1064/building_decision_trees.html#external-regressor-configuration & https://tedana--1064.org.readthedocs.build/en/1064/included_decision_trees.html#demo-external-regressors-single-model) and try to run this to see if it both makes sense and runs as expected. (@tsalo @eurunuela @n-reddy @goodalse2019 @dowdlelt )

[handwerkerd]

I really don't want to hold this PR up any more. I say we merge as-is, make a release-candidate release, and deal with any bugs or typos before the next actual release.

I'd still like to have @n-reddy and/or @goodalse2019 at least look over the documentation and run before merging. It would be good to see this understandable and executable by someone who hasn't yet look deep into the code. If they have time in the next week or so, I'm happy to wait before merging.

[tsalo]

In that case, can we commit to merging on a specific date (e.g., next Friday- 08/09) if they don't provide feedback by then?

[goodalse2019]

@handwerkerd and @tsalo I will have finished looked over the documentation and run by tomorrow evening. I have tested it on one subject very basic to make sure it ran and all looked like it went well! I am going to try on a couple others expanding my regressors tsv file a bit more. For testing purposes - is there a lower end of # of subjects you prefer this to be run over (e.g., at least 10?) or will just a couple do? Also - is there a good way to provide a overview of results for people to look over just to see the tests? Like uploading pdf or ppt in mattermost with some of the output figures?

[handwerkerd]

I'm at a conference 8/6-9 so 8/12 would be a good day to merge if we don't get feedback requiring changes. (Earlier if we get feedback saying everything is fine.

@goodalse2019 At this point, we're not claiming the results will be better or worse, so I'd just like to see this functioning as expected on typical data. If you have more than one dataset sitting around, then trying this on 2-3 subjects with different acquisition parameters and making sure it runs as expected would be fine. Seeing the results will be useful for the next step of trying to figure out what's actually beneficial. If it's not too large, I think you can upload a pdf as a comment on this PR, which will mean we'd be able to find it when continuing work on these methods. Thank you.

[goodalse2019]

Sounds good, we have two - I can test on a couple subject from both and will upload some results asap

[goodalse2019]

I've run the single model demo tree and attached are some output summaries after running on n = 2 in two me datasets:

Dataset 1A - single model Dataset1_singlemodel.pdf Dataset 2 - single model Dataset2_singlemodel.pdf

Still working on the task_motion_models demo - dataset 1 has a task and I am trying to get that to work but first few attempts were no good on my end. I can run over a few more subs in each if you think helpful and update these documents. I will be out of town in a wedding 8/2-8/4 so won't be until 8/5 that I can pick up again.

[tsalo]

Thanks @goodalse2019. I think the results look reasonable, and the fact that it ran through is good news. @handwerkerd are you comfortable merging given her results?

[n-reddy]

I reviewed the documentation and it makes sense to me! I just made a few small edits for clarity. Still planning to test the code if I get a chance in the next couple days.

[handwerkerd]

Thank you all. I'm comfortable merging, but I'm currently traveling through Tuesday and wouldn't mind giving it one more quick check myself. If others are comfortable merging, then give approving reviews.

[n-reddy]

Just an update that both the single model and the task_motion_models worked well on two datasets I tested with different parameters! A small comment perhaps for future tweaks of the demo .json is that the classification tags on the demo_external_regressors_motion_task_models can be a little confusing. For example, see in the pdfs linked below that the highlighted component was tagged as 'Likely BOLD, External Regressors, Unlikely BOLD' in the motion and task model.

Minimal model: minimal.pdf

External regressors - Motion parameters only (single model): external_regressors_motion.pdf

External regressors - Motion parameters and task regressors: external_regressors_motionTask.pdf