Multivariate_Searchlight

Searchlight

Kriegeskorteらが開発したSearchlightは、各voxelからの半径範囲(radius=0,1,2,3,...で差し渡したvoxels数指定)内でMVPAを行い、その精度を各voxelに与えて脳画像の形で出力できるようにする方法です。

PyMVPA2.*で実行するためには、以下のようなコードを書き加える必要があります。PyMVPA2.6のマニュアルを参照してください。

dataset=M.fmri_dataset(...)の内部、末尾に次の1行を加える。これはgrey matter全体でsearchlightを行うためですが(前の行のmaskと同じ構造が総が指定されています)、特定のROI maskに絞っても可能です。

add_fa={'vt_thr_glm': os.path.join(sessionPath,'analyze/structural/l_gm.nii')})# added for searchlight

そしてsearchlight用のコードをfoldwiseCvedAnovaSelectedPLRの後に、たとえば以下のように加えます。その他の必要な実装は

PyMVPA2.6のマニュアルの6.2.3 Minimal Searchlight Exampleを読めば分かります。

center_ids = dataset.fa.vt_thr_glm.nonzero()[0]

for radius in [3]:#radius is now 3. You may use 0, 1, or 2 instead.

    sl = sphere_searchlight(foldwiseCvedAnovaSelectedPLR, radius=radius, space='voxel_indices',

                            center_ids=center_ids,

                            postproc=mean_sample())

    ds = dataset.copy(deep=False,

                      sa=['targets', 'chunks'],

                      fa=['voxel_indices'],

                      a=['mapper'])

    sl_map = sl(ds)

    sl_map.samples *= -1

    sl_map.samples += 1

    niftiresults = map2nifti(sl_map, imghdr=dataset.a.imghdr)

    niftiresults.to_filename(os.path.join(sessionPath,'analyze/functional/result-searchlight.nii'))

出力されたresult-searchlight.niiは脳画像のように出力することが可能です。

それらの画像の各voxelの輝度は、その点における分類モデルの精度(%)に対応しております。それをz-score(平均0、標準偏差1)に計算し直して脳画像出力する場合は、以下のようにします。

import numpy as np

import pylab

import csv

import sys

import scipy

import nibabel

img = nibabel.load('result-searchlight.nii')

imgdata=img.get_data()

flattened_imgdata=np.ndarray.flatten(imgdata)

non_zero_data=[x for x in flattened_imgdata if not x==0] # Extract the voxels of which the accuracy is not 0.

print len(non_zero_data)

mean_non_zero_data=np.mean(non_zero_data)

print 'The mean is', mean_non_zero_data

sd_non_zero_data=np.std(non_zero_data)

print 'The standard deviation is', sd_non_zero_data

imgdatanew=imgdata

for i in range(np.shape(imgdata)[0]):

    for j in range(np.shape(imgdata)[1]):

        for k in range(np.shape(imgdata)[2]):

            if (imgdata[i,j,k]==0):

                imgdatanew[i,j,k]=0;

            else:

                tmp=(float)(imgdata[i,j,k]-mean_non_zero_data)/sd_non_zero_data

                if tmp>3.08: #P(z>3.08)=0.001 on the standard normal curve. Used a z-score cummulative table.

                    imgdatanew[i,j,k]=tmp

                else:

                    imgdatanew[i,j,k]=0

filename='univariate_result-searchlight.nii’

newimg = nibabel.Nifti1Image(imgdatanew, None, img.get_header()) 

newimg.to_filename(filename)

print newimg