Nuisance Correction

Evaluation of White Matter PCs for Nuisance Regression

Goal: Determine if the nuisance correction method of removing the first five principal components of the white matter from the fMRI data is statistically supported.

To do this, an F-test was performed in order to determine a p-value that represents how well the white matter model fits the data. The HCP dataset was utilized for twenty subjects performing the emotion and motor tasks described here.

Why an F-test?

An F-test of overall significance in regression "compares the fits of different linear models." This test as we apply it compares the model from the white matter components to an "intercept-only model,... a regression model that contains no predictors." The null hypothesis is therefore that the white matter model and intercept-only model are equal with the alternative hypothesis being that the fit of the intercept-only model is "significantly reduced" compared to the white matter model. Thus if we obtain a p-value lower than the significance level we will specify as 0.05, then we reject the null hypothesis and conclude that the white matter model provides a better fit than the intercept-only model, which means the white matter model "predicts the response variable better than the mean of the response." More info

Description

We first calculate the beta values of the GLM of the white matter components. The GLM system of equations can be described as: $$y = Xb + e$$ where $X$ is the design matrix, $b$ is the beta value vector, and $e$ is the added error vector. The values predicted by the model are $\hat{y} = Xb$. We then calculate the sum of squares for the model as $SSM = \sum_{i=1}^n (\hat{y_i}-\bar{y})^2$ and the sum of squares for error as $SSE = \sum_{i=1}^n (y_i-\hat{y_i})^2$. We calculate an $F$ statistic with $DFM = p-1$ degrees of freedom in the numerator and $DFE = n-p$ degrees of freedom in the denominator where $n$ is the number of timesteps and $p$ is the number of regressors. The $F$ statistic is then calculated as follows: $$F_{DFM,DFE} = \frac{SSM \times DFE}{SSE \times DFM}$$ From here, we can use the value of the $F$ statistic and the degrees of freedom to calculate the p-value. A high F-value corresponds to a low p-value.

Implementation

After running the datasets through the pipeline in order to obtain the registered fMRI images, the following shell script was utilized in order to perform the F-tests on each of the subjects' task data.



In [ ]:

    
#!/bin/bash
# analysis.sh
#$1 = reg_fmri
#$2 = reg_struct

exec 4<$1
exec 5<$2

while read -r reg_fmri <&4 && read -r reg_struct<&5; do
    python analysis.py $reg_fmri $reg_struct
done 4<$1 5<$2

The script calls analysis.py which is as follows:



In [ ]:

    
from ndmg.utils import utils as mgu
import nibabel as nb
import numpy as np
import scipy.signal as signal
import scipy.stats as stats
import sys
fmri = sys.argv[1]
amri = sys.argv[2]
amask = '../../FNGS_server/atlases/mask/MNI152_T1-2mm_brain_mask.nii.gz'
gmask = '../avg152T1_gray_bin.nii'
an = 1
outdir = './'
lvmask = '../../FNGS_server/atlases/mask/HarvOx_lv_thr25-2mm.nii.gz'
amask_im = nb.load(amask)
amm = amask_im.get_data()
gmask_im = nb.load(gmask)
gmm = gmask_im.get_data()
lv_im = nb.load(lvmask)
lvm = lv_im.get_data()
fmri_name = mgu().get_filename(fmri)
anat_name = mgu().get_filename(amri)
nuisname = "".join([anat_name, "_nuis"])
map_path = mgu().name_tmps(outdir, nuisname, "_map")
cmd1 = " ".join(["fast -t", str(int(an)), "-n 3 -o", map_path, amri])
mgu().execute_cmd(cmd1)
wm_prob = map_path + "_pve_2.nii.gz"
prob = nb.load(wm_prob)
prob_dat = prob.get_data()
mask = (prob_dat > .99).astype(int)
for i in range(0, 2):
    erode_mask = np.zeros(mask.shape)
    x, y, z = np.where(mask != 0)
    if (x.shape == y.shape and y.shape == z.shape):
        for j in range(0, x.shape[0]):
            md = mask.shape
            if (mask[x[j],y[j],z[j]] and
                mask[np.min((x[j]+1, md[0]-1)),y[j],z[j]] and
                mask[x[j],np.min((y[j]+1, md[1]-1)),z[j]] and
                mask[x[j],y[j],np.min((z[j]+1, md[2]-1))] and
                mask[np.max((x[j]-1, 0)),y[j],z[j]] and
                mask[x[j],np.max((y[j]-1, 0)),z[j]] and
                mask[x[j],y[j],np.max((z[j]-1, 0))]):
                erode_mask[x[j],y[j],z[j]] = 1
    else:
        raise ValueError('Your mask erosion has an invalid shape.')
    mask = erode_mask

wmm = mask
mask = gmm
for i in range(0, 1):
    erode_mask = np.zeros(mask.shape)
    x, y, z = np.where(mask != 0)
    if (x.shape == y.shape and y.shape == z.shape):
        for j in range(0, x.shape[0]):
            md = mask.shape
            if (mask[x[j],y[j],z[j]] and
                mask[np.min((x[j]+1, md[0]-1)),y[j],z[j]] and
                mask[x[j],np.min((y[j]+1, md[1]-1)),z[j]] and
                mask[x[j],y[j],np.min((z[j]+1, md[2]-1))] and
                mask[np.max((x[j]-1, 0)),y[j],z[j]] and
                mask[x[j],np.max((y[j]-1, 0)),z[j]] and
                mask[x[j],y[j],np.max((z[j]-1, 0))]):
                erode_mask[x[j],y[j],z[j]] = 1
    else:
        raise ValueError('Your mask erosion has an invalid shape.')
    mask = erode_mask

gmm = mask
fmri_im = nb.load(fmri)
amri_im = nb.load(amri)
fmri_dat = fmri_im.get_data()
voxel = fmri_dat[fmri_dat.sum(axis=3) > 0, :].T
wm_ts = fmri_dat[wmm != 0, :].T
lv_ts = fmri_dat[lvm != 0, :].T
t = voxel.shape[0]
lin_reg = np.array(range(0, t))
quad_reg = np.array(range(0, t))**2
csf_reg = lv_ts.mean(axis=1)
wm_ts = wm_ts[:, wm_ts.std(axis=0) != 0]
wm_ts = signal.detrend(wm_ts, axis=0, type='linear')
wm_ts = wm_ts - wm_ts.mean(axis=0)
wm_ts = np.divide(wm_ts, wm_ts.std(axis=0))
U, s, V = np.linalg.svd(wm_ts)
wm_reg = U[:, 0:5]
R = np.column_stack((np.ones(t), lin_reg, quad_reg, wm_reg, csf_reg))
coefs = np.linalg.inv(R.T.dot(R)).dot(R.T).dot(voxel)
yhat = R.dot(coefs)
res = voxel - yhat
ssm = np.sum((yhat-voxel.mean(0))**2, axis=0)
sse = np.sum(res**2, axis=0)
dfm = coefs.shape[0]-1
dfe = voxel.shape[0]-coefs.shape[0]
msm = ssm/dfm
mse = sse/dfe
f = msm/mse
f_brain = np.zeros(fmri_dat.shape[:3])
f_brain[fmri_dat.sum(axis=3) > 0] = f.T
img = nb.Nifti1Image(f_brain, header=fmri_im.header, affine=fmri_im.affine)
nb.save(img, fmri[11:17] + '_fvals.nii.gz')
p = 1 - stats.f.cdf(f_brain[f_brain != 0], dfm, dfe)
p_white = 1 - stats.f.cdf(f_brain[(f_brain != 0) & (wmm != 0)], dfm, dfe)
p_grey = 1 - stats.f.cdf(f_brain[(f_brain != 0) & (gmm != 0)], dfm, dfe)
p_other = 1 - stats.f.cdf(f_brain[(f_brain != 0) & (wmm != 0) & (gmm !=0)], dfm, dfe)
fi = open('pvals_all.txt', 'a+')
for val in p.flat:
    fi.write(repr(val) + '\n')

fi.close()
fi = open('pwhite_all.txt', 'a+')
fm = open('pwhite_mean.txt', 'a+')
fp = open('pwhite_' + fmri[11:17] + '.txt', 'w')
for val in p_white.flat:
    fi.write(repr(val) + '\n')
    fp.write(repr(val) + '\n')

fm.write(repr(np.mean(p_white)) + '\n')
fi.close()
fm.close()
fp.close()
fi = open('pgrey_all.txt', 'a+')
fm = open('pgrey_mean.txt', 'a+')
fp = open('pgrey_' + fmri[11:17] + '.txt', 'w')
for val in p_grey.flat:
    fi.write(repr(val) + '\n')
    fp.write(repr(val) + '\n')

fm.write(repr(np.mean(p_grey)) + '\n')
fi.close()
fm.close()
fp.close()
fi = open('pother_all.txt', 'a+')
fm = open('pother_mean.txt', 'a+')
fp = open('pother_' + fmri[11:17] + '.txt', 'w')
for val in p_other.flat:
    fi.write(repr(val) + '\n')
    fp.write(repr(val) + '\n')

fm.write(repr(np.mean(p_other)) + '\n')
fi.close()
fm.close()
fp.close()

Simulation

We show a basic simulation below to prove the validity of our method. The first simulation shows a 2x2x2 "brain" with 100 timesteps where the voxel values for the fMRI and white matter are independently randomly chosen. In this simulation we expect little correlation between the timeseries and therefore the p-values should be rather high with little to none below the 0.05 threshold. The second simulation shows a 2x2x2 "brain" with 100 timesteps where the voxel values for the fMRI and white matter are set as the same randomly chosen values. In this simulation we expect total correlation between the timeseries and therefore the p-values should all be below our threshold. As we can see, our simulated results are as expected.



In [23]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal
import scipy.stats as stats
fmri_dat = np.random.rand(2,2,2,100)
wm = np.random.rand(2,2,2,100)
voxel = fmri_dat[fmri_dat.sum(axis=3) > 0, :].T
wm_ts = wm.T
t = voxel.shape[0]
lin_reg = np.array(range(0, t))
quad_reg = np.array(range(0, t))**2
wm_ts = wm_ts[:, wm_ts.std(axis=0) != 0]
wm_ts = signal.detrend(wm_ts, axis=0, type='linear')
wm_ts = wm_ts - wm_ts.mean(axis=0)
wm_ts = np.divide(wm_ts, wm_ts.std(axis=0))
U, s, V = np.linalg.svd(wm_ts)
wm_reg = U[:, 0:5]
R = np.column_stack((np.ones(t), lin_reg, quad_reg, wm_reg))
coefs = np.linalg.inv(R.T.dot(R)).dot(R.T).dot(voxel)
yhat = R.dot(coefs)
res = voxel - yhat
ssm = np.sum((yhat-voxel.mean(0))**2, axis=0)
sse = np.sum(res**2, axis=0)
dfm = coefs.shape[0]-1
dfe = voxel.shape[0]-coefs.shape[0]
msm = ssm/dfm
mse = sse/dfe
f = msm/mse
f_brain = np.zeros(fmri_dat.shape[:3])
f_brain[fmri_dat.sum(axis=3) > 0] = f.T
p = 1 - stats.f.cdf(f_brain[f_brain != 0], dfm, dfe)
print p
plt.figure(0, figsize=(8,8))
plt.hist(p, 20)
plt.title('Histogram of all p-values')
plt.xlabel('p-value')
plt.ylabel('Number of voxels');









    



[ 0.10880807  0.48019167  0.28623981  0.79338569  0.81082436  0.86541309
  0.20617917  0.85899707]



In [18]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal
import scipy.stats as stats
fmri_dat = np.random.rand(2,2,2,100)
voxel = fmri_dat[fmri_dat.sum(axis=3) > 0, :].T
wm_ts = voxel
t = voxel.shape[0]
lin_reg = np.array(range(0, t))
quad_reg = np.array(range(0, t))**2
wm_ts = wm_ts[:, wm_ts.std(axis=0) != 0]
wm_ts = signal.detrend(wm_ts, axis=0, type='linear')
wm_ts = wm_ts - wm_ts.mean(axis=0)
wm_ts = np.divide(wm_ts, wm_ts.std(axis=0))
U, s, V = np.linalg.svd(wm_ts)
wm_reg = U[:, 0:5]
R = np.column_stack((np.ones(t), lin_reg, quad_reg, wm_reg))
coefs = np.linalg.inv(R.T.dot(R)).dot(R.T).dot(voxel)
yhat = R.dot(coefs)
res = voxel - yhat
ssm = np.sum((yhat-voxel.mean(0))**2, axis=0)
sse = np.sum(res**2, axis=0)
dfm = coefs.shape[0]-1
dfe = voxel.shape[0]-coefs.shape[0]
msm = ssm/dfm
mse = sse/dfe
f = msm/mse
f_brain = np.zeros(fmri_dat.shape[:3])
f_brain[fmri_dat.sum(axis=3) > 0] = f.T
p = 1 - stats.f.cdf(f_brain[f_brain != 0], dfm, dfe)
print p
plt.figure(0, figsize=(8,8))
plt.hist(p, 20)
plt.title('Histogram of all p-values')
plt.xlabel('p-value')
plt.ylabel('Number of voxels');









    



[  1.11022302e-16   1.11022302e-16   1.11022302e-16   1.11022302e-16
   1.11022302e-16   1.11022302e-16   1.11022302e-16   1.11022302e-16]

Results

For the emotion task, we obtain the following histograms for the p values:



In [78]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
with open('./emotion/pvals_all.txt') as f:
    pvals = f.readlines()
pvals = [float(x.strip()) for x in pvals]
plt.figure(0, figsize=(8,8))
plt.hist(pvals, 20)
plt.title('Histogram of all p-values for Emotion Task')
plt.xlabel('p-value')
plt.ylabel('Number of voxels');
with open('./emotion/pwhite_mean.txt') as f:
    pwmean = f.readlines()
pwmean = [float(x.strip()) for x in pwmean]
with open('./emotion/pgrey_mean.txt') as f:
    pgmean = f.readlines()
pgmean = [float(x.strip()) for x in pgmean]
with open('./emotion/pother_mean.txt') as f:
    pomean = f.readlines()
pomean = [float(x.strip()) for x in pomean]
plt.figure(1, figsize=(8,8))
bins = np.linspace(0, 0.225, 10)
plt.hist(pwmean, bins, alpha=0.5, label='White')
plt.hist(pgmean, bins, alpha=0.5, label='Grey')
plt.hist(pomean, bins, alpha=0.5, label='Other')
plt.legend(loc='upper right')
plt.title('Histogram of mean white matter p-values for Emotion Task')
plt.xlabel('p-value')
plt.ylabel('Number of subjects');
subjects = ['101107', '102311', '102816', '103515', '106521', '108828', '111312', '111413', '113619', '116524', \
            '118932', '120515', '123925', '130013', '166438', '180129', '204016', '255639', '293748']
plt.figure(2, figsize=(16,80))
for i in range(len(subjects)):
    with open('./emotion/pwhite_' + subjects[i] + '.txt') as f:
        pw = f.readlines()
    pw = [float(x.strip()) for x in pw]
    with open('./emotion/pgrey_' + subjects[i] + '.txt') as f:
        pg = f.readlines()
    pg = [float(x.strip()) for x in pg]
    with open('./emotion/pother_' + subjects[i] + '.txt') as f:
        po = f.readlines()
    po = [float(x.strip()) for x in po]
    plt.subplot(len(subjects)/2+1, 2, i+1)
    bins = np.linspace(0, 1, 20)
    plt.hist(pw, bins, alpha=0.5, label='White')
    plt.hist(pg, bins, alpha=0.5, label='Grey')
    plt.hist(po, bins, alpha=0.5, label='Other')
    plt.legend(loc='upper right')
    title = 'Histogram of p-values for Subject ' + subjects[i]
    plt.title(title)
    plt.xlabel('p-value')
    plt.ylabel('Number of voxels');

By saving the F values as a nifti image, we can view them as a brain. Here is an example of Subject 101107's F values viewed as a brain:

For the motor task, we obtain the following histograms for the p values:



In [79]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
with open('./motor/pvals_all.txt') as f:
    pvals = f.readlines()
pvals = [float(x.strip()) for x in pvals]
plt.figure(0, figsize=(8,8))
plt.hist(pvals, 20)
plt.title('Histogram of all p-values for Motor Task')
plt.xlabel('p-value')
plt.ylabel('Number of voxels');
with open('./motor/pwhite_mean.txt') as f:
    pwmean = f.readlines()
pwmean = [float(x.strip()) for x in pwmean]
with open('./motor/pgrey_mean.txt') as f:
    pgmean = f.readlines()
pgmean = [float(x.strip()) for x in pgmean]
with open('./motor/pother_mean.txt') as f:
    pomean = f.readlines()
pomean = [float(x.strip()) for x in pomean]
plt.figure(1, figsize=(8,8))
bins = np.linspace(0, 0.06, 13)
plt.hist(pwmean, bins, alpha=0.5, label='White')
plt.hist(pgmean, bins, alpha=0.5, label='Grey')
plt.hist(pomean, bins, alpha=0.5, label='Other')
plt.legend(loc='upper right')
plt.title('Histogram of mean white matter p-values for Motor Task')
plt.xlabel('p-value')
plt.ylabel('Number of subjects');
subjects = ['101107', '102311', '102816', '103515', '106521', '108828', '111312', '111413', '113619', '116524', \
            '118932', '120515', '123925', '130013', '166438', '171431', '180129', '204016', '255639', '293748']
plt.figure(2, figsize=(16,80))
for i in range(len(subjects)):
    with open('./motor/pwhite_' + subjects[i] + '.txt') as f:
        pw = f.readlines()
    pw = [float(x.strip()) for x in pw]
    with open('./motor/pgrey_' + subjects[i] + '.txt') as f:
        pg = f.readlines()
    pg = [float(x.strip()) for x in pg]
    with open('./motor/pother_' + subjects[i] + '.txt') as f:
        po = f.readlines()
    po = [float(x.strip()) for x in po]
    plt.subplot(len(subjects)/2+1, 2, i+1)
    bins = np.linspace(0, 1, 20)
    plt.hist(pw, bins, alpha=0.5, label='White')
    plt.hist(pg, bins, alpha=0.5, label='Grey')
    plt.hist(po, bins, alpha=0.5, label='Other')
    plt.legend(loc='upper right')
    title = 'Histogram of p-values for Subject ' + subjects[i]
    plt.title(title)
    plt.xlabel('p-value')
    plt.ylabel('Number of voxels');

Here is an example of Subject 101107's F values viewed as a brain:

Discussion

While most of the mean p-values for the emotion task are above the specified significance of 0.05, when we look at the p-values voxelwise, the vast majority of the p-values are below the specified significance. For the motor task, both the majority of the mean p-values and the majority of the voxelwise p-values are below the specified significance. Since most values are below the significance, we should reject the null hypothesis in favor of the alternative hypothesis that the white matter model fits the data better than the intercept-only model. This suggests that many of the voxels in the data can be explained by white matter coefficients, which is problematic for the purposes of nuisance correction. However, some of the p-values were still above the significance, so using white matter principal components for nuisance correction could work acceptably in some cases.

Discriminability

We further evaluated the emotion task to compare the discriminability of the white matter method and a high pass filtering method for nuisance correction. We calculated discriminability in two different ways: the first in which all six task blocks were treated as separate scans of the same subject and the second in which every task block was split into two and each block was treated as a separate subject with two scans. We see the results below for the white matter method followed by the high pass filtering method. Note that an older version of the registration portion of the pipeline was used to obtain the processed data as the newer version was incompatible with the white matter compcor method.

White Matter Six Blocks:

White Matter Split Blocks:

High Pass Filtering Six Blocks:

High Pass Filtering Split Blocks:

We can see that the high pass filtering method of nuisance correction achieves better discriminability scores than the white matter compcor method.

Signal Correlation

White Matter CompCor Strategy

Below is code used to generate the HRF:



In [ ]:

    
import numpy as np
from nipy.modalities.fmri import hrf, utils
hrf_func = utils.lambdify_t(hrf.glover(utils.T))
tsteps = 136
frames = 176
e = np.true_divide(tsteps,frames)
e = np.multiply(e, frames-1)
t = np.linspace(0,e,frames)
h = hrf_func(t)*100
s = np.zeros(frames)
s[4]=1    # 3 sec
s[31]=1   # 24 sec
s[58]=1   # 45 sec
s[85]=1   # 66 sec
s[113]=1  # 87 sec
s[140]=1  # 108 sec
x = np.convolve(h,s)
x = x[:frames]
np.savetxt('./hs_emotion_wm.txt', x)

To determine correlation:



In [77]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
tsteps = 176
rois = 70
fnames = ['../ewok_week_0508/emotion_wm/101107_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/102311_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/102816_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_wm/103515_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/106521_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/108828_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_wm/111312_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/111413_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/113619_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_wm/116524_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/118932_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/120515_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_wm/123925_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/130013_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/166438_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_wm/180129_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/204016_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_wm/255639_3T_tfMRI_EMOTION_LR_desikan-2mm.txt',
          '../ewok_week_0508/emotion_wm/293748_3T_tfMRI_EMOTION_LR_desikan-2mm.txt']
x = np.zeros((len(fnames),tsteps,rois))
for i in range(x.shape[0]):
    x[i,:,:] = np.transpose(np.loadtxt(fnames[i]))
y = np.loadtxt('../ewok_week_0508/hs_emotion_wm.txt')
rsq = np.zeros((x.shape[0],rois))
for i in range(x.shape[0]):
    for j in range(x.shape[2]):
        coef = np.corrcoef(x[i,:,j],y)[0,1]
        rsq[i,j] = coef**2
rsq_avg_wm = np.zeros((1,rois))
for k in range(rsq.shape[1]):
    rsq_avg_wm[0,k] = np.mean(rsq[:,k])
print "Average R^2 over each ROI:"
print rsq_avg_wm
print "\nAverage R^2:"
print np.mean(rsq_avg_wm)









    



Average R^2 over each ROI:
[[ 0.01063135  0.00880775  0.01503879  0.01421322  0.01122874  0.00828013
   0.00580757  0.01659324  0.02011401  0.01693436  0.01611977  0.01445611
   0.0088344   0.01415333  0.00851566  0.01158036  0.02196428  0.01484749
   0.00847402  0.01162333  0.01500538  0.00710696  0.01219527  0.01594302
   0.01331847  0.02398743  0.01452153  0.00846787  0.01150786  0.0155653
   0.01178817  0.00749967  0.0202306   0.00969559  0.01119373  0.02047915
   0.01016982  0.0155005   0.01762364  0.00910465  0.01369283  0.02349894
   0.01318633  0.01965251  0.01163995  0.01578193  0.02246949  0.00617977
   0.01637042  0.00538566  0.01850068  0.0200847   0.01825059  0.01303813
   0.01374323  0.01149065  0.00820823  0.01855627  0.01528554  0.02853751
   0.02132405  0.00854699  0.01132234  0.02136579  0.01649634  0.01361173
   0.02007257  0.00783491  0.01005197  0.01251295]]

Average R^2:
0.0140830791465

ROIs and HRFs for example subject 101107:



In [69]:

    
tsteps = 136
frames = 176
e = np.true_divide(tsteps,frames)
e = np.multiply(e, frames-1)
t = np.linspace(0,e,frames)
plt.figure(figsize=(16,45))
for i in range(x.shape[2]):
    plt.subplot(14,5,i+1)
    plt.plot(t,y)
    plt.plot(t,x[0,:,i])
    plt.xlim(0,tsteps)

High Pass Filtering Strategy

Below is code used to generate the HRF:



In [ ]:

    
import numpy as np
from nipy.modalities.fmri import hrf, utils
hrf_func = utils.lambdify_t(hrf.glover(utils.T))
tsteps = 136
frames = 176
e = np.true_divide(tsteps,frames)
e = np.multiply(e, frames-1)
t = np.linspace(0,e,frames)
h = hrf_func(t)*100
h = h[2:]
s = np.zeros(frames-2)
s[2]=1    # 3 sec
s[29]=1   # 24 sec
s[56]=1   # 45 sec
s[83]=1   # 66 sec
s[111]=1  # 87 sec
s[138]=1  # 108 sec
x = np.convolve(h,s)
x = x[:frames-2]
np.savetxt('./hs_emotion_hpf.txt', x)

To determine correlation:



In [75]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
tsteps = 174
rois = 70
fnames = ['../ewok_week_0508/emotion_hpf/101107_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/102311_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/102816_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_hpf/103515_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/106521_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/108828_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_hpf/111312_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/111413_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/113619_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_hpf/116524_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/118932_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/120515_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_hpf/123925_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/130013_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/166438_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', \
          '../ewok_week_0508/emotion_hpf/180129_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/204016_3T_tfMRI_EMOTION_LR_desikan-2mm.txt', '../ewok_week_0508/emotion_hpf/255639_3T_tfMRI_EMOTION_LR_desikan-2mm.txt',
          '../ewok_week_0508/emotion_hpf/293748_3T_tfMRI_EMOTION_LR_desikan-2mm.txt']
x = np.zeros((len(fnames),tsteps,rois))
for i in range(x.shape[0]):
    x[i,:,:] = np.transpose(np.loadtxt(fnames[i]))
y = np.loadtxt('../ewok_week_0508/hs_emotion_hpf.txt')
rsq = np.zeros((x.shape[0],rois))
for i in range(x.shape[0]):
    for j in range(x.shape[2]):
        coef = np.corrcoef(x[i,:,j],y)[0,1]
        rsq[i,j] = coef**2
rsq_avg = np.zeros((1,rois))
for k in range(rsq.shape[1]):
    rsq_avg[0,k] = np.mean(rsq[:,k])
print "Average R^2 over each ROI:"
print rsq_avg
print "\nAverage R^2:"
print np.mean(rsq_avg)









    



Average R^2 over each ROI:
[[ 0.00978246  0.01138539  0.01202619  0.01572711  0.01102334  0.02106373
   0.01197759  0.01219902  0.02011171  0.01840148  0.01714934  0.02370187
   0.01641629  0.01271637  0.0135069   0.01688999  0.01456494  0.01817852
   0.01451592  0.01687498  0.02067427  0.01137753  0.01614305  0.01702569
   0.01587241  0.02416854  0.01572209  0.02234469  0.02267174  0.019989
   0.01420698  0.01627403  0.01674093  0.01695103  0.0111629   0.02575995
   0.01662949  0.01494446  0.02014949  0.01516628  0.02442706  0.01506209
   0.02072433  0.02514193  0.0117251   0.0111236   0.02798868  0.00768364
   0.01676975  0.01176701  0.01490756  0.02023952  0.02100153  0.01591012
   0.01179671  0.01430744  0.01444871  0.01687819  0.02181816  0.02696451
   0.02681924  0.00931821  0.02154081  0.02372     0.02304167  0.01354124
   0.02280658  0.01200017  0.00823826  0.01270696]]

Average R^2:
0.0169515213024

ROIs and HRFs for example subject 101107:



In [76]:

    
tsteps = 136
frames = 176
e = np.true_divide(tsteps,frames)
e = np.multiply(e, frames-1)
t = np.linspace(0,e,frames-2)
plt.figure(figsize=(16,45))
for i in range(x.shape[2]):
    plt.subplot(14,5,i+1)
    plt.plot(t,y)
    plt.plot(t,x[0,:,i])
    plt.xlim(0,tsteps)



In [79]:

    
print "Increase in R^2 value per ROI:"
print rsq_avg-rsq_avg_wm
print "\nAverage increase in R^2:"
print np.mean(rsq_avg-rsq_avg_wm)









    



Increase in R^2 value per ROI:
[[ -8.48887024e-04   2.57764662e-03  -3.01259703e-03   1.51388353e-03
   -2.05397163e-04   1.27836036e-02   6.17002412e-03  -4.39421819e-03
   -2.29308539e-06   1.46712190e-03   1.02956990e-03   9.24575793e-03
    7.58188824e-03  -1.43695929e-03   4.99123674e-03   5.30962669e-03
   -7.39933384e-03   3.33102922e-03   6.04190054e-03   5.25165025e-03
    5.66888971e-03   4.27057089e-03   3.94777999e-03   1.08266446e-03
    2.55393123e-03   1.81111873e-04   1.20056354e-03   1.38768161e-02
    1.11638788e-02   4.42369914e-03   2.41881420e-03   8.77436162e-03
   -3.48967488e-03   7.25544536e-03  -3.08344860e-05   5.28079757e-03
    6.45966390e-03  -5.56044202e-04   2.52585596e-03   6.06163377e-03
    1.07342305e-02  -8.43684473e-03   7.53799717e-03   5.48942166e-03
    8.51510122e-05  -4.65833687e-03   5.51918521e-03   1.50386495e-03
    3.99331275e-04   6.38135741e-03  -3.59311412e-03   1.54813731e-04
    2.75093456e-03   2.87199112e-03  -1.94652447e-03   2.81678943e-03
    6.24048255e-03  -1.67807922e-03   6.53261537e-03  -1.57300138e-03
    5.49518891e-03   7.71222801e-04   1.02184744e-02   2.35421446e-03
    6.54533424e-03  -7.04858304e-05   2.73400374e-03   4.16525470e-03
   -1.81370968e-03   1.94009755e-04]]

Average increase in R^2:
0.00286844215597

Based on the results of all of the above tests, we found it best to use high pass filtering instead of removal of white matter components for our pipeline's nuisance correction.