If you’re using aggressive denoising, it’s certainly valid to include the mean WM, CSF, and global time series (as well as discrete cosine frequencies) in the nuisance model alongside the components identified as noise by ICA-AROMA and perform the regression in a single step. In this case, the output time series is orthogonalised with respect to all regressors in the model.
The situation becomes more complicated for the case of nonaggressive denoising. In particular (as discussed in the issue that you linked):
- any further regression steps will reduce the overall variance in the data, as measured by the residual sum of squares (less relevant here), and
- it’s possible that the orthogonality constraint on the residuals will result in the appearance of structured noise in the data that is related to previously removed regressors.
Intuitively, we might expect that – in order to prevent reintroduction of such structured noise – we would want to re-extract the mean CSF, WM, and global time series after performing the nonaggressive denoising, since they will already be orthogonal to any sources of variance removed during the nonaggressive denoising step.
Interestingly, however, previous simulations indicate that using the signal extracted prior to ICA-AROMA denoising confers a slight benefit in recovering a known ground-truth signal. I can’t justify this result theoretically, and it defies intuition to some extent. I’m sorry that I don’t have a more conclusive answer for you, but this remains an active area of research at this time.