pymvpa_basics

MVPAの計算をPyMVPA2.＊で行う時重要なポイントがいくつかある。(執筆中)

これはNeurodebian上のPyMVPA2.*をipythonで動かすことを全体に書かれています。マニュアルはここからダウンロードできます。

1)機能画像は4Dのファイルで1個であり、その中に時間的順序を追って全volumesが含まれること。

SPMの3D to 4D Conversion機能による4D化についてはpreprocessing(SPM12_CONN)の最後を良く読んでください。Smoothingをかけるデータの前の方がMVPAの精度は普通良いです。

2)機能画像とcoregisterする解剖学的画像は1次maskとして空間的解像度(voxelサイズ)が同じであること。低解像度だとメモリを要さないこと。

resliceに関しては、Nibabelはとても便利(その1: Dipyとともに脳画像の解像度変更)(Nibabel1)

を読んでください。初心者はMATLAB-SPMを使うと良いです(メモリも節約されます)。voxel sizeを構造・機能の双方の画像で[2 2 2]から[3 3 4]に落とした例をここに挙げておきます。

ただし、構造画像(低解像度のgrey matter maskであるl_gm.nii)は、CONNが出力したwc1cs20161020ys-0003-00001-000001-01.niiからipython, Nibabel, Dipyを用いて作りました(reslice.py)。これは中級者のテクニックです。ROIをまず考えない時はGrey matter maskを使うと良いでしょう。

3)chunksとtargetsが1行1volumeに対応する条件の時系列的変化に対応するテキストファイル(chunksTargets.txt)で、行数は全volumes数に一致すること。

これは以下のようなもので(抜粋)、条件とrun番号の組が1行となって、1 volumeに対応する。BOLD delayを考慮して時系列的に作っておき(これもPythonでスクリプトを作っておくと良い)、<Patient ID>/behaviouralに置く。この抜粋では1 という数字はfirst runを意味します。m-elefantのようなchunkの書き方は、それぞれのアイテムを区別するのと同時に、-(ハイフン)の前にある頭文字で、条件のカテゴリーを区別するような文字列操作を可能にします。

base 1

t-tape_measurec 1

rest 1

m-elephant 1

MVPAスクリプトの中の、データをロードする重要な部分は以下の通りである。これを抜粋する。

このWebだとPythonに必要な行のindentが復元されないので注意してください。

import mvpa2.suite as M

import numpy as N

import os # file path processing, directory contents

sessionPath = '/home/brain/work/20161020ys'#Path in Neurodebian

chunksTargets_boldDelay="chunksTargets.txt"

volAttribrutes = M.SampleAttributes(os.path.join(sessionPath,'behavioural',chunksTargets_boldDelay))

dataset = M.fmri_dataset(samples=os.path.join(sessionPath,'analyze/functional/l_4D.nii'),

targets=volAttribrutes.targets,

chunks=volAttribrutes.chunks,

mask=os.path.join(sessionPath,'analyze/structural/l_gm.nii'))

その後の計算の流れは以下の通りになります。

#あとは、データのdetrendをします。

M.poly_detrend(dataset, polyord=1, chunks_attr='chunks')

#z-scoreを計算します。

M.zscore(dataset, chunks_attr='chunks', param_est=('targets', ['base']))

という流れになります。

#次は1試行の中の複数のvolumesの平均化をします。

dataset = dataset.get_mapped(M.mean_group_sample(['chunks','targets']))

dataset.chunks = N.mod(N.arange(0,dataset.shape[0]),5)

#条件カテゴリーのラベルを抽出します。

dataset.targets = [t[0] for t in dataset.targets]

dataset = dataset[N.array([l in ['m', 't'] for l in dataset.sa.targets], dtype='bool')]

#以下が機械学習のアルゴリズムを集約している実装で、ANOVAで500個のvoxelsを素性選択(feature selection)し、PLR(penalized logistic regression)を使用している例です。

anovaSelectedPLR = M.FeatureSelectionClassifier(

M.PLR(),

M.SensitivityBasedFeatureSelection(

M.OneWayAnova(),

M.FixedNElementTailSelector(500, mode='select', tail='upper')

)

#以下が交差評価(cross-validation)の部分になります。leave one run outしています。

foldwiseCvedAnovaSelectedPLR = M.CrossValidation(

anovaSelectedPLR,

M.NFoldPartitioner(),

enable_ca=['samples_error','stats', 'calling_time','confusion']

)

#分類子を走らせ、精度を計算します。

results = foldwiseCvedAnovaSelectedPLR(dataset)

print 'accuracy',N.round(100-N.mean(results)*100,1),'%'

# 素性選択にANOVAを利用した場合、その情報を改めて分析し直す時のコードです。

numFeatures=500

sensana = anovaSelectedPLR.get_sensitivity_analyzer(postproc=M.maxofabs_sample())

cv_sensana = M.RepeatedMeasure(sensana, M.NFoldPartitioner())

sens = cv_sensana(dataset)

print sens.shape

M.map2nifti(dataset, N.mean(sens,0)).to_filename(os.path.join(sessionPath,'analyze/functional/',"anovaSensitivity_"+str(numFeatures)+'.nii'))

learnerWeightsAmalgFolds = sensana(dataset)

sum(learnerWeightsAmalgFolds.samples > 0)

learnerWeightedVoxels = dataset.fa.voxel_indices[(learnerWeightsAmalgFolds.samples > 0)[0]]

print 'learnerWeightedVoxels', learnerWeightedVoxels

N.savetxt('PLR-' + str(numFeatures) + '-learnerWeightedVoxels' + '.csv',learnerWeightedVoxels, delimiter=',', fmt='%d')

N.save('PLR-' + str(numFeatures) + '-learnerWeightedVoxels' + '.npy',learnerWeightedVoxels)

N.savetxt(os.path.join(sessionPath,'PLR-' + str(numFeatures) +'-learnerWeightedVoxels' + '.csv'),learnerWeightedVoxels, delimiter=',', fmt='%d')

N.save(os.path.join(sessionPath,'PLR-' + str(numFeatures) + '-learnerWeightedVoxels' + '.npy'),learnerWeightedVoxels)

Page updated

Google Sites

Report abuse