[FIX] Include ignored components in ME-DN T1c time series

[tsalo]

The ME-DN time series should include the ignored ICA components, but does not. The current approach is to regress out all ICA components from the data to get resid, and then to combine that with the BOLD-T1c (GS-controlled accepted components) time series as the ME-DN T1c time series. With this PR, the residuals time series will not have the ignored components regressed out, and thus the ME-DN time series will include those components.

Changes proposed in this pull request:

• Add rej (index of rejected components) argument to gscontrol_mmix.
• Compose ICA residuals time series (resid) of the data minus all of the ICA components except the ignored components.
• There are some minor improvements to documentation as well.

[tsalo]

As expected, this breaks the integration test on the dn_ts_OC_T1c.nii file, but just on that file.

[emdupre]

This is a great catch, @tsalo ! I'd really love @handwerkerd's review and to hear his opinion on the T1c series, generally, if he has time.

I'm just not sure I see the benefit of regressing the global signal instead of, say, the rejected comments. So, I'm thinking it might make sense to repurpose this into a kind of aggressive denoising, as discussed in #39. WDYT ?

[tsalo]

Isn't the removal of T1 global signal largely there based on the Power paper? I know that this particular approach, while recommended, hasn't been evaluated in any publications, which might be one reason to avoid making it mandatory. When we add in GODEC in #120, we'll be adding an argument (currently ws_denoise) to allow users to specify what kind of post-tedana denoising, if any, they'd like to implement in order to remove "widespread noise." We can make None the default at that point, I think.

I know that it's probably a good idea to orthogonalize the bad components with respect to the good ones, and we can do that with lstsq in writeresults pretty easily, right? I don't know how that interacts with the T1c GSR, although I wasn't involved in #39 and couldn't figure out where GSR was discussed. Could you elaborate on that a bit more?

[tsalo]

Regarding the orthogonalization, would it work to do the following in the tedana workflow, after selcomps?

if orth:  # a new bool argument
    a = mmix[:, acc]
    midkrej = np.hstack((midk, rej))
    b = mmix[:, midkrej]
    betas = np.linalg.lstsq(a, b, rcond=None)[0]
    pred_b = np.dot(a, betas)
    resid = b - pred_b
    mmix[:, midkrej] = resid

And then mmix (orthogonalized) can be fed into writeresults and gscontrol_mmix.

[handwerkerd]

I've been a bit behind at looking into code, but please please please do not regress out the whole brain global signal across the brain as a default setting. I'm familiar with the Power papers, but there are serious issues with setting GSR as a default. I can talk more about this at some point, but not this week.

[tsalo]

Just to be clear, are you referring to the standard GSR that is done to the concatenated and optimally combined data here or to the minimal image regression (aka mmix GSR, aka T1-correction, aka who-knows-what because there's no canonical name afaik) done here (or both)?

Both are done by default, but the T1-correction outputs additional files (and so doesn't have to affect results, depending on which files people use), while the standard GSR impacts all downstream outputs.

I'm hoping to re-integrate GODEC in #120, which will (in my mind) involve modularizing the post-TEDICA processing and will offer an alternative to the T1-correction. We'll need to add an argument for the post-TEDICA processing to perform, but I do have the default as None, which I think is fair.

[handwerkerd]

I think both. I don't remember GSR as a standard step in earlier pipelines, but I'm not sure I specifically looked for it. If I'm following the first one correctly regressing out the global signal before ICA seems particularly problematic.

[tsalo]

As best I can tell, both GSR steps were first added here, on 2017/05/23, in the initial shift from v2.5 to v3.2, so I don't think either step is used in any "official" releases of the pipeline.

Given that we're planning to roll back to v2.5 in #119, perhaps we can use that PR to change the defaults to what they were at that time. We can still retain the T1-correction, but can turn it off by default (as is the plan in #120).

[codecov]

Codecov Report

Merging #125 into master will decrease coverage by 0.3%.
The diff coverage is 0%.

@@            Coverage Diff             @@
##           master     #125      +/-   ##
==========================================
- Coverage   52.96%   52.65%   -0.31%     
==========================================
  Files          32       32              
  Lines        1903     1899       -4     
==========================================
- Hits         1008     1000       -8     
- Misses        895      899       +4

Impacted Files	Coverage Δ
tedana/io.py	35.55% <0%> (-1.09%)	⬇️
tedana/utils.py	98.01% <0%> (-0.08%)	⬇️
tedana/decay.py	92% <0%> (ø)	⬆️
tedana/tests/test_utils.py	100% <0%> (ø)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 63d643a...5c79624. Read the comment docs.

[tsalo]

Unit tests and the first integration test (running the pipeline) are passing, so I think this is ready to merge. I would like to keep this PR concise, so I'd prefer to leave changing the defaults for standard GSR and T1c-GSR until another PR.

[emdupre]

It looks like this is ready -- to be clear, this is just updating the ME-DN-T1c to (theoretically) match the creation of the ME-DN, no ?

I know we had a lot of changes since this was first created, so I wanted to make sure it's not intended to do anything else !

[tsalo]

That's right, that's all the changes are meant to do.

I just realized that I forgot to update the docstring for gscontrol_mmix, so there will be one more commit in just a second to fix that.

[emdupre]

OK, this looks good to me ! I'll merge this in 🎉