核磁的机器学习——Nilearn包的教程

最近在学习Nilearn包的使用，现将学习笔记和思考发布在这里，供大家参考。
以下训练的数据均来自于官网的练习数据 haxby2001的数据

一、导入核磁数据

这个没啥特别的，指定图像地址就好了，唯一要注意的是，需要是4D的nii.gz文件格式，这个可以通过其他软件转换。

#导入数据就是命名nii.gz的数据地址
file_name = r'D:\jupyter\nilearn\haxby2001\subj2\bold.nii.gz'

可选教程：如何查看图像

#查看数据的方法，用 plotting.view_img函数
from nilearn import plotting
plotting.view_img(file_name)

二、定义一个mask

这个定义的mask，就是ROI。
需要注意的是：是不是一定需要mask？
答案是：不是必须，但是建议。
定义mask的好处在于，能容易出结果，范围小了，也更适合讲故事了。
如果不要mask，那就是全脑的体素进行分析，那矩阵量就太大了，但是也可以分析，这是另一种方法，后面会介绍。这里先介绍简单的。

mask_filename = r'D:\jupyter\nilearn\haxby2001\subj2\mask4_vt.nii.gz'

可选教程：如何查看ROI

#如何查看ROI
anatomical_image = r'D:\jupyter\nilearn\haxby2001\subj2\anat.nii.gz'
# cmap：The colormap to use
plotting.plot_roi(mask_filename, bg_img=anatomical_image, cmap='Paired')

mask_filename的图像：
plotting.view_img(mask_filename)

anatomical_image的图像：
plotting.view_img(anatomical_image)

补充：
1、这个anatomical image图像，其实就是T1像，就是患者的结构像
2、plot_roi函数的用法：plot_roi(roi_img, bg_img=, cut_coords=None, output_file=None, display_mode=‘ortho’, figure=None, axes=None, title=None, annotate=True, draw_cross=True, black_bg=‘auto’, threshold=0.5, alpha=0.7, cmap=<matplotlib.colors.LinearSegmentedColormap object at 0x000002144F1BD490>, dim=‘auto’, vmin=None, vmax=None, resampling_interpolation=‘nearest’, view_type=‘continuous’, linewidths=2.5, **kwargs)

3、加载标签 Label

把标签存储在 CSV 文件中，以空格分隔，然后使用 Pandas 将它们加载到数组中。

import pandas as pd
# Load behavioral information
target = r'D:\jupyter\nilearn\haxby2001\subj2\labels.txt'
behavioral = pd.read_csv(target, delimiter=' ')
print(behavioral)

labels.txt文件：

conditions = behavioral['labels']

输出的conditions是panda的序列：pandas.core.series.Series，一共有1452个数字，刚好对应bold.nii.gz的1452个时间点。
这里我需要跟新手提一嘴，一般常规的功能像分析，是要去除前10个时间点的，所以使用自己的数据的时候，就需要和去除后的时间点对应上。

4、将分析限制为猫和人脸

上面的labels里面包含了很多conditions。因此，数据相当大。并非所有这些数据都对我们感兴趣进行解码，因此我们将仅保留与人脸或猫对应的fmri信号。我们创建属于条件的样本的掩码；然后将此掩码应用于fmri数据，以将分类限制为面部与猫的区分。

condition_mask = conditions.isin(['face', 'cat'])

pandas库Series的函数的isin函数（A的元素是否在B当中），输出True和False.
由于数据在一张大的4D图像中，我们需要使用index_img来轻松进行拆分。

from nilearn.image import index_img
fmri_niimgs = index_img(fmri_filename, condition_mask)

label_mask输出的是True和False，需要得到具体的label

conditions = conditions[condition_mask] #conditions剩下216行
# Convert to numpy array
conditions = conditions.values #把conditions的值提取出来，就是216个face或cat的字符串，numpy.ndarray格式

5、定义评估器

前面我们得到了目标label，得到了label对应的激活时间点，也得到了局部的ROI体素。首先，有一个很基础知识需要知道，就是在bold信号里，大脑某个体素值增高，就代表该区域被激活，反之就是抑制。
因此，事实上，以上这么多步骤，其目的就是：把某区域的体素的灰度值增高与降低同标签label进行建模预测。如果能预测成功，那这块区域激活或者抑制就与这类标签具有很强的关系，也就达到了我们的目的。
上面是废话，下面开始继续：
1.定义评估器：

from nilearn.decoding import Decoder
decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True)

standardize=True , 对每个体素上的时间序列数据做0-1归一化。
2.拟合数据

decoder.fit(fmri_niimgs, conditions)

3.检测预测

prediction = decoder.predict(fmri_niimgs)
print(prediction)

当然，这种流程的预测是没有意义的，原因是没有划分测试集和验证集，相当于提高告诉答案，然后再考试，准确率当然是1了。

6、使用交叉验证测量预测分数

让我们忽略训练期间的最后 30 个数据点，并测试对这最后 30 个点的预测：

fmri_niimgs_train = index_img(fmri_niimgs, slice(0, -30))
fmri_niimgs_test = index_img(fmri_niimgs, slice(-30, None))
conditions_train = conditions[:-30]
conditions_test = conditions[-30:]
decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True)
decoder.fit(fmri_niimgs_train, conditions_train)
prediction = decoder.predict(fmri_niimgs_test)
# The prediction accuracy is calculated on the test data: this is the accuracy
# of our model on examples it hasn't seen to examine how well the model perform
# in general.
print("Prediction Accuracy: {:.3f}".format((prediction == conditions_test).sum() / float(len(conditions_test))))

0.767的准确率看起来还不错，接下来需要进行手动的K折交叉验证，然后输出AUC曲线下面积来正确评估这个ROI mask来预测行为的准确率。

from sklearn.model_selection import KFold
cv = KFold(n_splits=5)
# The "cv" object's split method can now accept data and create a
# generator which can yield the splits.
fold = 0
for train, test in cv.split(conditions):fold += 1decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True)decoder.fit(index_img(fmri_niimgs, train), conditions[train])prediction = decoder.predict(index_img(fmri_niimgs, test))print("CV Fold {:01d} | Prediction Accuracy: {:.3f}".format(fold,(prediction == conditions[test]).sum() / float(len(conditions[test]))))

AUC曲线下面积的代码日后补充进来。

6.2 与解码器的交叉验证

这一节还是说交叉验证的事儿，上一节是手动进行K折交叉验证，这里是用评估器默认的功能来实现交叉验证。
解码器还默认实现交叉验证循环并返回一个形状数组（交叉验证参数，n_folds），其中将准确度定义为评分参数来使用准确度分数来衡量其性能。

n_folds = 5
decoder = Decoder(estimator='svc', mask=mask_filename,standardize=True, cv=n_folds,scoring='accuracy'
)
decoder.fit(fmri_niimgs, conditions)

以上是根据时间点来计算的，事实上，静息态功能磁共振在设计上，经常会有block设计，block就是组块的意思，通俗来讲，就是静息-任务-静息-任务-静息。。。这样的实验。
那么来说，对于这种组块设计的实验，其实应该加上组块的信息进行多个事件组的分析，也就是group（在这里应该是block）分析，其结果会更好点。
详细解释：
试想一下，把bold.nii.gz里面单个时间点当成1个人，那么就有216个人，我们去预测哪些人是cat，哪些人是face。
在传统临床上，一般是把216个人使用随机抽样的方法，变成训练集和验证集（有人会提到应该分3个集合，但是216的数据量还达不到那个程度），然后在训练集上使用K折交叉验证，最后用验证集去验证。
但是在核磁里，时间点是连续的，我们没办法确定具体这个时间点，这个患者的大脑状态就对应着cat或者face，时间分辨率就不是核磁的长处，所以这里就产生了误差。
所以更好的方法是，因为时间点是连续的，静息-任务-静息的block是设计好的，那么研究组间差异，会比研究单个时间点更有价值。即这里的216个人是按号码排好的，face和cat是固定出现的，那么我们去用单个block去机器学习，然后统计单个block的预测准确率，这种方法会更适用于核磁。所以这里的交叉验证方法也就相应的使用留一法，即留下一个block。但是要记住，对单个时间点进行留一法是非常错误的。
官网数据的block重复出现了12次，每次出现的时候face和cat包含了9个时间点。如果是自己的数据的话，block重复的次数肯定要2和更高，不然就会报错。
block留一法：

session_label = behavioral['chunks'][condition_mask]
from sklearn.model_selection import LeaveOneGroupOut
cv = LeaveOneGroupOut()
decoder = Decoder(estimator='svc', mask=mask_filename, standardize=True, cv=cv)
decoder.fit(fmri_niimgs, conditions, groups=session_label)
print(decoder.cv_scores_)

可以看出，在第八次的block中，预测有点问题，其他block都非常棒。
下面是12个block中，cat和face的分布：

暂时没看出第八次block有啥特别的地方，不负责盲猜可能和头动有一定关系。

8、检查模型权重，并转换成nii图像

coef_ = decoder.coef_
print(coef_)

这是一个 numpy 数组，每个体素只有一个系数，所以就是[1, 464]的shape。

coef_img = decoder.coef_img_['face']
decoder.coef_img_['face'].to_filename('haxby_svc_weights.nii.gz')

保存带有权重的nii图像

plotting.view_img(decoder.coef_img_['face'], bg_img=anatomical_image,title="SVM weights", dim=-1
)

使用虚拟分类器（dummy_classification）来作为底线（假设基线）检验上面的SVC分类器性能。

dummy_decoder = Decoder(estimator='dummy_classifier', mask=mask_filename, cv=cv)
dummy_decoder.fit(fmri_niimgs, conditions, groups=session_label)
# Now, we can compare these scores by simply taking a mean over folds
print(dummy_decoder.cv_scores_)

可以看出，虚拟分类器的预测率很低，这就是随便使用一个分类器能达到的底线预测率，而SVC确实提高了很多。
关于虚拟分类器（dummy_classification）的解释：
虚拟分类器为您提供了一种“基准”性能测量 - 即，即使只是猜测，也应该达到预期的成功率。
假设您想确定某个给定的对象是拥有还是不具有某种属性。如果您已经分析了大量这些对象并且发现90％包含目标属性，那么猜测该对象的每个未来实例都拥有目标属性，这使您有90％的正确猜测的可能性。以这种方式构建您的猜测等同于您在引用的文档中使用most_frequent方法。由于许多机器学习任务试图增加（例如）分类任务的成功率，所以评估基线成功率可以为分类器应该超出的最小值提供底限值。在上面讨论的假设中，您会希望分类器的准确率达到90％以上，因为即使是“虚拟”分类器，成功率也可达到90％。
如果使用上面讨论的数据训练带有stratified参数的虚拟分类器，则该分类器将预测其遇到的每个对象都有90％的概率拥有目标属性。这与训练具有most_frequent参数的虚拟分类器不同，因为后者会猜测未来对象具有目标属性。

理解以上，简单的Nilearn就算入门了，后续会继续更新我的学习笔记。我认为，这个方法还是很不错的，简单不复杂，而且契合现在机器学习的热潮，具有应用潜力，比如在我的领域，我准备用这个方法找分类针刺刺激的真假穴的脑区。

我学习完这个初步教程，我还有两个问题：
1、对于机器学习方法来说，使用原始未经处理的数据还是使用转换后的数据，对于结果预测率来说没有区别，重要的是，输入的变量有没有包含决定性变量。那么，在核磁领域，是经过slice timing、head motion、realign、segmentation等等步骤后的数据用来预测好点，还是原始数据好点？
个人见解：其实没太大区别，未处理的数据和处理后的数据，只要他们都是nii.gz格式的图像，那就能用上述方法分析，然后看那种分析准确率更高。
2、需不需要做slice timing，把时间点基线统一？
个人见解：这个等我通过学习了更高阶的nilearn方法，阅读了他们团队的文章后，再回答这个问题。
3、对于非block设计，如病人干预前后，那么1个病人有2个核磁状态，假如有100个病人，对于这种实验，该怎么使用nilearn去分类不同干预措施呢？
个人见解：等我后续学完更高阶的nilearn，学习完后，我再来更新这个问题。
提供一个思路：对于严格遵循随机对照方法的随机对照试验，因为其基线是差异很小的，如果这个时候，基线的差异小到可以忽略，那么就可以把干预后的结果作为label，然后抽样就在100个病人里随机抽，这样就回到了传统的机器学习方法。值得一试。

Over.
clancy wu