Hi Alfonso,

This is my-own follow-up. Actually I created a matlab script that shows how the signal changes by additionally adding more regressors into the design matrix for denoising.

In the attached image, the results with the threshold of 0.55 (left) and 0.99 (right) were compared with distribution of correlation coefficients (upper) and standardized timeseries and correlation map of a single slice (lower). The timeseries were standardized using Z-score only to see changes from signed values (residuals) instead of percentage of change. I still agree with you on the risk of misleading in using Z-score for timeseires; an isolated Z-value (even if it's high or low) doesn't mean much here. Instead, the question is whether or not there is a globally (over all voxels) synchronized (and also synchronized with the head motion estimation) and abrupt perturbation, which is highly likely due to head motion. 

Although I bit surprising, the threshold for WM/CSF voxels did not make a huge difference by adding the WM/CSF regressors, but thus increased p-value (or goodness-of-fit) of K-S test, i.e., it made the distribution of correlation coefficients more Gaussian after the denoising regression with global signal (from 0.25 to 0.37).

I guess that's probably because low threshold definition of WM and CSF includes some effects from GM voxels, thus along with the GM mean signal the collinearity, CompCor regressors with low threshold would decreased. Although I don't know how exactly it affected.

But without inclusion of global signal, the difference was very small; next to negligible, I'd say.

Best,