Repository: SIMEXP/load_confounds (language: python)

Issue:
`nilearn` behaviour with PSC changed
Found when handling PR #96 
Nilearn version 0.7.0
This issue doesn't concert the main package, but should be noted for our related test `test_nilearn_standardize_psc`. The PSC strategy caused confounds not being regressed out in the data when confound variable has mean of 0 or less. Warning message:
```
nilearn/signal.py:81: UserWarning: psc standardization strategy is meaningless for features that have a mean of 0 or less. These time series are set to 0.
```
This issue will be reported to nilearn, and we should expect future changes to this test. 
#52 ADD three ICA-AROMA related strategy
Three new ICA-AROMA related strategy added to address issue #52:
1. (non-aggressive) ICA-AROMA from Ciric 2017
2. (non-aggressive) AROMA-GSR from Ciric 2017
3. aggressive ICA-AROMA from Pruim et al. 2015

We can decide if we want strategy 3 or not since it's not really in our main reference.

Strategy 3 was the most straight forward one to implement. There are a lot of debates on Strategy 1 and Strategy 2. I followed the discussion on [this neurostars thread](https://neurostars.org/t/ica-aroma-agg-vs-non-agg/3708) to implement the first two. According to the thread, the current implementation is supported by the simulation results of fmriprep notebooks `05`: https://github.com/nipreps/fmriprep-notebooks/
I need more pairs of eyes to double-check this part!

I added new test files for ICA-AROMA just in case I will break the existing tests.
The test confound file was chopped from the one of `sub-01` in [the test dataset of `fitlins`](https://datasets.datalad.org/?dir=/labs/poldrack/ds003_fmriprep/)
BUT I broke the test anyway :see_no_evil: 
I kept all the old tests with the original file, but used the new one for my strategies. I found out that having two sets of confound files of different length broke `_simu_img` and some related utilities/tests. The fix was to add confound file as an input variable to the simulation was generated based on related files.

On a similar note, the tests on existing strategies can also be improved a bit. I will describe it in a different issue tomorrow.

Finally, sorry for the messy commit history!! Accidentally committed some large file so I did some clumsy git-gymnastics :sweat_smile: 


Discussion / hints:
thanks for pointing this out! 

@all-contributors add @htwangtw for bug
@pbellec 

I've put up [a pull request](https://github.com/SIMEXP/load_confounds/pull/107) to add @htwangtw! :tada:
My PR that allow negative mean signal has now been merged to nilearn master branch. It will be in the next bug fix release of 0.7. Will continue to watch this spot and fix up the test!
awesome! thanks for the update!
Also I suggest to ignore the codeclimate warnings. I don't think length is an issue for strategies.py, because of the ultra-modular nature of the code in there. We are working on `_load_single` as part of #94 
## Some notes related to ICA AROMA implementation in fMRIPrep.
- ICA-AROMA denoising is performed in MNI152NLin6Asym space only.
- Extra files from fMRIPrep when ICA aroma is on:
```*_AROMAnoiseICs.csv```: ICA component number that are classified as noise; should match relevant headers in `*desc-confounds_regressors.tsv`
```*_desc-MELODIC_mixing.tsv```: MELODIC mixing matrix, all ICA components; noise components specified in ```*_AROMAnoiseICs.csv``` are saved in `*desc-confounds_regressors.tsv` with headers.
```*_desc-smoothAROMAnonaggr_bold.nii.gz```: non-aggressive ICA-AROMA output produced by fMRIprep
- one can manually recreate `*space-MNI152NLin6Asym_desc-smoothAROMAnonaggr_bold.nii.gz` with the following files and `fsl_regfilt`:
```
fsl_regfilt -i sub-<subject_label>_task-<task_id>_**space-T1w_desc-preproc_bold**.nii.gz \
    -f $(cat sub-<subject_label>_task-<task_id>_**AROMAnoiseICs**.csv) \
    -d sub-<subject_label>_task-<task_id>_**desc-MELODIC_mixing**.tsv \
    -o sub-<subject_label>_task-<task_id>_space-T1w_desc-AROMAnonaggr_bold.nii.gz
```

## In Ciric 2017, the ICA-AROMA related models are performed as follow:
>ICA-AROMA includes two de-noising steps. The first de-noising step occurs immediately after classification. All component time series (signal and noise) are included as predictors in the linear model, and the residual BOLD time series is obtained via partial regression of only the noise time series. 

I believe this generates ```*_desc-smoothAROMAnonaggr_bold.nii.gz``` 

> A second confound regression step occurs after temporal filtering, wherein mean signals from WM and CSF (models 13 and 14) and the global signal (model 14 only) were regressed from the data.

This is what I have for `load_confound` to grab using `ICAAROMA` (model 13) and `AROMAGSR` (model 14).

## "Aggressive ICA-AROMA" from Pruim 2015
https://doi.org/10.1016/j.neuroimage.2015.02.064
> ... ICs can be fully regressed out of the data, analogous to RP-regression (Friston et al., 1996, Satterthwaite et al., 2013a). In this case, all variance associated with these noise ICs will be removed, including shared-variance with ‘good’ ICs. Therefore, such a component regression method can be regarded as ‘aggressive denoising’. 

This super vague description is implemented as `AggrICAAROMA`, presumably performed on the regular `*_desc-preproc_bold.nii.gz` only.
This is looking great and is almost ready to be merged! 

Obviously we need to fix the tests. I have reviewed how these tests work and I think it's pretty logical (pfiou) 
The problem here is I think the very limited number of degrees of freedom (30 time points - number of confounds). We  could increase the tolerance threshold to something like 0.4. I like better the idea of artificially increasing the number of degrees of freedom. This could be achieved by stacking the confounds on itself in that line:
https://github.com/htwangtw/load_confounds/blob/ica_aroma/load_confounds/tests/test_parser.py#L25 

By adding a single line that duplicates the rows twice or thrice, this should work. The simulation data is built based on the confounds so there should be no need to change anything else in the code. 
@pbellec thank you for finding the solution! It did cause the assertion in `_regression` to fail, but it's an easy fix.
The concerning part is `test_nilearn_standardize_psc` continues to fail, and I noticed a warning from `nilearn`:
```
nilearn/signal.py:81: UserWarning: psc standardization strategy is meaningless for features that have a mean of 0 or less. These time series are set to 0.
    warnings.warn('psc standardization strategy is meaningless '
```
~I did check `demean` was off in that test so I am not sure what is going on...~
Edit:
Among the six basic motion parameters in the new test, there happens to five out of six with negative mean, hence the warning. Nilearn would then set the invalid values to 0 in that case:
https://github.com/nilearn/nilearn/blob/d2ad41a9fa0aa7041c7e826518c6c3692a92a2e1/nilearn/signal.py#L81

Not sure if this should be seen as an issue on nilearn side or an unfortunate test data?
So it looks like nilearn's behaviour with PSC has changed. I suggest to change the test so it expects to fail. Then the test will pass. Please open a separate issue on the PSC behaviour, and we will deal with it separately. 
# [Codecov](https://codecov.io/gh/SIMEXP/load_confounds/pull/96?src=pr&el=h1) Report
> Merging [#96](https://codecov.io/gh/SIMEXP/load_confounds/pull/96?src=pr&el=desc) (c2366f9) into [master](https://codecov.io/gh/SIMEXP/load_confounds/commit/a13646d5353909e9a039f0185433f8dce182bb03?el=desc) (a13646d) will **increase** coverage by `0.31%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/SIMEXP/load_confounds/pull/96/graphs/tree.svg?width=650&height=150&src=pr&token=7KD817T4LJ)](https://codecov.io/gh/SIMEXP/load_confounds/pull/96?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master      #96      +/-   ##
==========================================
+ Coverage   97.58%   97.90%   +0.31%     
==========================================
  Files           6        6              
  Lines         414      477      +63     
==========================================
+ Hits          404      467      +63     
  Misses         10       10              
```


| [Impacted Files](https://codecov.io/gh/SIMEXP/load_confounds/pull/96?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [load\_confounds/parser.py](https://codecov.io/gh/SIMEXP/load_confounds/pull/96/diff?src=pr&el=tree#diff-bG9hZF9jb25mb3VuZHMvcGFyc2VyLnB5) | `97.45% <100.00%> (+0.10%)` | :arrow_up: |
| [load\_confounds/strategies.py](https://codecov.io/gh/SIMEXP/load_confounds/pull/96/diff?src=pr&el=tree#diff-bG9hZF9jb25mb3VuZHMvc3RyYXRlZ2llcy5weQ==) | `100.00% <100.00%> (ø)` | |
| [load\_confounds/tests/test\_parser.py](https://codecov.io/gh/SIMEXP/load_confounds/pull/96/diff?src=pr&el=tree#diff-bG9hZF9jb25mb3VuZHMvdGVzdHMvdGVzdF9wYXJzZXIucHk=) | `100.00% <100.00%> (ø)` | |
| [load\_confounds/tests/test\_strategies.py](https://codecov.io/gh/SIMEXP/load_confounds/pull/96/diff?src=pr&el=tree#diff-bG9hZF9jb25mb3VuZHMvdGVzdHMvdGVzdF9zdHJhdGVnaWVzLnB5) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/SIMEXP/load_confounds/pull/96?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/SIMEXP/load_confounds/pull/96?src=pr&el=footer). Last update [a13646d...c2366f9](https://codecov.io/gh/SIMEXP/load_confounds/pull/96?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

I now have made the change to reflect the suggestions, with a small extra to check the User warning of psc come from the line we are expecting. Now all tests should pass. I will open an issue related to PSC change, and we can discuss how to move this forward. 
@all-contributors add @htwangtw for code, data, doc, test
@pbellec 

I've put up [a pull request](https://github.com/SIMEXP/load_confounds/pull/106) to add @htwangtw! :tada:
This is is great, huge thanks @htwangtw !! It's a really valuable feature for `load_confounds`. Merging.

Fix patch:
diff --git a/requirements.txt b/requirements.txt
index 750c66a..4aec963 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,4 +2,4 @@ numpy>=1.17.4
 pandas>=0.25.3
 scikit-learn>=0.21.3
 scipy>=1.3.2
-nilearn>=0.6.2
+nilearn>=0.7.1
diff --git a/setup.py b/setup.py
index 2581e7f..abea2c3 100644
--- a/setup.py
+++ b/setup.py
@@ -29,7 +29,7 @@
         "pandas>=0.25.3",
         "scikit-learn>=0.21.3",
         "scipy>=1.3.2",
-        "nilearn>=0.6.2",
+        "nilearn>=0.7.1",
     ],  # external packages as dependencies
     classifiers=[
         "Development Status :: 3 - Alpha",


Test patch:
diff --git a/load_confounds/tests/test_parser.py b/load_confounds/tests/test_parser.py
index f385b49..13e4963 100644
--- a/load_confounds/tests/test_parser.py
+++ b/load_confounds/tests/test_parser.py
@@ -170,14 +170,10 @@ def test_nilearn_standardize_psc():
 
     img, mask_conf, mask_rand, X = _simu_img(demean=False)
 
-    # Areas with
-    with pytest.warns(UserWarning) as records:
-        tseries_raw, tseries_clean = _denoise(img, mask_conf, X, "psc")
-
-    nilear_warning = sum('psc standardization strategy' in str(r.message) for r in records)
-    assert nilear_warning == 1
+    # Areas with confound
+    tseries_raw, tseries_clean = _denoise(img, mask_conf, X, "psc")
     corr = _corr_tseries(tseries_raw, tseries_clean)
-    assert not corr.mean() < 0.2  # expect failure
+    assert corr.mean() < 0.2
 
     # Areas with random noise
     tseries_raw, tseries_clean = _denoise(img, mask_rand, X, "psc")