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 研究的分步指南