Nature Protocols：为解释神经成像中的机器学习模型建立一个统一的框架（下）

请点击上面“思影科技”四个字，选择关注我们，思影科技专注于脑影像数据处理，涵盖（fMRI,结构像,DTI,ASL,EEG/ERP,FNIRS,眼动）等，希望专业的内容可以给关注者带来帮助，欢迎留言讨论及转发推荐，也欢迎了解思影科技的课程及数据处理服务，可添加微信号siyingyxf或18983979082咨询(咨询电话[phone]，杨晓飞）。(文末点击浏览）

机器学习是一种用于创建连接大脑功能与行为之间关系的强大的计算工具，它在神经科学领域的应用越来越广泛。然而，这些模型非常复杂，通常很难解释，因此很难评估它们在神经科学上的有效性，以及它们对理解大脑的贡献。为了使基于神经影像的机器学习模型可以解释，它们应该（i）对人类是可理解的，（ii）提供有关特定大脑路径或区域中代表了哪些心理或行为构造的有用信息，并且（iii）证明它们是基于与之相关的神经生物学信号，而不是由干扰变量或伪影。因此，研究人员引入了一个统一的框架，该框架由模型级别、功能级别和生物学级别的评估组成，以提供补充结果，用来支持对模型工作方式和原因的理解。在昨天的分享中，我们已经理解了本文中所述的框架所针对的三个不同阶段的评估，今天我们带大家一起来看一下作者在该框架下进行的一个具体的示例，本次的示例方案详细的展现了具体操作过程以及代码，大家可以在阅读后到GitHub地址下载对应的代码，按照文中的步骤分步骤进行推演，以方便将这样一个系统的和有效的方案应用在自己的研究中。本文发表在Nature  Protocols杂志（可添加微信号siyingyxf或18983979082获取原文）。

流程

模型构建 耗时20分钟到几个小时

      1. 此步骤将构建基于fMRI的ML模型，该模型可预测目标结果。这是模型解释框架工作流程的前提步骤。尽管此步骤的细节不是当前方案的主要重点，但是本文简要描述了模型构建的过程，以提供有关此方案中使用的两种类型的模型的信息。选项A描述了广泛使用的线性算法SVM的训练。本文选择SVM是因为它是当前神经影像学文献中最受欢迎的ML算法之一，例如，根据作者在参考文献中（Woo C, Chang L J, Lindquist M A, et al. Building better biomarkers: brain models in translational neuroimaging[J]. Nature Neuroscience, 2017, 20(3): 365-377.）进行的调查，在神经影像学研究的481个ML模型中，有46.4％使用了SVM，然后是判别分析（discriminant analysis，占比12.7%）和逻辑回归（logistic regression，占比7.5%）。选项A中的步骤描述了如何使用CanlabCore工具（box 1）构建SVM模型。选项B是建立非线性模型的过程。作为非线性模型的示例，本文选择了卷积神经网络（CNN），这是一种成功的深度学习模型，适用于各种预测应用。选项B中的步骤描述了如何使用Keras构建CNN模型。

A. 训练SVMS

a) 准备一个数据矩阵。您可以使用以下几行代码从带有灰质mask的图像数据中导入功能磁共振成像数据。

b) (可选）使用mask。在训练模型之前，可以使用mask选择所需要的大脑特征。 例如，以下代码使用了Woo等人（Woo C, Koban L, Kross E, et al. Separate neural representations for physical pain and social rejection[J]. Nature Communications, 2014, 5(1): 5380-5380.）使用的mask。

c) 准备结果变量并将其添加到dat.Y。对于分类任务，使用分类变量（例如，不同条件的二进制值）定义结果向量，而对于回归分析，结果向量是连续变量（例如，trail-by-trail评分）。以下代码在示例分析中为SVM二进制分类器定义了结果向量。 具有[1，-1]值的编码值与Spider包兼容，应使用：

d) 定义训练集和测试集以进行交叉验证。为了获得对模型性能的无偏估计，应该使用交叉验证在内部对模型进行验证。交叉验证的关键步骤是选择一种将数据分为训练集和测试集的合适策略。可以将一定百分比的数据定义为测试集（k倍交叉验证）（例如，对10倍交叉验证使用10％的数据），也可以将一个被试的数据作为每一折的测试数据（留一交叉验证）。因为随着k的增加，预测精度的偏差会减少，方差会增加，因此需要在考虑可用数据集大小的情况下找到一种合适的折衷解决方案。留一交叉验证可用于小型数据集，而在大型数据集中，5折或10折交叉验证会更有效。注意不要在训练和测试数据之间引入任何依赖关系（见下文）。

       predict函数的可选输入参数'nfolds'可以指定将使用哪种类型的交叉验证。 输入k时，它将运行k倍交叉验证。可选输入也可以是用于使用自定义交叉验证（例如留一，leave-two-trials-out）的向量。如果可选输入等于1，则该函数将不会运行交叉验证。如果未提供任何可选输入，则默认情况下该函数将进行五折交叉验证。例如，作者可以在示例数据集中使用留一交叉验证。也就是说，这里将在每一折验证中保留一个被试的数据作为测试数据，而将其余58个受试者的数据用作训练数据。输入参数“ nfolds”定义如下：

关键步骤！保持训练集和测试集之间相互独立至关重要。在训练集和测试集之间执行一些分析（例如降噪，特征选择或缩放）可能会在训练集和测试数据集之间产生依赖性，从而导致预测性能的估算值出现偏差。如果被试之间呈现相关性，那么训练集和测试数据集之间也可能存在依赖关系，这在双胞胎研究中很常见，在其他一些研究中也存在，例如Human Connectome Project，ABCD和UK Biobank。此外，其他常见情况可能会引入某种依赖性。这些问题并没有明确地否定交叉验证，但是它们使得在独立的测试组中测试模型泛化性变得更加重要，并且最终使得在与训练样本在人群和采集程序上明显不同的数据群中进行泛化性测试变得更加重要。

关键步骤！如果将k折交叉验证与回归方法一起使用，作者建议使用分层交叉验证。要查看CanlabCore工具如何实现分层交叉验证，请参见stratified_holdout_set.m。

e) 用SVM模型拟合训练数据 作者使用predict函数，它是CanlabCore工具中fmri_data对象的一种方法（请参见box 2）。此函数使作者能够通过交叉验证轻松地对fMRI数据运行许多不同的ML算法。以下代码行将使SVM分类器拟合数据并测试分类器的交叉验证错误率。本文将在步骤2中说明交叉验证的详细信息。

输出变量stats提供多种信息，包括模型权重，基于交叉验证的y预测值（\hat{y}）和模型性能数值（有关更多详细信息，请参见box 2）。

       注意，可以使用嵌套的交叉验证来找到模型的最佳超参数。在嵌套交叉验证中，对训练集执行附加的嵌套交叉验证循环以调整超参数，而外部交叉验证循环用于估计模型性能。在本文的示例中，为简单起见，作者将函数默认设置为SVM算法，其中C = 1。

       如果在单个训练集上测试了具有多个不同参数的多个模型，则必须进行多重比较校正（例如错误发现率False Discovery Rate（FDR）校正）。请注意，针对多个比较的校正旨在防止出现假阳性特征识别，但不能避免由于在交叉验证循环之外测试多个模型而导致的过高的准确率。因此，一个更重要的步骤是在其他验证或独立测试数据集上测试最终的模型。

       在以下步骤中，本文将分析线性SVM模型（分类器模型）。尽管如此，预测函数还可以建立基于回归的模型来预测连续的结果。例如，下面的代码运行主成分回归（principal component regression，PCR）：

B. 训练CNN

a) 准备数据矩阵。作者使用pandas包中的DataFrame创建数据模版，使用Nilearn包从NIfTI文件中读取数据。这些在示例Jupyter Notebook中的‘Part 1: Initializing Data Matrix’和‘Part 2: Loading Data Function’均已实现。示例代码可以从https://github.com/cocoanlab/interpret_ml_neuroimaging/blob/master/scripts/cnn_lrp.ipynb获得。

b)定义一个CNN模型。使用Keras，可以使用以下代码定义CNN模型，该代码定义了四个卷积层和两个全连接层，然后是softmax输出层。 人们还可以使用Adam优化器定义损失函数以进行模型训练（cnn_lrp.ipynb的‘Part 3: CNN Model’）。

c) 将数据分为训练和测试集。如上所述，在线性SVM模型的情况下，定义用作交叉验证的训练数据和测试数据。作者使用scikit-learn创建了交叉验证测试框架，该框架在cnn_lrp.ipynb的‘Part 4: Model Training example’和‘Part 5: Leave-One-Subject-Out Cross-validation’ 中实现。

d) 在训练数据上拟合CNN模型。 可以使用小批量随机梯度下降方法（mini-batch stochastic gradient descent method）（例如Adam）的变体来训练CNN模型，该方法通过Keras软件包中的fit函数实现。在示例代码中，每次迭代使用的批量大小为32。运行此处定义的train_model函数时，作者使用Keras中的evaluate函数评估每个epoch的训练损失。以下代码还定义了训练CNN模型的epoch数（cnn_lrp.ipynb的‘Part 3. CNN Model’的第二节）：

e) 根据测试数据测试模型。CNN的训练完成后，作者可以按照下面的代码在单独的测试数据集上测试模型，这将返回测试集的预测准确率：

评估交叉验证性能 耗时5-20分钟关键！

       此步骤从准确性、敏感性和特异性方面评估新模型的预测性能。使用交叉验证或保留的测试数据来获得模型性能的无偏估计是至关重要的。下面，将提供留一交叉验证和八折交叉验证的分析示例。请注意，为方便起见，作者为k折交叉验证的折数k做出了任意选择，因为交叉验证通常对k的选择不敏感，但研究人员应考虑选择其他的交叉验证策略。

       2.使用训练集迭代拟合模型，在每一折的测试集上测试模型，并汇总每一折的预测结果值。在第1步中，作者将数据分为训练集和测试集。例如，将拆分数据集以执行留一交叉验证。因此，作者保留一个被试的数据作为测试集，并使用模型拟合其余58个被试的数据（在本文的示例中为SVM）。随后，使用该模型对留出来的数据进行分类，并为每个数据点保存其与超平面的距离。将这个过程重复59次。该过程在predict函数中是自动化的。

       此外，以下代码显示了如何运行八折分层交叉验证：

关键步骤！对于交叉验证的所有迭代，模型的超参数应保持相同。但是在嵌套交叉验证（一种选择最佳超参数的常用方法，具体流程如下图所示）中，可以针对每一折使用不同的超参数。

嵌套交叉验证流程图：内层交叉验证（inner loop）：用于模型选择，可以进行模型选择；外层交叉验证（outer loop）：k折交叉验证将数据分为训练集和测试集（注：图片来自https://ljalphabeta.gitbooks.io/python-/content/nested.html）

如果您对脑影像磁共振数据处理及机器学习感兴趣，欢迎浏览思影科技课程及数据处理服务（可添加微信号siyingyxf或18983979082咨询）：

3. 通过比较模型预测和实际结果来计算预测准确性、敏感性和特异性。在分类任务中，可以将精度定义为1减去分类错误率，分类错误率即分类错误的数据数量与数据总数之比。灵敏度定义为正确分类的阳性与阳性数据点数量之比，特异性定义为正确分类的阴性与阴性数据点数量之比。图3a示出了交叉验证的分类的准确性估计的结果。CanlabCore工具中的roc_plot函数可以计算所有这些性能指标：

示例输出：

注意！检查结果度量和模型响应的可靠性可以帮助确保模型性能无偏。通常，模型性能不能超过结果度量的可靠性。但是，这并不是严格意义上的，因为可靠性度量本身是错误估计的，并且估计的可靠性可能会受到除大脑与结果关系以外的变量的影响。尽管如此，结果的可靠性通常为任何预测模型提供了上界，因此可以认为是健全性检查。

图3 | SVM模型的预测性能（步骤2和3）以及线性模型的特征级评估结果（步骤7，选项A和B）。

a.图中显示了通过留一交叉验证验证的SVM模型的分类性能，错误分类的阈值设置为0(步骤2和步骤3)。Top panel显示了与超平面的交叉验证距离和决策阈值，bottom panel展示了ROC图。黄色的点表示分类正确，灰色的点表示分类错误。支持向量机模型的准确率达到92%±0.9%。

b.权重分布图显示了bootstrap测试识别出的SVM模型的显著特征权值，并以q（即FDR校正后的p） < 0.05的FDR进行阈值处理(步骤7，选项A)。c.权重分布图显示了在RFE过程后的最终预测的SVM特征，最终特征数= 20,000，每一步中去除的特征数= 5,000(步骤7，选项B)。

      分析模型是否受干扰变量或是伪影的影响 耗时5-20分钟关键！为了确保模型的响应不能被干扰变量或伪影所解释，研究者可以使用相关的干扰变量(例如，扫描仪头部的运动)来预测模型的响应。如果干扰变量和其他有害变量可以预测模型响应，则模型很可能受到它们的影响。同样，也可以检查感兴趣的结果是否与潜在的干扰因素和其他有害变量有关。以下步骤描述了如何执行此类分析：

4.创建一组有害变量（nuisance variables）（例如，六个trail-level平均纵向位移的运动参数(roll, pitch, yaw, x, y and z)）。作者也可以使用不同类型的干扰变量(例如，生理学干扰、来自独立成分分析的噪音成分)或TR-level或subject-level的数据。在这个示例分析中，使用了trail-level的移动参数。

5. 准备一个fmri_data对象变量（object variable）。应该将有害变量分配给obj.dat field，并将模型响应(在示例分析中，作者使用了数据点到超平面的距离)分配给obj.Y（输出） filed。

6. 使用有害数据（nuisance data）预测模型响应。在示例分析中，作者使用了多元回归，因为在这个示例中只有几个特征。在本例中，作者没有使用交叉验证，但是也可以使用交叉验证(例如，5折或10折交叉验证;请参阅第1步(iv))。模型的预测性能可以通过预测值(\hat{y})与真实值(y)之间的相关性来检验。

识别显著特征 耗时20分钟到几天

7. 识别显著特征。在这里，作者提出了四种评估特征重要性的方法：以自举测试为例评估线性模型中特征稳定性的方法；以RFE和“虚拟病变”分析为例评估线性模型中特征重要性的方法；以LRP（分层相关传播）为例评估非线性模型中特征重要性的方法。

自举测试识别对跨单元预测做出可重复(稳定)贡献的特征(例如，本例中的被试)。在以前的研究中已经成功地使用了自举测试。选项A中的步骤更详细地描述了自举测试。

与自举测试不同，RFE主要是一种嵌套的特征选择方法。模型被反复训练，同时通过一个消除步骤来定义一定数量的特征，在每次迭代中删除这些特征，直到满足一个停止条件为止(例如，当模型中到达一定数量的特征时)。研究者也可以根据最高的交叉验证预测精度选择最终的模型。选项B中的步骤描述了如何在线性模型中执行RFE。

在“虚拟损伤”分析中，可以通过以下方法来评估特征的重要性：

(i)当从模型中删除一组特征(区域或网络)时，预测性能下降了多少;

(ii)仅使用一组特征时，模型的性能如何。分析结果如图4所示，实例结果如表4所示。可以按照选项C中的步骤进行分析。

LRP(如选项D所述)的目标是分解最终的预测得分，并突出模型的每个输入特征对预测的贡献。LRP（分层相关传播）对于层次化的非线性模型如CNN特别有用，因为它可以有效地将分解的贡献从模型的输出层传播到输入特征层。LRP相比其他评估方法的一个独特之处在于，它可以对模型的每个预测类分别获得每个特征的重要性评估。此外，LRP计算的评估值可以是正的，也可以是负的，具体地表示每个特征是支持还是反对增加评估所计算的预测类的预测分数。尽管此方法最初用于解释对单个实例的预测，但在当前方案中，作者专注于组水平的解释，并以最终的解释作为单个解释的平均估计（图5）。

       A.自举测试

i) 数据重采样。自举测试的必要步骤是创建足够数量的新数据样本（即自举样本）。通过采样n个数据点并替换大小为n的数据集来创建m个自举样本。 例如，作者的数据集包含944个trail（即n = 944），作者随机抽取10,000个自举样本（m = 10,000），每个样本都包含944个数据点。

ii)使用每个自举样本训练模型。每个自举样本均用作新模型的训练数据。模型的超参数应与原始模型相同。此步骤在计算复杂度很高，因为它需要大量的迭代，具体取决于自举样本的数量。对于1,000个样本，最小P值为1 / 1,000或0.001。因此，当希望在分布的尾部具有合理的数值精度时，≥10,000个样本是理想的，这在神经成像中通常是这样。

iii)计算预测权重的P值。计算模型中two-tailed uncorrected P value，即权重在零上下附近的权重比例。另一种选择是使用样本分布的平均值和标准偏差来计算z和P值。

iv)进行多重比较校正。例如，用于多重比较的一种常见的校正选择是FDR校正，作者使用它来实现图3b所示的结果。校正后的阈值确定哪些特征始终有助于自举样本的预测。

      所有这些分析步骤都可以使用predict函数的“ bootweights”可选输入参数来执行。可选输入“ bootsamples”可用于设置自举样本的数量。默认值为100，但以前的研究中使用了5,000或10,000。

自举测试的结果存储在stats_boot.WTS中。另外，stats_boot.weight_obj还包含statistic_image对象中的自举测试结果，包括模型权重和P值。 statistic_image对象的threshold方法可用于执行多重比较校正（例如，q <0.05，FDR校正）。

关键步骤！如果数据量很大，则自举测试可能需要几天的时间。为了缩短时间，研究人员可以使用“useparallel”选项并行处理，或者，如果有多台计算机或节点可用，则可以使用可选输入“savebootweights”来保存来自多个自举样本的自举权重，并结合使用之后的结果。

B.RFE

i)选择RFE的参数。RFE的基本参数包括停止策略和消除步长。停止策略确定了重复消除和训练过程的终止时间。这取决于最终模型的期望特性。消除步长是每次迭代消除的特征数量。 消除步长的大小通常定义为其余特征的百分比，在先前的研究中有所不同。为了实现图3c所示的结果，作者使用了固定大小（5,000个特征）的消除步长。

ii)  使用预定义的模型参数在逐渐减少的特征数量上迭代训练模型。在每次迭代中，将一个新模型拟合到尚存的特征集中，并通过交叉验证评估其性能。交叉验证的结果应保存以供以后评估。按模型的权重的绝对值对它们进行排序，并消除与绝对值最小的权重相对应的某些特征。消除特征的数量由消除步长定义。在下一次迭代中，使用简化的特征集将新模型拟合到数据集。重复此过程，直到满足停止条件为止。

       所有这些分析步骤都可以通过svm_rfe函数执行，该函数对SVM执行RFE。在下面的示例中，输入data是一个包含训练数据的fmri_data对象，“ n_removal”选项指定了每次迭代中消除的特征数量，“ n_finalfeat”选项定义了停止策略。其余选项与predict函数相同。

C.“虚拟损伤”分析

i)准备划分mask（parcellation mask）。选择感兴趣的parcellation，并为parcellation准备mask。对于本方案，作者准备了一个mask。

ii)反复删除一组特征并测试模型。在每次迭代中，从完整模型中删除一部分网络或区域，并在测试简化模型时记录预测性能。

关键步骤！如果将训练数据集用于“虚拟损伤”分析，则该分析也应进行交叉验证。在每次迭代中，应从交叉验证中将mask应用于完整模型，并在每一折留出来的数据上测试简化模型。在收集所有交叉验证的预测后，可以计算简化特征后的模型的预测性能。

关键步骤！为了获得简化特征后的模型的预测，这里使用点积计算表示pattern expression values，但也可以使用其他相似性指标（例如Pearson相关性，Spearman相关性，余弦相似性）：

       其中，n是模型的体素数量，w是体素级别的预测权重，x是测试数据。 预测模型由跨所有体素的预测权重（\vec{w}）组成，用于指定激活的位置和模式。权重告诉作者如何将fMRI数据整合到单个预测中，然后将其用于分类测试或回归分析。

关键步骤！由于在删除模型的一部分后并未针对零阈值对模型进行优化，因此可以使用roc_plot函数的balanced precision threshold选项，该函数查找并使用最佳的阈值，该阈值可最大化平衡分类精度（各类别的平均准确性）。根据平衡精度输入关键参数“Optimal balanced error rate”到阈值。

iii) 迭代地保留一组用于预测的特征，并测试简化的模型。“虚拟病变”分析的第二个选项是删除目标特征集以外的所有特征（例如，属于特定静息态网络的特征集），然后共同测试保留的特征集的预测性能。

      其中，a_i^{l}是第i个结点的激活值。

      这个过程一直持续到输入层，然后计算每个输入特征的相关分数。通过对所有可用的输入数据进行上述过程的迭代，可以通过对单个解释进行平均得到特征级别相关性的组水平解释。此外，可以使用单样本t检验来检查哪些特征的相关性得分始终高于零或低于零。从t检验中获得P值后，需要使用FDR进行多重比较校正。

      要实现LRP（分层相关传播），可以使用与Keras兼容的innvate包（cnn_lrp.ipynb中的‘Part 7. Layer-wise Relevance Propagation’）。

      一旦重要的特征被识别出来，研究人员需要将它们可视化，看看它们是否有意义。用于模型可视化的方法和工具有很多，您可以选择适合自己需要的任何可视化工具。一般来说，对于线性模型，重要的特征可以在一个标准化的大脑空间中被可视化，例如MNI空间。在作者的例子中，作者使用以下MATLAB代码来可视化bootstrap测试(图3b)和RFE方法(图3c)识别的重要权重：

       使用orthviews函数，可以在标准MNI大脑图像上可视化图像。此外，使用带有可选输入“thresh”的write函数，可以将结果图像保存为NIfTI图像文件，然后可以与其他可视化工具一起使用。对于使用CNN模型的LRP结果，在图5a中显示了阈值化的平均相关分数。图5b展示了当输入是在拒绝条件（rejection condition）下获得的图像时预测拒绝条件的相关分数的可视化。

       为了能够创建具有印刷质量的脑图，CanlabCore工具箱提供了fmridisplay对象。例如,canlab_results_fmridisplay.m文件中的fmridisplay对象是非常有用的。此外，CanlabCore工具箱还提供了为fmri_data和region对象类创建surface图像、cluster_surf.m和surface()方法。在Python中，可以使用许多不同的可视化选项，包括plotting.plot_stat_map, plotting.plot_glass_brain和plotting.plot_surf_stat_map。这些函数的示例在Jupyter notebook上都有（见cnn_lrp.ipynb），可以在https://github.com/cocoanlab/interpret_ml_neuroimaging 查看相关代码。 

注：LRP（分层相关传播）理解起来可能有点困难，由于神经网络逐层传播的特性，不方便像其他线性模型那样直接可视化特征。这里采用的方法是逐层迭代地反向计算相关神经元的重要性得分来计算最终在输出层的特征的重要性得分，类似的案例可以查看相关文献（Hazlett H C, Gu H, Munsell B C, et al. Early brain development in infants at high risk for autism spectrum disorder[J]. Nature, 2017, 542(7641): 348-351.）。在这篇论文中，作者从婴儿的结构像数据中提取了若干个特征送入全连接网络中进行分类，然后通过类似LRP的方法计算得到对预测儿童自闭症比较重要的几个脑区特征。但是注意，文献中的计算方法可能和本文所讨论的方法不一致。如下图所示：首先通过神经元之间的连边权重计算神经元的贡献度，首先计算中间层中每个神经元的贡献度，然后选出这一层中比较重要的神经元F和G，然后再计算上一层中每个神经元与F、G之间的贡献度，选出其中比较重要的神经元A、C。如此，用户可以计算出在输入层中比较重要的特征输入。

（图片来自Hazlett H C, Gu H, Munsell B C, et al. Early brain development in infants at high risk for autism spectrum disorder[J]. Nature, 2017, 542(7641): 348-351.）

图4 |“虚拟损伤”分析示意图(步骤7，选项C)。“虚拟损伤”分析通过每次从模型中移除或使用一个区域或网络来调查单个区域或网络如何对模型的最终预测做出贡献。根据所选择的区域划分，对原始模型中的区域进行遮盖(要么删除一个区域，要么只使用一个区域进行预测)，并对遮盖后的模型性能进行评估。

表4“虚拟病灶”分析的一个例子

图5 | 逐层相关性传播结果（步骤7，选项D）。

a和b，作者运行LRP来解释每个被试和实验条件下每次trial的预测（程序的步骤7，选项D）。然后，作者计算了两种情况下被试的平均相关性得分。在a中，作者展示了受热痛苦的实验条件下的平均相关评分图，在未校正的P < 0.001处进行阈值处理，这相当于在q < 0.029处进行FDR处理。类似地，在b中，作者展示了在社会拒绝实验条件下的平均相关评分图，在未校正的P < 0.001处进行阈值处理(相当于FDR在q < 0.019处进行阈值处理)。

在独立数据集上进行泛化性测试 耗时1-2小时

关键步骤！此阶段测试了新模型对新个体、来自不同实验室、不同扫描仪以及不同的实验和测试条件的多个数据集上的泛化能力。泛化性测试是评估一个模型在各种不同环境变化下的鲁棒性的过程，因此这个阶段应该在不同的环境和测试条件下，对许多不同的数据集进行多次测试。为了说明这一点，作者在这里提供了一个例子来测试两个在先前研究中用于预测疼痛强度的先验预测模型的泛化性：NPS和SIIPS(图6a,c)。示例数据集用作独立的测试数据集。

8. 准备要测试的预测模型和测试数据。确保mask在MATLAB路径中。作者使用fmri_data对象读取contrast image数据。

9. 计算patter expression values。作者将apply_mask函数与“ pattern_expression”选项一起使用，以计算模型权重与大脑volume数据之间的点积。作者为每个被试和实验条件计算了一个pattern expression values，从而为每个实验条件获得59个pattern expression value。

10. 评估模型的预测性能，包括特异性和敏感性。在本例分析中，作者尝试使用NPS和SIIPS1区分热实验条件（heat condition）和温暖条件（warmth condition）以及排斥（rejection condition）和友好条件（friend condition）。假设每个人有一个值（pattern expression value），作者使用2AFC测试，它比较同一被试中两个对比条件下的pattern expression value，而不是使用阈值来判定。被试内数值更高的value被分类为positive condition。可以使用带有“twochoice”选项的roc_plot函数来执行2AFC测试。函数的输入是pattern expression value的向量和相应的二进制结果。结果如图6b,d所示。

图6 |泛化性测试(步骤8-10)。

a，超过FDR阈值q < 0.05的NPS预测权重图。

b，为了评估NPS(步骤8-10)的泛化性，作者使用不同实验条件下的信号模式权重（signature pattern weights）和激活图之间的点积来计算pattern expression values。然后，作者进行了2AFC测试，比较热（heat）与温暖（warmth）、排斥（rejection）与友善（friend）的对比实验。方框图显示了NPS对这四种实验条件的响应。方框之间的线描述了正确的(粉红色的线)和错误的(蓝色的线)分类。ROC曲线显示了区分热（heat）和温暖（warmth）条件(黄色)以及拒绝（rejection）和友善（friend）条件(蓝色)的敏感性和特异性。

c，超过FDR阈值q < 0.05的SIIPS预测权重图。

d，方框图显示了SIIPS1对这四种实验条件的响应。方框之间的线描述了正确的(粉红色的线)和错误的(蓝色的线)分类。ROC曲线显示了区分热（heat）和温暖（warmth）条件(黄色)以及拒绝（rejection）和友善（friend）条件(蓝色)的敏感性和特异性。

评价该模型在神经生物学上的合理性耗时1-2小时

11. 评价该模型在神经生物学上的合理性。该步骤旨在通过检查来自先前文献和更具侵入性的研究的证据来评估模型在神经生物学上的合理性和有效性。如上所述，可能无法对此生物学级别的评估提供明确的答案，而应将其视为一项开放性调查，需要来自不同学科、多模式和多学科的长期合作。检验模型在神经生物学上的合理性的一种方法是根据神经科学文献评估模型在特征级别和模型级别上的评估结果。因此，在当前的方案中，作者提供了两个生物级别评估的示例: 使用大规模的静息态功能网络评估模型的重叠部分，以及基于大规模元分析数据库Neurosynth的term-based解码。

选项A中的分析通过计算阈值化map或非阈值化map(预测权重或特征重要性)与功能网络之间的重叠(或模式相似度)来验证是哪个静息态功能网络在预测模型中扮演重要角色。对于图7a所示的分析结果，作者使用了来自Step 7a (bootstrap测试)的阈值预测权重图，并检查了它与七个大型功能网络的重叠部分。

按照选项B中的步骤执行基于元分析数据库的解码分析。

(A). 与大规模静息态功能网络的重叠

i)准备一个感兴趣的功能划分模版。在这里，作者使用fmri_data对象加载了具有七个静息态功能网络的唯一值的mask图像文件。在下一步的分析中，作者创建了一个指标矩阵，其中的维数等于体素的数量×网络的数量。

ii)根据特征显著性，准备阈值化的图像向量。对于本例，作者准备了基于步骤7A自举测试结果的阈值map。

iii)计算阈值模式map与每个网络之间的重叠比例。对于这种分析，可以使用canlab_pattern_similarity.m函数的“posterior_overlap”选项，它提供了观察给定每个网络的阈值map的后验概率。作者分别计算了正预测权重和负预测权重的基于重叠的相似度。

(B). 基于元分析数据库的解码分析

      i. 使用Neurosynth提供的大型元分析解码框架，可以识别与预测权重或特征重要性值的阈值化或非阈值map相关的心理术语。Neurosynth解码器使用为各种心理术语生成的元分析map，并评估它们与输入的大脑map的相似性。它返回一个具有输入和元分析map之间的相关系数的术语列表。研究人员可以使用Neurosynth Python包中的解码器(https://github.com/neurosynth/neurosynth/blob/master/neurosynth/analysis/decode.py)或Neurosynth解码器web应用程序(http://neurosynth.org/decode/)。这里，作者提供了一个用于解码分析的Python代码，图7b给出了非阈值化SVM模型的解码分析的示例结果。

图7 | 生物学级别评估（步骤11）和表征分析（步骤12-15）的示例。

a，雷达图描述了阈值化支持向量机模型与静息态功能网络之间观察重叠的后验概率(步骤11，选项a)。粉红色的图表示与模型正预测权值的重叠部分，蓝色的图表示与模型负预测权值的重叠部分。

b，柱状图显示了将Neurosynth解码器得到的函数项应用于非阈值化SVM模型(步骤11，选项b)，粉红色的柱状图表示正权值的解码结果，蓝色的柱状图表示负权值的解码结果。在表征分析(步骤12-15)中，作者比较了NPS和SIIPS对四种刺激条件的反应。作者首先使用相关系数比较两个准确率矩阵，然后使用网络和层次聚类方法可视化不同实验条件之间的关系。

执行表征相似性分析 耗时1-2小时

关键！这一步的目的是通过检查和比较多个实例和条件下的模型决策来阐明模型的表征和决策原则。有许多不同的方法可以实现这一步，但这里作者提供了一个关于两个疼痛的先天预测模型的表征相似性分析的例子，即NPS和SIIPS1。作者在示例数据集上测试了这些预测模型，该数据集包含四个实验条件的数据:疼痛（heat，pain）、温暖（warmth）、拒绝（rejection）和友善（friend）。然后，作者对每一对条件进行强制选择测试，以获得分类精度矩阵，然后将其作为距离度量(即，对于[A vs. B]，更高的分类精度意味着对于预测模型来说，实验条件A和实验条件B之间的表征距离相差甚远)。实例分析结果如图7c所示。

12. 得到不同模型的分类精度矩阵。使用步骤9获得的pattern expression value，对不同的实验条件pair之间执行强制选择分类测试。

13. 比较两个模型之间的分类精度模式(步骤13-15)。可以使用许多不同的方法来检查和比较多个条件下的模型表征。在这个方案中，作者首先使用两个向量化精度矩阵之间的皮尔逊相关。要做到这一点，运行以下代码:

在示例数据中，精度矩阵之间的相关性为0.71，表明模型之间的关系非常紧密（图7c）。

14. 使用MATLAB的图形分析工具可视化不同实验条件之间的关系。通过将分类精度从1减去，然后除以向量和，将精度矩阵转换为权重。用加权邻接矩阵绘制无向加权网络。网络图揭示了两种预测模型之间的一些差异——对于NPS，热实验条件（heat condition）与其他实验条件相差甚远，但是对于SIIPS1，热实验条件（heat condition）和拒绝实验条件（rejection condition）比NPS更接近(图7c)。要执行此步骤，使用以下代码：

15. 用single linkage method绘制树状图，以验证这些实验条件是如何分层聚类的(图7c)。考虑到示例数据集中实验条件的数量很少，聚类分析在这里可能不是很有用。然而，如果有更多的实验条件，聚类分析将变得更有用。要绘制这两个模型的树状图，请运行以下代码：

耗时

步骤1: 建立模型:20分钟到几个小时。A，SVMs: 20分钟至几个小时;B, CNN:几个小时

步骤2-3: 交叉验证的性能评估:5-20分钟

步骤4-6: 干扰变量分析：5-20分钟

步骤7: 确定重要特征:20分钟到几天。A，自举测试:1小时到几天; B, RFE: 1-6小时; C，“虚拟病灶”分析:20min ~ 1h; D, LRP: 20分钟到1小时

步骤8-10: 使用独立数据集测试泛化性:1-2小时

步骤11: 评价模型在神经生物学上的合理性:1-2 h

步骤12-15: 表征分析：1-2h

预期结果

当前方案提出了一种工作流程，该工作流程有望产生支持神经影像ML模型验证和解释的互补结果。该工作流是基于此处介绍的统一模型解释框架。

除了本过程中此处说明的方法外，作者还基于示例数据集提供了示例图和预期结果的可视化效果。图3a展示了使用留一交叉验证评估模型性能的示例结果，作者在步骤2和3中进行了介绍。在图3b，c和5中，通过自举测试，RFE和LRP（步骤7，选项A，B和D）分别识别出的重要特征并将其可视化。表4模拟了在步骤7（选项C）中描述的“虚拟病变”分析的示例结果。在步骤8-10中执行的泛化性测试的示例结果如图6b，d所示。生物学级别评估（步骤11）的结果如图7a，b所示。最后，图7c显示了在步骤12-15中所描述的表征相似性分析的可视化。

总结：

虽然基于神经成像的ML模型的解释是一个开放的过程，但本文中提出的方案将作为重要的一步，以发展可解释的和神经科学上可信的神经成像ML模型和生物标记。使用所提出的模型解释框架和精心设计的用于评估不同级别的模型测试。本方法将模型产生的若干证据汇聚在一起，从而构成对神经成像ML模型的整体解释和理解。对本文方案的具体学习以后，你会发现，在本框架下，以往分散在不同的优秀研究中被大量使用且被证明有效的方案都具体化为了详细的操作流程，并且这些流程的实现都已经有了现有且被验证过的工具的支持，希望类似的研究框架都可以形成一套系统的工作流程，从而帮助进行科学探索的研究者快速掌握可以解决复杂问题的有效工具。你还在等什么呢？快去下载了代码，在你的电脑跑起来吧。

如需原文及补充材料请加思影科技微信：siyingyxf 或者18983979082获取,如对思影课程及服务感兴趣也可加此微信号咨询。觉得对您的研究有帮助，请给个转发，以及右下角点击一下在看，是对思影科技莫大的支持。

微信扫码或者长按选择识别关注思影

非常感谢转发支持与推荐

欢迎浏览思影的数据处理业务及课程介绍。（请直接点击下文文字即可浏览思影科技所有的课程，欢迎添加微信号siyingyxf或18983979082进行咨询，所有课程均开放报名，报名后我们会第一时间联系，并保留已报名学员名额）：

第七届脑电信号数据处理提高班（南京，11.18-23）

脑磁图（MEG）数据处理学习班（预报名）

数据处理业务介绍：

招聘及产品：

招聘：脑影像数据处理工程师（重庆&南京）

此处点击在看，让更多朋友关注思影