SPM/参数化经验贝叶斯 (PEB)

概述

一个常见的实验目标是检验有效连接在受试者群体之间是否不同，或者在群体内部根据行为测量（例如测试分数）是否不同。一种方法是采用 DCM 连接参数并应用经典检验（例如 t 检验、CVA）。这种方法的缺点是它丢弃了关于连接强度的估计不确定性（方差）。或者，人们可以在参数之上构建一个贝叶斯层次模型 - 描述组水平效应如何在受试者基础上约束参数估计。SPM12 包含一个参数化经验贝叶斯 (PEB) 模型，它使评估组效应和受试者间连接参数的可变性成为可能。

本页描述了执行 DCM+PEB 分析的主要步骤。有关更详细的解释，请参阅教程论文：第 1 部分关于 fMRI 的 DCM 和第 2 部分关于 PEB，它们附带一个示例 fMRI 数据集。

第一级分析

指定一个示例 DCM 以供一个受试者使用

首先使用 SPM 的图形用户界面为一个示例受试者指定一个单一 DCM（您可以在 SPM 手册的第 35 章或上述教程论文中找到一个已完成的 fMRI 示例）。在这个模型中，确保打开所有感兴趣的连接。我们将称之为“完整”模型。

在受试者之间复制 DCM 并将它们组装到一个组 DCM (GCM) 文件中

接下来，在受试者之间复制此 DCM。您可以使用批处理来实现。

对于 fMRI

从主 SPM 窗口单击批处理以打开批处理编辑器。
单击 SPM -> DCM -> DCM 规范 -> fMRI 的 DCM -> 指定组

对于 M/EEG

从主 SPM 窗口单击批处理以打开批处理编辑器。
单击 SPM -> DCM -> DCM 规范 -> M/EEG 的 DCM

输出将是一个单元格数组，包含一列和每行一个受试者 (GCM_name.mat)。该数组中的条目可以是单个 DCM 的文件名，也可以是 DCM 结构本身。

估计模型

估计模型（即将它们拟合到数据，以获得参数估计）。这可以通过批处理轻松完成。

从主 SPM 窗口单击批处理以打开批处理编辑器。
单击 SPM -> DCM -> DCM 估计。
选择您在上面创建的 GCM 文件，其中包含一个包含所有受试者 DCM（或其文件名）的单元格数组。
单击输出，然后单击“覆盖现有 GCM/DCM 文件”。
单击绿色播放按钮。

以下是使用底层 Matlab 函数实现相同结果的方法

% Collate DCMs into a GCM file
GCM = {'DCM_subject1_full_model.mat'
       'DCM_subject2_full_model.mat'};

% Estimate each subject's model
GCM = spm_dcm_fit(GCM); 

% Alternatively, replace the line above with the following code to alternate between estimating
% DCMs and estimating group effects. This is slower, but can draw subjects out 
% of local optima towards the group mean.
% GCM = spm_dcm_peb_fit(GCM);

% Write results
save('GCM_example.mat','GCM');

第二级分析

估计一个二级 PEB (参数化经验贝叶斯) 模型

完成第一级分析后，我们现在在参数上创建一个二级（组）一般线性模型。

在批处理编辑器中选择 SPM -> DCM -> 二级 -> 指定/估计 PEB。
为分析命名并选择上面创建的 GCM 文件。
在协变量下，添加您想建模的任何组间效应。（这类似于在常规二级 SPM 分析中指定设计矩阵）。第一个回归量将建模组均值，并自动添加到设计矩阵中。一些 PEB 函数假设您添加的第一个协变量将具有实验意义，而后续协变量是无意义的干扰变量。如果您不填写协变量选项，则只建模组均值。
在字段下，选择在组级别建模 DCM 的哪些字段。我们建议将其限制在少量参数。例如，如果您对 A 矩阵和 B 矩阵感兴趣，请对每个矩阵运行一个单独的 PEB 分析。

这将创建一个名为 PEB_name.mat 的组级别模型。它包含代表每个组间效应对每个 DCM 连接的影响的参数。以下是使用底层 SPM 函数的等效 Matlab 代码

% Specify PEB model settings
% The 'all' option means the between-subject variability of each connection will 
% be estimated individually
M   = struct();
M.Q = 'all'; 

% Specify design matrix for N subjects. It should start with a constant column
M.X = ones(N,1);

% Choose field
field = {'A'};

% Estimate model
PEB     = spm_dcm_peb(GCM,M,field);

save('PEB_example.mat','PEB');

将完整的 PEB 模型与嵌套 PEB 模型进行比较以检验特定假设

上面创建的 PEB 模型将具有大量参数（组效应的数量乘以 DCM 连接的数量）。要检验假设，您可以将此完整 PEB 模型与某些参数关闭（固定在其先验均值为零）的嵌套 PEB 模型进行比较。例如，您可以关闭从区域 R1 到区域 R2 的连接上所有组水平效应，包括组均值。完整 PEB 与此嵌套 PEB 之间的模型证据差异，其中 R1->R2 参数已关闭，将告诉您 R1->R2 在组水平上的重要性。

要执行此分析，首先将每个假设实现为一个 DCM，其中某些连接已打开或关闭。例如，如果您有五个不同的连接假设，则应指定五个不同的 DCM。这些 DCM 不需要进行估计 - 它们只是为了告诉软件尝试关闭哪些参数。您可以使用 GUI 或编写脚本修改现有模型来实现，例如

% Get an existing model. We'll use the first subject's DCM as a template
DCM_full = GCM{1};

% IMPORTANT: If the model has already been estimated, clear out the old priors, or changes to DCM.a,b,c will be ignored
if isfield(DCM_full,'M')
    DCM_full = rmfield(DCM_full ,'M');
end

% Specify candidate models that differ in particular A-matrix connections, e.g.
DCM_model1 = DCM_full;
DCM_model1.a(1,2) = 0; % Switching off the connection from region 2 to region 1

DCM_model2 = DCM_full;
DCM_model2.a(3,4) = 0; % Switching off the connection from region 4 to region 3

无论你是使用图形界面还是脚本指定上述候选模型，都需要编写一个简单的脚本将它们整理成一个单元格数组，每个模型对应一列。

GCM = {DCM_full, DCM_model1, DCM_model2};
save('GCM_templates.mat', 'GCM');

回顾一下，现在你拥有一个包含组级参数的 PEB 文件，以及包含每个连接假设的 GCM 文件。现在，你可以比较不同模型的证据。使用批处理

在批处理编辑器中选择 SPM -> DCM -> 二级 -> 比较/平均 PEB 模型。
选择你在上一步中创建的 PEB 模型。
对于 DCM，选择 GCM_templates.mat 文件。
点击绿色播放按钮。

以下是等效的 Matlab 代码

% Compare nested PEB models. Decide which connections to switch off based on the 
% structure of each template DCM in the GCM cell array. This should have one row 
% and one column per candidate model (it doesn't need to be estimated).
BMA = spm_dcm_peb_bmc(PEB, GCM);

将创建两个窗口。一个标题为 BMC (“贝叶斯模型比较”)，显示比较完整 PEB 模型和嵌套 PEB 模型的结果。另一个窗口标题为 “PEB - 审查参数”，显示组级估计的连接强度，在 PEB（贝叶斯模型平均）上平均。

搜索嵌套 PEB 模型

除了比较特定假设，如上所述，你可能希望简单地从 PEB 中删除任何对模型证据没有贡献的参数。这种方法（以前称为事后搜索）可以通过以下步骤执行。

在批处理编辑器中选择 SPM -> DCM -> 二级 -> 搜索嵌套 PEB 模型。
选择之前创建的 PEB 文件和 GCM 文件，该文件包含一级 DCM。
按下绿色播放按钮。

以下是等效的 Matlab 代码

% Search over nested PEB models.
BMA = spm_dcm_peb_bmc(PEB);

同样，将创建两个窗口。一个标题为 BMC (“贝叶斯模型比较”)，显示被关闭的连接 - 包括组均值（跨受试者的共性）和前两个受试者间效应。另一个窗口标题为 “PEB - 审查参数”，显示估计的组级连接强度，在搜索识别的 PEB 模型（贝叶斯模型平均）上平均。

审查结果

要审查 BMA 分析（或 PEB 模型）的结果，请使用 PEB 审查功能

% Review results
spm_dcm_peb_review(BMA,GCM)

请注意，在这行代码中，我们为审查功能提供了组 DCM (GCM)，它提供了一些有用的信息，例如连接的名称。

解释输出

模型

回顾一下，PEB 模型是一个通用线性模型

$\theta ^{(1)}=X\theta ^{(2)}+\epsilon ^{(2)}$

其中 $\theta ^{(1)}$ 是所有受试者的所有 DCM 连接的所有连接参数的向量。设计矩阵 $X$ 包含每个连接的组平均强度的回归器（列），以及每个协变量对每个连接的影响。例如，对于特定连接，可能有一个回归器用于组平均，另一个回归器用于受试者年龄对该连接的影响。相应的参数 $\theta ^{(2)}$ 量化了每个协变量对每个连接的影响，这些影响是从数据中估计出来的。受试者间残差变异性（随机效应）构成向量 $\epsilon ^{(2)}$ 。

参数

估计的 PEB 参数 $\theta ^{(2)}$ 在 spm_dcm_peb_review 图形界面中显示为条形图，每个协变量对应一个条形图。更具体地说，当估计 PEB 模型时，将计算参数上的多元正态概率密度，期望值存储在字段 PEB.Ep 中，协方差矩阵存储在字段 PEB.Cp 中。条形的高度对应于期望值 (Ep)，误差线是 90% 贝叶斯置信区间（可信区间），从协方差矩阵 (Cp) 的主对角线计算得出。

中的一些参数 $\theta ^{(2)}$ 编码了跨受试者的共性。这些是受试者 DCM 连接的组平均值，或者基线连接，取决于设计矩阵中的协变量是否以均值为中心（见下面的示例设计矩阵）。正参数估计表明连接与协变量正相关。相反，负值表明连接与协变量之间存在负相关关系。

PEB 参数 $\theta ^{(2)}$ 的单位取决于底层 DCM 连接参数的单位。一些连接参数的单位为赫兹 (Hz)，因为它们是变化率（时间常数的倒数）。其他连接参数 - 那些具有正或负约束的参数，例如 fMRI 的 DCM 中的自连接 - 是无量纲的对数缩放参数。这些参数缩放了一些默认连接强度，该强度以 Hz 为单位。有关更多详细信息，请参阅正在使用的特定 DCM 模型的相关论文。对于 fMRI 的 DCM，可以在教程论文 - fMRI 的 DCM 第一部分中找到有关单位的详细解释。

后验概率

在 PEB 审查 GUI (spm_dcm_peb_review.m) 中，可以通过单击条形图中相关的条形来找到每个参数的概率。下拉菜单提供了两种计算这些概率的不同方法。

如果选择选项参数概率> 0，则会单独考虑每个参数，并且概率是根据（边缘）后验概率分布下有多少区域大于零来计算的。换句话说，概率与条形图上显示的粉红色误差条有关。请注意，这没有考虑参数之间的协方差。

或者，在审查贝叶斯模型比较或搜索的结果时，会显示选项自由能（有与无），如果可用，则应使用该选项。在对简化的 PEB 模型进行自动搜索时，这些概率是通过比较包含特定连接的所有模型的证据与不包含该连接的所有模型的证据（从搜索中获得的最佳 256 个模型中）来计算的。如果提供预定义模型，则使用相同的方法来计算概率，只是概率仅针对与共性和第一个组差异（协变量）相关的参数计算。

GUI 还提供了根据这些概率进行阈值的选项，下面将对此进行更详细的讨论。

显著性和阈值

在贝叶斯分析中没有显著性的概念 - 每个效应或模型都有一个后验概率。例如，您完全有权报告后验概率为 53%、86% 或 99% 的效应。

Kass 和 Raftery (1995) 建议对证据水平进行描述性标签，如下表所示。例如，考虑两个模型 m1 和 m2，它们的自由能分别为 F1=4 和 F2=7（自由能是对数模型证据的近似值）。然后，对数证据的比率，称为对数贝叶斯因子，就是自由能的差值：log BF = F1-F2 = 3。假设模型的先验分布是平坦的，这将对应于模型 m2 的后验概率为 95%（这可以使用每个模型的自由能的 softmax 函数来计算）。从表格中可以看出，对数贝叶斯因子为 3 被描述为“强证据”。因此，在论文中，您可以写“模型 m2 的后验概率为 95%，这对应于支持该模型的强证据”。

来自 Kass 和 Raftery (1995) 的证据水平类别
贝叶斯因子	对数（以 e 为底）贝叶斯因子	后验概率	证据水平
1 到 3	0 到 1	0.5 到 0.73	不值得一提
3 到 20	1 到 3	0.73 到 0.95	积极的
20 到 150	3 到 5	0.95 到 0.99	强
> 150	> 5	> 0.99	非常强

不需要根据效应的概率对效应进行阈值处理，因为贝叶斯分析没有显著性的概念。然而，如果有很多参数，重点关注最有可能的效应可能会有所帮助，例如，只绘制具有强证据的效应。PEB 审查 GUI 中提供了执行此操作的选项。然而，请记住，这仅仅是为了显示目的，而低于阈值的参数仍然可以对模型做出重要贡献。

留一法交叉验证

在 PEB 分析结果中识别出一个或多个组效应后，您可能希望询问对特定连接的效应是否足够大，足以预测新的受试者是否属于特定组，或者预测连续回归量，例如临床评分。您可以按照以下步骤执行此操作。

在批处理编辑器中选择 SPM -> DCM -> 二级 -> 预测（交叉验证）。
与步骤 3 一样指定 PEB 模型。这次，在字段下，您可能希望包含您之前发现表现出显著的组间效应的一个或两个连接。
按绿色播放按钮。

或者，使用 Matlab 代码

% Perform leave-one-out cross validation (GCM,M,field are as before)
spm_dcm_loo(GCM,M,field);

现在将估计 PEB 模型，同时排除一个受试者，并使用该模型根据所选择的特定连接来预测设计矩阵中的第一个组间效应（在常数列之后）。生成的图显示了每个受试者的预测组效应以及预测评分与已知评分之间的相关性。如果要预测的效应是二元的（例如患者或对照），则底部的图显示了在受试者基础上，可以对预测的组成员身份有多自信。如果它是连续变量，例如临床评分，则该图显示了预测准确率。

示例设计矩阵

这里，我们提供一些关于如何为一些典型实验设计定义组间设计矩阵 M.X 的示例。通用线性模型 (GLM) 的所有相同原则也适用于 PEB 方案，唯一的限制是 PEB 软件期望设计矩阵以一列 1 开头。对于每种类型的实验设计，有多种方法来编码设计矩阵。决定使用哪种设计的部分因素将基于可解释性 - 例如，您是否希望参数编码组之间的差异，还是分别编码每个组的平均值？此外，不同的设计将具有不同的统计效率。例如，如果回归量是正交的（统计独立），则效率将最大化。如果您对哪种设计最佳有任何疑问，请通过尝试不同的选项并选择给出最积极的自由能 PEB.F 的选项来执行贝叶斯模型比较。

一个组间因素（两组）

这可以用两个回归量来建模，分别编码 1) 组均值和 2) 相对于均值的组差异

${\begin{bmatrix}1&-1\\1&-1\\1&1\\1&1\end{bmatrix}}$

这里，行是四个受试者，列是回归量。前两个受试者属于组 1（在第二个回归量中用 -1 表示），后两个受试者属于组 2（在第二个回归量中用 1 表示）。请注意，为了使这种解释成立，第二个回归量必须具有零均值。如果组差异回归量没有零均值，则可以手动设置它

X(:,2:end) = X(:,2:end) - mean(X(:,2:end))

随后的 PEB 模型将包含参数，这些参数编码每个连接强度的组均值（由于第一个回归量）以及编码组差异导致的均值偏差的参数（由于回归量 2）。对于组差异，正估计参数表示组 2 的连接性强于组 1，负参数表示组 1 的连接性强于组 2。

或者，如果第二个回归量没有居中，则回归量为：1) 组 1 的连接性，以及 2) 组之间的差异

${\begin{bmatrix}1&0\\1&0\\1&1\\1&1\end{bmatrix}}$

组间因素（三组）

假设有三个被试组，每组两个被试。如果假设组间存在线性效应 - 即第1组 > 第2组 > 第3组 - 则可以使用两个回归变量来建模：1）总体均值，2）第1组和第3组之间的差异。

${\begin{bmatrix}1&1\\1&1\\1&0\\1&0\\1&-1\\1&-1\end{bmatrix}}$

或者，如果目标是测试共同点和组间差异（不指定线性效应），可以使用三个回归变量来编码：1）总体均值，2）第一组和第二组之间的差异，3）第二组和第三组之间的差异。

${\begin{bmatrix}1&-1&0\\1&-1&0\\1&1&-1\\1&1&-1\\1&0&1\\1&0&1\end{bmatrix}}$

或者，为了单独建模每一组，可以选择一组作为基线（这里，第1组），回归变量编码：1）第一组，2）相对于第一组，属于第二组的附加效应，3）相对于第一组，属于第三组的附加效应。

${\begin{bmatrix}1&0&0\\1&0&0\\1&1&0\\1&1&0\\1&0&1\\1&0&1\end{bmatrix}}$

非均衡设计 - 两个组间因素（三组）

考虑一个有三个组别的实验：A）接受药物的患者，B）接受安慰剂的患者，C）一个单一的控制组。这是一个非均衡设计的例子。虽然之前部分列出的设计可以使用，但使用以下三个回归变量可能更直观：1）基线（对照），2）成为患者的附加效应，3）药物和安慰剂之间的差异。每组 2 个被试。

${\begin{bmatrix}1&1&1\\1&1&1\\1&1&-1\\1&1&-1\\1&0&0\\1&0&0\end{bmatrix}}$

请注意，这些回归变量不是正交的，因此这将比接下来描述的完全均衡的析因设计效率更低。

均衡析因设计 - 两个组间因素（四组）

考虑一个均衡析因设计，有四个被试组：1）接受治疗的患者，2）接受安慰剂的患者，3）接受治疗的对照组，4）接受安慰剂的对照组。回归变量将编码：1）总体均值，2）成为患者的主要效应，3）治疗的主要效应，4）交互作用，即患者和对照组之间药物效应的差异。每组 2 个被试。

${\begin{bmatrix}1&1&1&1\\1&1&1&1\\1&1&-1&-1\\1&1&-1&-1\\1&-1&1&-1\\1&-1&1&-1\\1&-1&-1&1\\1&-1&-1&1\end{bmatrix}}$

这里，最后一列 - 交互作用 - 是通过两个主要效应的逐元素相乘生成的（两个主要效应必须先进行均值化）。请注意，此设计导致正交的回归变量 - 使析因设计达到最佳效率。

协变量

除了上述回归量外，通常还会有协变量，例如年龄或临床评分。这些可以添加到设计矩阵中，通常将它们（跨所有受试者）进行均值中心化是一个好主意，以使第一个回归量可以解释为均值。请注意，对于添加到设计矩阵中的每个协变量，每个 DCM 连接都会添加一个参数，因此最终可能会获得大量参数。因此，最好保持谨慎，将协变量数量保持在最低限度。

分层实验设计

考虑一个在运行间层级具有因子设计的实验。例如，我们将想象一项研究，其中两组受试者，一组接受药物，另一组接受安慰剂，每个受试者扫描两次：治疗前和治疗后。这是一个平衡的 2x2 设计，具有单元格或实验条件

第 1 组治疗前 (GCM1_pre)
第 1 组治疗后 (GCM1_post)
第 2 组治疗前 (GCM2_pre)
第 2 组治疗后 (GCM2_post)

我们假设为每个受试者的这四种实验条件中的每一种拟合了单独的 DCM。括号中的名称是存储每个条件下 DCM 的示例变量名称。每个变量都是一个单元格数组，每行一个受试者。以下是使用 PEB 对这种设计进行建模的两种方法。

选项 1：2 级设计（第一级：DCM，第二级：PEB）

创建一个 PEB 模型，其中设计矩阵 X 具有以下回归量

总体均值
组的主要效应（第 1 组为 -1，第 2 组为 1，然后进行均值校正）
治疗的主要效应（治疗前为 -1，治疗后为 1，然后进行均值校正）
组和治疗的交互作用（两个均值校正后的主要效应逐元素相乘）

并将此 PEB 模型拟合到跨受试者和时间点的所有 DCM（以匹配您创建的回归量的正确顺序），例如

GCM = {GCM1_pre; GCM1_post; GCM2_pre; GCM2_post};
PEB = spm_dcm_peb(GCM, X);

如果您希望比较预定义的替代 PEB 模型的证据，则会出现一个额外的技术考量，因为每个模型的证据只会在前两个回归量（在本例中，总体均值和组的主要效应）方面进行比较。要比较所有这 4 个因素的预定义模型，只需重新运行 PEB 分析，重新排序设计矩阵，使感兴趣的效应成为第二个回归量。例如，要调查哪个模型最能解释治疗的效果，请将列重新排序为：均值、治疗、组、交互，然后运行 spm_dcm_peb 和 spm_dcm_peb_bmc（或使用批处理）。

选项 2：3 级设计（第一级：DCM，第二级和第三级：PEB）

或者，可以使用 PEB-of-PEBs 方法，在 3 级层级中对相同设计进行建模，其中为每个组分别创建一个 PEB 模型。对于第 1 组，这将具有一个设计矩阵 X1，其中包含回归量

第 1 组的均值
治疗对第 1 组的影响（治疗前为 -1，治疗后为 1，然后进行均值校正）

这将被拟合到来自第 1 组的所有 DCM 中的一个列向量中

GCM1 = [GCM1_pre; GCM1_post];
PEB1 = spm_dcm_peb(GCM1, X1);

对于第 2 组，PEB 模型将具有一个设计矩阵 X2，其中包含回归量

第 2 组的均值
治疗对第 2 组的影响（治疗前为 -1，治疗后为 1，然后进行均值校正）

这将被拟合到第 2 组的 DCM 中

GCM2 = [GCM2_pre; GCM2_post];
PEB2 = spm_dcm_peb(GCM2, X2);

随后的 PEB 模型 PEB1 和 PEB2 将具有参数，这些参数编码每个 DCM 连接的两个效应（均值和治疗）。现在，将 PEB 模型的参数提升到层级的第三级，以识别组之间的共性和差异。设计矩阵 X3 将包含回归量

所有受试者的均值
组差异（第 1 组为 -1，第 2 组为 1，然后进行均值校正）

也就是说，设计矩阵的维数为 2x2，值为：[1 -1; 1 1]。这被拟合到第二级 PEB 参数中

PEBs = {PEB1; PEB2};
PEB3 = spm_dcm_peb(PEBs,X3);

最终模型第 3 级 PEB 模型，名为 PEB3，将包括每个第二级效应的组均值和组差异的参数。即，对于每个连接，将存在一个参数用于

总体平均连接性
治疗的主要效应
组的主要效应
组和治疗的交互作用

要执行贝叶斯模型简化并查看结果，请键入

BMA3 = spm_dcm_peb_bmc(PEB3);
spm_dcm_peb_review(BMA3);

注意：当提示 GCM 数组时，请单击取消按钮，以便 PEB-of-PEB 结果正确显示。在 SPM 的未来版本中，将解决此步骤的需要。

选择哪个选项？

选项 1 更简单，但选项 2 可能提供受试者间变异性的更好模型（随机效应）。为了评估哪个选项更好，您可以尝试两种选项，并比较自由能。为此，请减去每个模型的自由能：PEB.F - PEB3.F。正数表示选项 1 更好，负数表示选项 2 更好。请确保将您的发现反馈到 SPM 邮件列表中 - 这将有助于开发最佳实践。

进一步阅读

教程论文和使用 PEB 进行组 DCM 研究的分步指南