This example will take a trivial case of having two regressors that are highly correlated and see if regressing the same component twice results in any re-introduction of artefacts.
The simplifying assumption is the noise components are identical between what is pulled out from ICA-AROMA as from another method (e.g. DVARS), while in reality they would be correlated, but not the same.
H0: If I apply the same regressor twice sequentially the "cleaned" data will be exactly the same.
H1: If I apply the same regressor twice sequentially the "cleaned" data will be different from the first to second application
This is taken from Chris Gorgolewski's investigation of ICA-AROMA
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
num_samples = 400
agg_signal_shared_variance = []
nonagg_signal_shared_variance = []
ortho_signal_shared_variance_deltas = []
orig_signal_shared_variance = []
agg_noise_shared_variance = []
nonagg_noise_shared_variance = []
ortho_noise_shared_variance = []
orig_noise_shared_variance = []
agg_shared_variance = []
nonagg_shared_variance = []
ortho_shared_variance = []
for _ in range(200):
# Generate the random samples.
unique_noise = np.random.randn(num_samples,1)
unique_signal = np.random.randn(num_samples,1)
shared_variance = np.random.randn(num_samples,1)
unexplained_bold_data = np.random.randn(num_samples,1)
noise_component = 0.2*unique_noise + 0.2*shared_variance
signal_component = 0.2*unique_signal + 0.2*shared_variance
bold_data = 0.4*unexplained_bold_data + 0.2*unique_noise + 0.2*shared_variance + 0.2*unique_signal
# JK: doing the first pass of regression
agg_res = np.dot(noise_component, np.dot(bold_data.T, np.linalg.pinv(noise_component).T))
newData_agg = bold_data - agg_res
# JK: regressing the denoised data with the same noise regressor
# this is supposed to simulate a situation where FD or DVARs shares a lot of variance with a noise component
# from ICA AROMA, in this most extreme version, I'm just using the same noise component as the regressor again.
agg_res2 = np.dot(noise_component, np.dot(newData_agg.T, np.linalg.pinv(noise_component).T))
newData_agg2 = newData_agg - agg_res2
# what is the correlation between the bold data, the noise, and the signal?
np.corrcoef(np.concatenate([bold_data, noise_component, signal_component], axis=1).T)
array([[ 1. , 0.50217411, 0.48331507], [ 0.50217411, 1. , 0.49249272], [ 0.48331507, 0.49249272, 1. ]])
# Plot the data on top of each other, supports H0
plt.plot(newData_agg, np.arange(0, 400))
plt.plot(newData_agg2, np.arange(0, 400))
[<matplotlib.lines.Line2D at 0x7f117d2dab38>]
# observe the difference in between the datasets (if it exists), supports H0
newData_agg - newData_agg2
array([[ -2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 8.32667268e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 8.32667268e-17], [ -5.55111512e-17], [ -1.34441069e-17], [ 0.00000000e+00], [ 8.32667268e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 1.38777878e-17], [ 5.55111512e-17], [ -5.55111512e-17], [ -2.77555756e-17], [ 2.77555756e-17], [ 0.00000000e+00], [ -2.77555756e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 2.77555756e-17], [ 5.55111512e-17], [ 0.00000000e+00], [ 2.94902991e-17], [ 0.00000000e+00], [ -1.38777878e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 2.77555756e-17], [ -1.11022302e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 2.77555756e-17], [ 0.00000000e+00], [ -5.55111512e-17], [ -3.46944695e-18], [ -1.11022302e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 1.38777878e-17], [ 5.55111512e-17], [ 0.00000000e+00], [ -5.55111512e-17], [ -2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ -1.38777878e-17], [ 2.77555756e-17], [ -1.11022302e-16], [ 2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 3.98986399e-17], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 4.33680869e-17], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ -8.32667268e-17], [ 0.00000000e+00], [ -2.77555756e-17], [ 2.77555756e-17], [ 0.00000000e+00], [ 8.32667268e-17], [ 5.55111512e-17], [ -2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -2.77555756e-17], [ 5.55111512e-17], [ -2.77555756e-17], [ -2.77555756e-17], [ -4.16333634e-17], [ 9.02056208e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 1.11022302e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -1.11022302e-16], [ 0.00000000e+00], [ 2.77555756e-17], [ 0.00000000e+00], [ 3.46944695e-17], [ 0.00000000e+00], [ 1.38777878e-17], [ 0.00000000e+00], [ -2.77555756e-17], [ 7.63278329e-17], [ -8.32667268e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 1.11022302e-16], [ 0.00000000e+00], [ -5.55111512e-17], [ -1.04083409e-16], [ -1.35308431e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 8.32667268e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ -1.11022302e-16], [ 0.00000000e+00], [ 5.55111512e-17], [ 2.77555756e-17], [ -4.16333634e-17], [ 0.00000000e+00], [ -3.12250226e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ -1.11022302e-16], [ -1.11022302e-16], [ 5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 1.11022302e-16], [ -2.08166817e-17], [ 4.16333634e-17], [ 0.00000000e+00], [ -5.55111512e-17], [ -4.16333634e-17], [ 5.55111512e-17], [ 5.55111512e-17], [ 2.77555756e-17], [ -2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ -8.32667268e-17], [ 5.55111512e-17], [ 0.00000000e+00], [ 9.71445147e-17], [ 1.73472348e-17], [ -1.11022302e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 8.32667268e-17], [ -5.55111512e-17], [ 5.55111512e-17], [ -2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 4.51028104e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ -5.55111512e-17], [ 1.11022302e-16], [ -2.77555756e-17], [ 9.02056208e-17], [ 0.00000000e+00], [ -1.38777878e-16], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 2.77555756e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ -4.16333634e-17], [ -5.55111512e-17], [ -1.38777878e-16], [ 0.00000000e+00], [ -5.55111512e-17], [ 5.46437895e-17], [ 2.77555756e-17], [ 5.55111512e-17], [ 8.32667268e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 2.77555756e-17], [ 0.00000000e+00], [ 1.11022302e-16], [ 0.00000000e+00], [ -2.77555756e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 1.38777878e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 8.32667268e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ -2.77555756e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 2.77555756e-17], [ -5.55111512e-17], [ 5.55111512e-17], [ -2.77555756e-17], [ 2.77555756e-17], [ 0.00000000e+00], [ 1.11022302e-16], [ 0.00000000e+00], [ -5.55111512e-17], [ -5.55111512e-17], [ 5.55111512e-17], [ 5.55111512e-17], [ -5.55111512e-17], [ 5.55111512e-17], [ -1.11022302e-16], [ -5.55111512e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ -5.63785130e-17], [ 5.55111512e-17], [ -3.03576608e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ -1.11022302e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ 1.11022302e-16], [ -5.55111512e-17], [ 9.02056208e-17], [ -1.38777878e-16], [ 2.77555756e-17], [ -5.55111512e-17], [ -5.55111512e-17], [ 3.46944695e-18], [ -5.55111512e-17], [ 0.00000000e+00], [ 2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ -1.38777878e-17], [ 5.55111512e-17], [ -1.38777878e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -2.77555756e-17], [ 5.55111512e-17], [ 2.77555756e-17], [ 1.11022302e-16], [ -3.81639165e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ -2.77555756e-17], [ 0.00000000e+00], [ 2.60208521e-18], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ -2.08166817e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -8.32667268e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 7.63278329e-17], [ 0.00000000e+00], [ 2.77555756e-17], [ -5.55111512e-17], [ 5.55111512e-17], [ 0.00000000e+00], [ 5.55111512e-17], [ 4.16333634e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 1.38777878e-17], [ 5.55111512e-17], [ 1.11022302e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 1.11022302e-16], [ 5.55111512e-17], [ 0.00000000e+00], [ -4.46691295e-17], [ 0.00000000e+00], [ 3.12250226e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -1.38777878e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ -8.32667268e-17], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 4.16333634e-17], [ -1.11022302e-16], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -2.77555756e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -1.11022302e-16], [ 5.55111512e-17], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 5.55111512e-17], [ 0.00000000e+00], [ -5.55111512e-17], [ 0.00000000e+00], [ 0.00000000e+00], [ 0.00000000e+00], [ -8.32667268e-17], [ 0.00000000e+00], [ 1.11022302e-16], [ -2.08166817e-17], [ 0.00000000e+00], [ 2.77555756e-17], [ -1.11022302e-16], [ 0.00000000e+00], [ -1.11022302e-16], [ 8.32667268e-17], [ 0.00000000e+00], [ 2.77555756e-17], [ 2.77555756e-17], [ -2.77555756e-17]])
# see if the correlations between components are
display(np.corrcoef(np.concatenate([newData_agg, noise_component, signal_component], axis=1).T))
display(np.corrcoef(np.concatenate([newData_agg2, noise_component, signal_component], axis=1).T))
array([[ 1.00000000e+00, 2.92960157e-04, 3.33706379e-01], [ 2.92960157e-04, 1.00000000e+00, 5.05627895e-01], [ 3.33706379e-01, 5.05627895e-01, 1.00000000e+00]])
array([[ 1.00000000e+00, 2.92960157e-04, 3.33706379e-01], [ 2.92960157e-04, 1.00000000e+00, 5.05627895e-01], [ 3.33706379e-01, 5.05627895e-01, 1.00000000e+00]])