稀有二分類數據的 Meta 分析是醫學研究領域的難點，其方法學目前仍不成熟。傳統的 Meta 分析是基于正態-正態模型的固定效應分析或隨機效應分析，但這種方法仍然存在方法學上的問題。Stijnen 等提出一種基于廣義線性混合效應模型（generalized linear mixed model，GLMM）的確切研究內似然模型（exact within-study likelihood models，EWLM）Meta 分析技術，包括二項式-正態模型（binomial-normal，BN）和超幾何正態模型（hypergeometric-normal model，HNM）等，可方便地實現稀有二分類數據的隨機效應 Meta 分析。本文結合實例介紹其在 SAS 軟件中的實現方法，并提供相關的 SAS 代碼。
引用本文： 吳敏, 鄭建清, 黃碧芬, 肖麗華. 在廣義線性混合模型框架下稀有二分類數據的隨機效應 Meta 分析. 中國循證醫學雜志, 2019, 19(7): 875-882. doi: 10.7507/1672-2531.201902062 復制
在醫學領域內，通常將發生概率小于 1/100（1%）的事件定義為罕見事件，例如某些藥物治療的少見或罕見的不良反應[1]。由于罕見事件發生率低，感興趣的結局很少，單個研究結果通常不足以發現不同醫療干預措施效果的差異[2]。對于罕見事件，通過 Meta 分析匯總來自多個臨床試驗的數據提供了一種增加統計效能的方法，也是獲取稀有數據高質量證據的唯一方法，但目前仍存在方法學上的挑戰，因為現有大多數 Meta 分析技術是基于大樣本研究獲得的，對罕見事件計算不適用[3]。罕見事件最常用的標準 Meta 分析方法是采用逆方差法，包括固定效應模型和 DerSimonian-Laird 隨機效應模型[4, 5]，該方法計算每項研究的治療效應和標準誤差（SE）。對于二分類變量，治療效應最常用的效應指標是比值比（odds ratios，OR），但在單零研究或雙零研究中，因為會出現分母為零的情況，根據 OR 計算治療效應是不可能的。在這種情況下，通常采用“連續性校正”對零事件進行修正，但實踐發現，這種處理方法對罕見事件的 Meta 分析會造成結果偏倚[3]。
因此，對罕見事件，研究者需要探索更先進的統計方法。Stijnen 等[6]提出一種基于廣義線性混合效應模型（generalized linear mixed model，GLMM）的確切研究內似然模型（exact within-study likelihood models，EWLM）Meta 分析技術，包括二項式-正態模型（binomial–normal models，BNM）和超幾何正態模型（hypergeometric-normal model，HNM）等，可方便地實現稀有二分類數據的隨機效應 Meta 分析。目前國內尚沒有學者介紹 EWLM 在稀有二分類數據的隨機效應 Meta 分析中的應用，本文將詳細介紹該模型，并結合實例介紹其在 SAS 軟件中的實現方法，提供相關的 SAS 代碼。
1 稀有二分類數據 Meta 分析的近似正態-正態模型（normal-normal model，NNM）及 EWLM 簡介
Efthimiou 發表了一篇關于罕見事件的 Meta 分析實用指南的文章[2]，詳細介紹了目前可用的稀有二分類數據的 Meta 分析方法，讀者可參閱此文了解相關的知識。但該文章并未涉二項式-正態模型。本文簡單介紹隨機效應 Meta 分析的標準方法（傳統 NNM），并詳細介紹 Stijnen 等提出的 EWLM 法。
1.1 隨機效應 Meta 分析的標準方法：NNM
假設 
個獨立研究，分別提供了真值 
的估計值 
及其標準誤 
，其中 
通常是兩種治療之間的效果差異的量度。標準隨機效應的建模方法[4, 7, 8]如下：
 |
通過計算相關的效應值后，總體效應估計值是每個獨立研究的加權平均值：
 |
Stijnen 等將模型（1）和（2）稱為 NNM。公式 1 中的參數 θ 是感興趣的主要參數，稱為整體效應，而 σ2 描述了不同研究的真實效應的變異。統計學推斷方法如下：
 |
 |
普通最大似然法（maximum likelihood，ML）或約束最大似然法（restricted maximum likelihood，REML）均可用于估計參數[4, 9, 10]。但在正態假設下，REML 在統計上最優。因此，在本文中，NNM 方法默認首先使用 REML 進行擬合。
公式 2 中的標準誤是在假設 Si 已知且固定的情況下得出的，是 
真實標準誤的低估值[11]。B?hning 等[12]討論了由于不考慮 Si 的不確定性而導致的一般推論問題。
1.2 EWLM
1.2.1 用于比例 Meta 分析的 BNM
假設研究人員關心的結局變量是某個事件的發生，并且感興趣的參數是比例 π。例如 π 是某種干預（治療）下某種不良事件的發生率。每個研究報告了患有不良事件的患者數量 Yi 和治療患者的總數 Ni。標準方法是使用對數比值（log odds）
作為效應參數。
由 
估計，其標準誤 
[13]。
對于事件發生率很低或事件發生為 0 的情況，Hamza 等[14]建議采用二項式分布代替正態模型：
 |
然后，公式（3）中的正態研究內似然被二項式似然 
代替。Stijnen 等將模型（2）和（4）稱為 BNM。這是一個廣義線性混合效應模型。更具體地說，它是隨機截距邏輯回歸模型。
在 Hamza 等[14]報道的擴展模擬中，BNM 模型總是比 NNM 模型表現得更好，對于事件發生率接近 0 時，BNM 方法產生無偏的估計和合理的覆蓋概率。
需要注意的是 Proc NLMIXED 的統計基礎是基于 t 分布的，其自由度 df=N?1。使用 t 分布的原因是標準 Wald 方法沒有考慮 σ2 估計的不確定性。或者可使用似然比檢驗來測試 θ，以及用于計算可信區間的剖面似然性。后者是通過為 θ 指定一個固定值并搜索 θ 的兩個值來確定的。
1.2.2 用于 OR 的 Meta 分析的 HNM
對于雙臂研究，（對數）比值比是最常用的效應指標。令 
和 
表示對照組 i（i=1，...，N）中的事件數和患者數，
和 
為治療組中的事件數和患者數。相關的效應值及標準誤計算方法是：
 |
對于事件數=0 的情況，采用增加 0.5 進行連續性校正[15]。在實踐中，對于納入研究包含太多事件數為 0 的研究，最常用的簡單方法是改用固定效應模型，包括（精確）Mantel-Haenszel 程序或（精確）條件邏輯回歸等。然而，固定效應模型 Meta 分析假設研究之間具有同質性，但這種假設一般是很難成立的，而且它可能導致無效推理[16]。因此，隨機效應分析仍是優選的。在一般情況下，近似 NNM 方法中偏倚產生的主要原因是估計 
與其方差 
之間的相關性以及連續性校正的使用。估計 
與其方差 
之間的協方差可以寫成
 |
由于使用連續性校正仍然可能產生偏差，Stijnen 等建議用給定研究中事件總數 
的精確條件似然來代替公式 2 中的正態似然，即非中心超幾何分布的似然性。
 |
Stijnen 等將模型（2）和（7）稱為 HNM。HNM 模型不再是線性的，而是 GLMM。實際上，它是一種混合效應的條件邏輯模型。
如果事件的總數相對于組樣本較小，則非中心超幾何分布可以通過二項分布很好地近似。給定事件的總數 Yi，組 1 中的事件數大致具有二項分布，即
 |
這意味著可以使用帶有偏移變量 
的隨機截距邏輯回歸模型來實現模型。Stijnen 等將模型（2）和（8）稱為 BNM。采用這種模型在實踐中可能比模型（2）和（7）更可行，因為隨機效應邏輯回歸的程序比處理非中心超幾何分布的程序更多見。
2 資料與方法
2.1 示例數據介紹
為便于說明上述兩種模型在 SAS 軟件中的具體操作方法，我們采用《實用循證醫學方法學》（第 2 版）第 10 章實例 7 的數據為實例進行演示[17]。該數據集來源于 Brophy 等發表的關于“β 受體阻滯藥在充血性心力衰竭的作用”的系統評價[18]。為擬合本文提供的模型，我們增加了兩個變量“ttal、ctl”，分別代表治療組和對照組的樣本量（表 1）。
表1
待分析數據集
在數據集中，存在 2 項單零試驗，分別是 study 7、study 13；5 項雙零試驗，分別是 study 3、study 4、study 5、study 11、study 12。
將表 1 的數據錄入 SAS 軟件，并存儲在名為 Example_data 的數據集中（因篇幅限制，數據步的 SAS 代碼簡要輸入）。本文 Meta 分析的目的有兩個：即模擬實現單臂研究的比例的 Meta 分析和雙臂研究的比值比的 Meta 分析。實現單臂研究的比例的 Meta 分析，我們僅對數據集中的治療組部分進行分析。
2.2 SAS 代碼介紹
本文提供的代碼見附錄。包括數據集的錄入、采用連續性校正的方法對數據集中的單零研究或雙零研究進行校正，獲得數據集 Example_data_CR，然后使用 PROC MIXED 程序或 PROC NLMIXED 程序，實現 NNM 模型、BNM 模型、HNM 模型的單臂研究的比例或雙臂研究的比值比的 Meta 分析。
3 結果
3.1 單臂研究的比例的 Meta 分析結果
單臂研究的比例 Meta 分析采用比值（odd）作為效應指標。表 2 是單臂研究比例的 Meta 分析結果，包括了 Proc Mixed 程序的固定效應模型結果、Proc Mixed 程序的隨機效應模型結果、NLMIXED 程序擬合的 NNM 模型結果和 NLMIXED 程序擬合的 BNM 模型結果。表 3 是將各種模型提供的平均對數比值及其可信區間轉換回以事件發生率為標度的結果，即中位發生率及其可信區間。比值與事件發生率 π 的關系是 π=odds/(1+odds)。表 4 中的效應估計 θ 是平均 odds 的 logit 轉換結果，即 odds=exp（θ）。因此 π=exp（θ）/（exp（θ）+1）。
表2
不同模型的 Meta 分析結果
表3
不同模型的 Meta 分析的中位事件發生率結果
表4
不同模型的協方差參數估計結果
從表 2 的結果可以看出，利用 Proc MIXED 程序和 Proc NLMIXED 程序獲得的 NNM 模型結果基本相同，說明本文提供的兩種代碼是等效的，但 Proc NLMIXED 程序所需要的代碼更為簡潔，提示 Proc NLMIXED 程序操作更為簡便。
我們側重查看 NLMIXED 隨機效應模型（NNM 模型）和 NLMIXED 隨機效應模型（BNM 模型）兩種結果的差異。從表 2 數據可以看出，通過 BNM 獲得的平均對數比值比及其可信區間，與 NNM 模型具有差異，但 BNM 模型在沒有連續性校正的情況下（注意附錄代碼直接采用原始數據集進行分析），仍然可以將 20 個原始研究納入分析，資料利用率為 100%。從表 3 數據可以看出，BNM 模型給出中位發病率估計值比 NNM 的要小，因為 NNM 并未采用連續性校正。表 5 的資料可以看出，與 BNM 方法相比，NNM 方法大大低估了跨研究方差。最后需要注意的是，BNM 模型的平均對數比值的標準誤較大，這可以通過 NNM 模型不考慮 
的不確定性和更大的跨研究方差估計的事實來解釋。
表5
不同模型的 Meta 分析結果
需要注意的是，無論是 Proc MIXED 程序還是 Proc NLMIXED 程序，在 NNM 模型背景下，均需要對單零研究和雙零研究進行連續性校正，否則程序會刪除單零研究和雙零研究，甚至可能報錯。
3.2 雙臂研究的 OR 的 Meta 分析結果
上述結果模擬的是單臂研究結果，但臨床中更為普遍的是研究人員希望看到某種干預措施對研究對象的效應，此時研究同時包含試驗組和對照組，即雙臂研究。在雙臂研究中，最常用的效應指標是比值比。表 6 給出了示例數據 NNM、BNM 和 HNM 模型雙臂研究的 OR 值 Meta 分析結果。
表6
不同模型的 Meta 分析的 OR 結果
需要注意的是，本文提供的示例擁有罕見事件和很大治療效應。這意味著公式（9）右側的第一項將占主導地位，導致負協方差，即 
的估計值為 0 甚至為負值。運行本文提供的附錄代碼可以看出，由于負協方差的存在，基于 NNM 模型的隨機效應模型變得不可行。這點在基于 MIXED 隨機效應模型的結果得到驗證，即隨機效應模型提供的固定效應的解和固定效應模型提供的解是一致的。在實踐中，
的估計值為 0 并不罕見。在這種情況下，NNM 模型和 BN 模型均會簡化為固定效應模型[16]。這個觀點在本文提供的示例得到驗證。
HNM 和 NNM、BNM 方法之間存在顯著差異。由于 NNM 模型和 BNM 模型均退化為固定效應模型，這些方法提供的平均對數比值比的標準誤太小，可能是由于研究間差異較大以及它忽略了 
的變異性這一事實。而 HNM 模型仍然能夠有效實現隨機效應模型，結果提供了較大的標準誤和較寬的可信區間，而且得出了與 NNM、BNM 模型完全相反的統計結果，其 OR 估計值為 0.670，標準誤 SE=0.146［95%CI（0.443，1.014），P=0.057 5］。這是因為 HNM 模型無需考慮納入的研究是否存在零事件，并且無需進行連續性校正。Stijnen 等通過小型模擬研究，驗證了 NNM 和 HNM 模型的性能，結果證實 HNM 模型實際上是無偏的，并且 95% 可信區間的覆蓋百分比令人滿意；而 NNM 模型在對數優勢比遠離 0 時，偏倚很大，95% 可信區間的覆蓋率很差。因此 Stijnen 等認為，HNM 模型運行良好，而 NNM 表現令人不滿意[6]。
4 討論
用于罕見事件 Meta 分析的不同方法（模型）可能導致不同的結論，并且數據越稀疏，不同方法結果之間差異越大。傳統上，事件結局數據的 Meta 分析是通過摘要度量（summary measures）方法完成的，該方法需要利用近似正態的研究內似然（即本文中的 NNM 模型）。然而對于稀有數據，NNM 模型的正態性假設有時候根本不能成立。因此，研究者仍然需要尋找適合這類數據的 Meta 分析模型[2]。
在本文中，我們采用 Stijnen 等提出基于 GLMM 的 EWLM Meta 分析技術，結合相關的示例已經表明，精確的研究內似然性在稀有數據的 Meta 分析中有幾個優點：① 避免了傳統方法中因估計和標準誤差之間的相關性引起的偏差；② 考慮了標準誤估計的不確定性，因此獲得整體效應參數的更現實的標準誤；③ 避免使用（任意）連續性校正；④ 使隨機效應模型 Meta 分析在稀疏事件數據集中更可行，而在實踐中，大多數情況下進行了固定效應模型分析；⑤ 多變量擴展與傳統的匯總度量方法一樣簡單[6]。
本文示例提供的數據分析結果也證實了這一點。利用 NNM 模型和 BNM 模型得出的 Meta 分析結果具有很大的統計學效應，而 HNM 模型得出的結果則是非常接近傳統的“統計顯著性”閾值（即 P<0.05）；說明三種模型對結果的統計顯著性有影響。造成這種原因是因為由于本文提供的示例之間的協方差極度靠近 0，NNM 模型和 BNM 模型均簡化為固定效應模型，而 HNM 則仍然可以保持隨機效應模型，而實踐證明，對同一組數據，采用固定效應模型或隨機效應模型可能得出完全相反的統計學結果[16]。Bradburn 等[3]發現，對于稀有數據，目前常用的 Meta 分析方法會產生偏倚，偏倚最大的方法是 I-V 法（Inverse-variance）及 D-L 法（DerSimonian-Laird）計算 OR 或 RD，M-H 法（Mantel-Haenszel）零單元格加 0.5 校正法計算 OR。當事件發生率低于 1% 時，Peto 一步法計算 OR 偏倚最小，且統計效能最高，可以提供最佳的可信區間，無需考慮研究內治療組和對照組樣本量的不平衡問題，或治療效應異常大的情況。張天嵩等[17]在其一書中介紹了本文提供的示例，并通過 Peto 一步法排除了 3、4、5、11、12 項試驗的數據，最終 OR 估計值為 0.651（95%CI：0.572，0.740）。
張天嵩等[17]的做法是目前多數統計軟件的默認選擇，正如 Cochrane 手冊建議，此類研究并未包含相關的 OR/RR 的信息，因此應予以排除。然而，一些研究人員指出，從倫理的角度來看，雙零研究中的患者應該被納入分析[19]；而另一些人認為這些研究可能通過樣本量來提供相對治療效果的信息[20]。對于罕見事件的另一個問題是，隨機效應 Meta 分析可能會使效應的估計值產生偏差，繼而導致錯誤的可信區間。故對于目前常用的 Meta 分析方法仍應持謹慎態度。Efthimiou 等[2]對 Peto 一步法、Mantel-Haenszel 法、Logistic 回歸法[3]、貝葉斯方法[21]、相關響應變量的 Beta 二項式模型[22]等模型進行了詳細的優劣比較。正常情況下，Peto 方法獲得的結果與沒有采用連續性校正的 MH 類似但不完全相同，而如果 MH 方法采用連續性校正的情況下，結果則產生明顯偏倚。因此相對于校正的 MH 法，Peto 方法可能是更優秀的方法。固定效應貝葉斯模型給出了與 MH 方法（沒有連續性校正）幾乎相同的結果。在貝葉斯模型和 MH 方法中，隨機效應建模的影響相當小，因為異質性的估計值很小（正如本文筆者提供的示例）。相關響應變量的 β-二項式模型在某些情況下可能無法收斂。與其他方法相比，收斂性問題是該模型的潛在缺點之一。總之，基于 Efthimiou 等的研究結果，罕見事件的不同 Meta 分析方法可能導致不同的結論，并且數據越稀疏，替代方法的結果之間的差異越大[23]。對于稀有數據，正如本文預期的，Stijnen 等提出確切似然方法是特別有利的。但是，這并不意味著確切似然方法適用于所有的稀疏事件數據，畢竟該方法仍然依賴于漸近參數。一般來說，在數據非常稀疏的情況下，應該對隨機效應 Meta 分析持謹慎態度，因為隨機效應的正態模型更難以證明。而貝葉斯方法因其不是基于漸近參數并且考慮了模型參數中的所有不確定性，因此可能是稀疏事件數據中的有價值的替代方法。讀者可以參見相關理論研究[24, 25]。Stijnen 等在其論文的討論部分也簡單提出貝葉斯方法存在的問題[6]。
針對罕見事件進行 Meta 分析的實踐問題，有幾個問題值得重視：① 稀有數據，應避免使用人工連續性校正，也應該避免刪除單零研究和雙零研究；② 當事件很少發生時，系統評價者應始終使用一系列替代模型進行廣泛的敏感性分析，以確保其結果的穩健性。當結果對模型的選擇非常敏感時，研究人員應該對他們如何呈現和解釋他們的發現特別謹慎。在本文中，我們側重在介紹 Stijnen 等提出的模型方法。結合示例和 Stijnen 等提供的示例，本文提出的模型被證明在 SAS 軟件中可以很容易擬合，如混合效應（條件）邏輯回歸可能比通常的方法更合適。在示例中，單變量情況中沒有遇到收斂問題，而我們只有治療效應參數和跨研究方差。如果給出足夠寬且精細的參數起始值，則實際上保證了收斂。
總之，基于廣義線性混合效應模型的 EWLM Meta 分析技術，包括二項式-正態模型和超幾何正態模型，可方便地實現稀有二分類數據的隨機效應 Meta 分析。
附錄：SAS 代碼
SAS 代碼 1：資料錄入
data Example_data;
input study td ta ttl cd ca ctl;
datalines;
1 5 20 25 6 19 25
2 1 8 9 2 14 16
……
/*連續性校正*/
data Example_data_CR;
set Example_data;
if (td=0) then do;
td = td+0.5;
ttl = ttl+0.5;
end;
if (cd=0) then do;
cd = cd+0.5;
ctl = ctl+0.5;
end;
run;
SAS 代碼 2：單臂研究中基于 PROC MIXED 程序的 NNM 模型（包含固定效應模型和隨機效應模型）
/*計算效應值及其標準誤（est）*/
data Example_data_CR_study;
set Example_data_CR;
logodds=log(td/(ttl-td));
est=1/td+1/(ttl-td);
run;
/*比例的固定效應模型*/
Proc mixed method =ml data=Example_data_CR_study;
class study;
model logodds =/ s cl;
repeated /group=study;
parms / parmsdata=Example_data_CR_study
eqcons=1 to 22;
ods output CovParms=CovP;
run;
/*隨機效應模型*/
/*根據固定效應模型的協方差參數，設置隨機效應模型的協方差*/
data covp;
set covp;
drop _LABEL_;
drop CovParm Group ;
rename Estimate=est;
run;
data start;
input est;
cards;
0.01
run;
data CovP;
set start CovP;
run;
proc print; run;
Proc mixed cl method=ml data=Example_data_CR_study;
class study;
model logodds= / s cl;
random int/ subject=study s;
repeated /group=study;
parms / parmsdata = CovP
eqcons = 2 to 23;
run;
SAS 代碼 3：單臂研究中基于 PROC NLMIXED 程序的 NNM 模型
proc nlmixed data=Example_data_CR qpoints=100;
parms theta= -2.5 thetai= -2.5 sigmasq=0.5;
logodds=log(td/(ttl-td));
sesq=1/td+1/(ttl-td);
model logodds ~ normal(thetai,sesq);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 4：單臂研究中基于 PROC NLMIXED 程序的 BNM 模型
proc nlmixed data=Example_data qpoints=100;
parms theta= -2.5 to -0.5 by 0.5 sigmasq=0.5;
pi=exp(logodds)/(1+exp(logodds));
model td ~ binomial(ttl,pi);
random logodds ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 5：雙臂研究中基于 PROC NLMIXED 程序的 NNM 模型
proc nlmixed data=Example_data_CR qpoints=100;
parms theta= -0.4 sigmasq=0.2;
logodds=log((td/(ttl-td))/(cd/(ctl-cd)));
sesq=1/td+1/(ttl-td)+1/cd+1/(ctl-cd);
model logodds ~ normal(thetai,sesq);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 6：雙臂研究中基于 PROC NLMIXED 程序的 BNM 模型
proc nlmixed data=Example_data qpoints=100;
parms theta= -2.5 to -0.5 by 0.5 sigmasq=0.5;
yi=td+cd;
pi=exp(log(ttl/ctl)+thetai)/(1+exp(log(ttl/ctl)+thetai));
model td ~ binomial(yi,pi);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 7：雙臂研究中基于 PROC NLMIXED 程序的 HNM 模型
proc nlmixed data=Example_data qpoints=100;
parms theta=-0.4 sigmasq=0 to 0.5 by 0.05;
y=td+cd;
num = lgamma(ttl+1) - lgamma(td+1) - lgamma(ttl-td+1)
+ lgamma(ctl+1) - lgamma(cd+1) - lgamma(ctl-cd+1)
+ td*thetai;
den = 0;
lo=max(0,y-ctl); hi=min(y,ttl);
do j=lo to hi;
z = lgamma(ttl+1) - lgamma(j+1) - lgamma(ttl-j+1)
+ lgamma(ctl+1) - lgamma(y-j+1) - lgamma(ctl-y+j+1)
+ j*thetai;
den = den + exp(z);
end;
ll = num - log(den);
model td ~ general(ll);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
在醫學領域內，通常將發生概率小于 1/100（1%）的事件定義為罕見事件，例如某些藥物治療的少見或罕見的不良反應[1]。由于罕見事件發生率低，感興趣的結局很少，單個研究結果通常不足以發現不同醫療干預措施效果的差異[2]。對于罕見事件，通過 Meta 分析匯總來自多個臨床試驗的數據提供了一種增加統計效能的方法，也是獲取稀有數據高質量證據的唯一方法，但目前仍存在方法學上的挑戰，因為現有大多數 Meta 分析技術是基于大樣本研究獲得的，對罕見事件計算不適用[3]。罕見事件最常用的標準 Meta 分析方法是采用逆方差法，包括固定效應模型和 DerSimonian-Laird 隨機效應模型[4, 5]，該方法計算每項研究的治療效應和標準誤差（SE）。對于二分類變量，治療效應最常用的效應指標是比值比（odds ratios，OR），但在單零研究或雙零研究中，因為會出現分母為零的情況，根據 OR 計算治療效應是不可能的。在這種情況下，通常采用“連續性校正”對零事件進行修正，但實踐發現，這種處理方法對罕見事件的 Meta 分析會造成結果偏倚[3]。
因此，對罕見事件，研究者需要探索更先進的統計方法。Stijnen 等[6]提出一種基于廣義線性混合效應模型（generalized linear mixed model，GLMM）的確切研究內似然模型（exact within-study likelihood models，EWLM）Meta 分析技術，包括二項式-正態模型（binomial–normal models，BNM）和超幾何正態模型（hypergeometric-normal model，HNM）等，可方便地實現稀有二分類數據的隨機效應 Meta 分析。目前國內尚沒有學者介紹 EWLM 在稀有二分類數據的隨機效應 Meta 分析中的應用，本文將詳細介紹該模型，并結合實例介紹其在 SAS 軟件中的實現方法，提供相關的 SAS 代碼。
1 稀有二分類數據 Meta 分析的近似正態-正態模型（normal-normal model，NNM）及 EWLM 簡介
Efthimiou 發表了一篇關于罕見事件的 Meta 分析實用指南的文章[2]，詳細介紹了目前可用的稀有二分類數據的 Meta 分析方法，讀者可參閱此文了解相關的知識。但該文章并未涉二項式-正態模型。本文簡單介紹隨機效應 Meta 分析的標準方法（傳統 NNM），并詳細介紹 Stijnen 等提出的 EWLM 法。
1.1 隨機效應 Meta 分析的標準方法：NNM
假設 個獨立研究，分別提供了真值 的估計值 及其標準誤 ，其中 通常是兩種治療之間的效果差異的量度。標準隨機效應的建模方法[4, 7, 8]如下：
|
通過計算相關的效應值后，總體效應估計值是每個獨立研究的加權平均值：
|
Stijnen 等將模型（1）和（2）稱為 NNM。公式 1 中的參數 θ 是感興趣的主要參數，稱為整體效應，而 σ2 描述了不同研究的真實效應的變異。統計學推斷方法如下：
|
|
普通最大似然法（maximum likelihood，ML）或約束最大似然法（restricted maximum likelihood，REML）均可用于估計參數[4, 9, 10]。但在正態假設下，REML 在統計上最優。因此，在本文中，NNM 方法默認首先使用 REML 進行擬合。
公式 2 中的標準誤是在假設 Si 已知且固定的情況下得出的，是 真實標準誤的低估值[11]。B?hning 等[12]討論了由于不考慮 Si 的不確定性而導致的一般推論問題。
1.2 EWLM
1.2.1 用于比例 Meta 分析的 BNM
假設研究人員關心的結局變量是某個事件的發生，并且感興趣的參數是比例 π。例如 π 是某種干預（治療）下某種不良事件的發生率。每個研究報告了患有不良事件的患者數量 Yi 和治療患者的總數 Ni。標準方法是使用對數比值（log odds） 作為效應參數。 由 估計，其標準誤 [13]。
對于事件發生率很低或事件發生為 0 的情況，Hamza 等[14]建議采用二項式分布代替正態模型：
|
然后，公式（3）中的正態研究內似然被二項式似然 代替。Stijnen 等將模型（2）和（4）稱為 BNM。這是一個廣義線性混合效應模型。更具體地說，它是隨機截距邏輯回歸模型。
在 Hamza 等[14]報道的擴展模擬中，BNM 模型總是比 NNM 模型表現得更好，對于事件發生率接近 0 時，BNM 方法產生無偏的估計和合理的覆蓋概率。
需要注意的是 Proc NLMIXED 的統計基礎是基于 t 分布的，其自由度 df=N?1。使用 t 分布的原因是標準 Wald 方法沒有考慮 σ2 估計的不確定性。或者可使用似然比檢驗來測試 θ，以及用于計算可信區間的剖面似然性。后者是通過為 θ 指定一個固定值并搜索 θ 的兩個值來確定的。
1.2.2 用于 OR 的 Meta 分析的 HNM
對于雙臂研究，（對數）比值比是最常用的效應指標。令 和 表示對照組 i（i=1，...，N）中的事件數和患者數， 和 為治療組中的事件數和患者數。相關的效應值及標準誤計算方法是：
|
對于事件數=0 的情況，采用增加 0.5 進行連續性校正[15]。在實踐中，對于納入研究包含太多事件數為 0 的研究，最常用的簡單方法是改用固定效應模型，包括（精確）Mantel-Haenszel 程序或（精確）條件邏輯回歸等。然而，固定效應模型 Meta 分析假設研究之間具有同質性，但這種假設一般是很難成立的，而且它可能導致無效推理[16]。因此，隨機效應分析仍是優選的。在一般情況下，近似 NNM 方法中偏倚產生的主要原因是估計 與其方差 之間的相關性以及連續性校正的使用。估計 與其方差 之間的協方差可以寫成
|
由于使用連續性校正仍然可能產生偏差，Stijnen 等建議用給定研究中事件總數 的精確條件似然來代替公式 2 中的正態似然，即非中心超幾何分布的似然性。
|
Stijnen 等將模型（2）和（7）稱為 HNM。HNM 模型不再是線性的，而是 GLMM。實際上，它是一種混合效應的條件邏輯模型。
如果事件的總數相對于組樣本較小，則非中心超幾何分布可以通過二項分布很好地近似。給定事件的總數 Yi，組 1 中的事件數大致具有二項分布，即
|
這意味著可以使用帶有偏移變量 的隨機截距邏輯回歸模型來實現模型。Stijnen 等將模型（2）和（8）稱為 BNM。采用這種模型在實踐中可能比模型（2）和（7）更可行，因為隨機效應邏輯回歸的程序比處理非中心超幾何分布的程序更多見。
2 資料與方法
2.1 示例數據介紹
為便于說明上述兩種模型在 SAS 軟件中的具體操作方法，我們采用《實用循證醫學方法學》（第 2 版）第 10 章實例 7 的數據為實例進行演示[17]。該數據集來源于 Brophy 等發表的關于“β 受體阻滯藥在充血性心力衰竭的作用”的系統評價[18]。為擬合本文提供的模型，我們增加了兩個變量“ttal、ctl”，分別代表治療組和對照組的樣本量（表 1）。
表1
待分析數據集
在數據集中，存在 2 項單零試驗，分別是 study 7、study 13；5 項雙零試驗，分別是 study 3、study 4、study 5、study 11、study 12。
將表 1 的數據錄入 SAS 軟件，并存儲在名為 Example_data 的數據集中（因篇幅限制，數據步的 SAS 代碼簡要輸入）。本文 Meta 分析的目的有兩個：即模擬實現單臂研究的比例的 Meta 分析和雙臂研究的比值比的 Meta 分析。實現單臂研究的比例的 Meta 分析，我們僅對數據集中的治療組部分進行分析。
2.2 SAS 代碼介紹
本文提供的代碼見附錄。包括數據集的錄入、采用連續性校正的方法對數據集中的單零研究或雙零研究進行校正，獲得數據集 Example_data_CR，然后使用 PROC MIXED 程序或 PROC NLMIXED 程序，實現 NNM 模型、BNM 模型、HNM 模型的單臂研究的比例或雙臂研究的比值比的 Meta 分析。
3 結果
3.1 單臂研究的比例的 Meta 分析結果
單臂研究的比例 Meta 分析采用比值（odd）作為效應指標。表 2 是單臂研究比例的 Meta 分析結果，包括了 Proc Mixed 程序的固定效應模型結果、Proc Mixed 程序的隨機效應模型結果、NLMIXED 程序擬合的 NNM 模型結果和 NLMIXED 程序擬合的 BNM 模型結果。表 3 是將各種模型提供的平均對數比值及其可信區間轉換回以事件發生率為標度的結果，即中位發生率及其可信區間。比值與事件發生率 π 的關系是 π=odds/(1+odds)。表 4 中的效應估計 θ 是平均 odds 的 logit 轉換結果，即 odds=exp（θ）。因此 π=exp（θ）/（exp（θ）+1）。
表2
不同模型的 Meta 分析結果
表3
不同模型的 Meta 分析的中位事件發生率結果
表4
不同模型的協方差參數估計結果
從表 2 的結果可以看出，利用 Proc MIXED 程序和 Proc NLMIXED 程序獲得的 NNM 模型結果基本相同，說明本文提供的兩種代碼是等效的，但 Proc NLMIXED 程序所需要的代碼更為簡潔，提示 Proc NLMIXED 程序操作更為簡便。
我們側重查看 NLMIXED 隨機效應模型（NNM 模型）和 NLMIXED 隨機效應模型（BNM 模型）兩種結果的差異。從表 2 數據可以看出，通過 BNM 獲得的平均對數比值比及其可信區間，與 NNM 模型具有差異，但 BNM 模型在沒有連續性校正的情況下（注意附錄代碼直接采用原始數據集進行分析），仍然可以將 20 個原始研究納入分析，資料利用率為 100%。從表 3 數據可以看出，BNM 模型給出中位發病率估計值比 NNM 的要小，因為 NNM 并未采用連續性校正。表 5 的資料可以看出，與 BNM 方法相比，NNM 方法大大低估了跨研究方差。最后需要注意的是，BNM 模型的平均對數比值的標準誤較大，這可以通過 NNM 模型不考慮 的不確定性和更大的跨研究方差估計的事實來解釋。
表5
不同模型的 Meta 分析結果
需要注意的是，無論是 Proc MIXED 程序還是 Proc NLMIXED 程序，在 NNM 模型背景下，均需要對單零研究和雙零研究進行連續性校正，否則程序會刪除單零研究和雙零研究，甚至可能報錯。
3.2 雙臂研究的 OR 的 Meta 分析結果
上述結果模擬的是單臂研究結果，但臨床中更為普遍的是研究人員希望看到某種干預措施對研究對象的效應，此時研究同時包含試驗組和對照組，即雙臂研究。在雙臂研究中，最常用的效應指標是比值比。表 6 給出了示例數據 NNM、BNM 和 HNM 模型雙臂研究的 OR 值 Meta 分析結果。
表6
不同模型的 Meta 分析的 OR 結果
需要注意的是，本文提供的示例擁有罕見事件和很大治療效應。這意味著公式（9）右側的第一項將占主導地位，導致負協方差，即 的估計值為 0 甚至為負值。運行本文提供的附錄代碼可以看出，由于負協方差的存在，基于 NNM 模型的隨機效應模型變得不可行。這點在基于 MIXED 隨機效應模型的結果得到驗證，即隨機效應模型提供的固定效應的解和固定效應模型提供的解是一致的。在實踐中， 的估計值為 0 并不罕見。在這種情況下，NNM 模型和 BN 模型均會簡化為固定效應模型[16]。這個觀點在本文提供的示例得到驗證。
HNM 和 NNM、BNM 方法之間存在顯著差異。由于 NNM 模型和 BNM 模型均退化為固定效應模型，這些方法提供的平均對數比值比的標準誤太小，可能是由于研究間差異較大以及它忽略了 的變異性這一事實。而 HNM 模型仍然能夠有效實現隨機效應模型，結果提供了較大的標準誤和較寬的可信區間，而且得出了與 NNM、BNM 模型完全相反的統計結果，其 OR 估計值為 0.670，標準誤 SE=0.146［95%CI（0.443，1.014），P=0.057 5］。這是因為 HNM 模型無需考慮納入的研究是否存在零事件，并且無需進行連續性校正。Stijnen 等通過小型模擬研究，驗證了 NNM 和 HNM 模型的性能，結果證實 HNM 模型實際上是無偏的，并且 95% 可信區間的覆蓋百分比令人滿意；而 NNM 模型在對數優勢比遠離 0 時，偏倚很大，95% 可信區間的覆蓋率很差。因此 Stijnen 等認為，HNM 模型運行良好，而 NNM 表現令人不滿意[6]。
4 討論
用于罕見事件 Meta 分析的不同方法（模型）可能導致不同的結論，并且數據越稀疏，不同方法結果之間差異越大。傳統上，事件結局數據的 Meta 分析是通過摘要度量（summary measures）方法完成的，該方法需要利用近似正態的研究內似然（即本文中的 NNM 模型）。然而對于稀有數據，NNM 模型的正態性假設有時候根本不能成立。因此，研究者仍然需要尋找適合這類數據的 Meta 分析模型[2]。
在本文中，我們采用 Stijnen 等提出基于 GLMM 的 EWLM Meta 分析技術，結合相關的示例已經表明，精確的研究內似然性在稀有數據的 Meta 分析中有幾個優點：① 避免了傳統方法中因估計和標準誤差之間的相關性引起的偏差；② 考慮了標準誤估計的不確定性，因此獲得整體效應參數的更現實的標準誤；③ 避免使用（任意）連續性校正；④ 使隨機效應模型 Meta 分析在稀疏事件數據集中更可行，而在實踐中，大多數情況下進行了固定效應模型分析；⑤ 多變量擴展與傳統的匯總度量方法一樣簡單[6]。
本文示例提供的數據分析結果也證實了這一點。利用 NNM 模型和 BNM 模型得出的 Meta 分析結果具有很大的統計學效應，而 HNM 模型得出的結果則是非常接近傳統的“統計顯著性”閾值（即 P<0.05）；說明三種模型對結果的統計顯著性有影響。造成這種原因是因為由于本文提供的示例之間的協方差極度靠近 0，NNM 模型和 BNM 模型均簡化為固定效應模型，而 HNM 則仍然可以保持隨機效應模型，而實踐證明，對同一組數據，采用固定效應模型或隨機效應模型可能得出完全相反的統計學結果[16]。Bradburn 等[3]發現，對于稀有數據，目前常用的 Meta 分析方法會產生偏倚，偏倚最大的方法是 I-V 法（Inverse-variance）及 D-L 法（DerSimonian-Laird）計算 OR 或 RD，M-H 法（Mantel-Haenszel）零單元格加 0.5 校正法計算 OR。當事件發生率低于 1% 時，Peto 一步法計算 OR 偏倚最小，且統計效能最高，可以提供最佳的可信區間，無需考慮研究內治療組和對照組樣本量的不平衡問題，或治療效應異常大的情況。張天嵩等[17]在其一書中介紹了本文提供的示例，并通過 Peto 一步法排除了 3、4、5、11、12 項試驗的數據，最終 OR 估計值為 0.651（95%CI：0.572，0.740）。
張天嵩等[17]的做法是目前多數統計軟件的默認選擇，正如 Cochrane 手冊建議，此類研究并未包含相關的 OR/RR 的信息，因此應予以排除。然而，一些研究人員指出，從倫理的角度來看，雙零研究中的患者應該被納入分析[19]；而另一些人認為這些研究可能通過樣本量來提供相對治療效果的信息[20]。對于罕見事件的另一個問題是，隨機效應 Meta 分析可能會使效應的估計值產生偏差，繼而導致錯誤的可信區間。故對于目前常用的 Meta 分析方法仍應持謹慎態度。Efthimiou 等[2]對 Peto 一步法、Mantel-Haenszel 法、Logistic 回歸法[3]、貝葉斯方法[21]、相關響應變量的 Beta 二項式模型[22]等模型進行了詳細的優劣比較。正常情況下，Peto 方法獲得的結果與沒有采用連續性校正的 MH 類似但不完全相同，而如果 MH 方法采用連續性校正的情況下，結果則產生明顯偏倚。因此相對于校正的 MH 法，Peto 方法可能是更優秀的方法。固定效應貝葉斯模型給出了與 MH 方法（沒有連續性校正）幾乎相同的結果。在貝葉斯模型和 MH 方法中，隨機效應建模的影響相當小，因為異質性的估計值很小（正如本文筆者提供的示例）。相關響應變量的 β-二項式模型在某些情況下可能無法收斂。與其他方法相比，收斂性問題是該模型的潛在缺點之一。總之，基于 Efthimiou 等的研究結果，罕見事件的不同 Meta 分析方法可能導致不同的結論，并且數據越稀疏，替代方法的結果之間的差異越大[23]。對于稀有數據，正如本文預期的，Stijnen 等提出確切似然方法是特別有利的。但是，這并不意味著確切似然方法適用于所有的稀疏事件數據，畢竟該方法仍然依賴于漸近參數。一般來說，在數據非常稀疏的情況下，應該對隨機效應 Meta 分析持謹慎態度，因為隨機效應的正態模型更難以證明。而貝葉斯方法因其不是基于漸近參數并且考慮了模型參數中的所有不確定性，因此可能是稀疏事件數據中的有價值的替代方法。讀者可以參見相關理論研究[24, 25]。Stijnen 等在其論文的討論部分也簡單提出貝葉斯方法存在的問題[6]。
針對罕見事件進行 Meta 分析的實踐問題，有幾個問題值得重視：① 稀有數據，應避免使用人工連續性校正，也應該避免刪除單零研究和雙零研究；② 當事件很少發生時，系統評價者應始終使用一系列替代模型進行廣泛的敏感性分析，以確保其結果的穩健性。當結果對模型的選擇非常敏感時，研究人員應該對他們如何呈現和解釋他們的發現特別謹慎。在本文中，我們側重在介紹 Stijnen 等提出的模型方法。結合示例和 Stijnen 等提供的示例，本文提出的模型被證明在 SAS 軟件中可以很容易擬合，如混合效應（條件）邏輯回歸可能比通常的方法更合適。在示例中，單變量情況中沒有遇到收斂問題，而我們只有治療效應參數和跨研究方差。如果給出足夠寬且精細的參數起始值，則實際上保證了收斂。
總之，基于廣義線性混合效應模型的 EWLM Meta 分析技術，包括二項式-正態模型和超幾何正態模型，可方便地實現稀有二分類數據的隨機效應 Meta 分析。
附錄：SAS 代碼
SAS 代碼 1：資料錄入
data Example_data;
input study td ta ttl cd ca ctl;
datalines;
1 5 20 25 6 19 25
2 1 8 9 2 14 16
……
/*連續性校正*/
data Example_data_CR;
set Example_data;
if (td=0) then do;
td = td+0.5;
ttl = ttl+0.5;
end;
if (cd=0) then do;
cd = cd+0.5;
ctl = ctl+0.5;
end;
run;
SAS 代碼 2：單臂研究中基于 PROC MIXED 程序的 NNM 模型（包含固定效應模型和隨機效應模型）
/*計算效應值及其標準誤（est）*/
data Example_data_CR_study;
set Example_data_CR;
logodds=log(td/(ttl-td));
est=1/td+1/(ttl-td);
run;
/*比例的固定效應模型*/
Proc mixed method =ml data=Example_data_CR_study;
class study;
model logodds =/ s cl;
repeated /group=study;
parms / parmsdata=Example_data_CR_study
eqcons=1 to 22;
ods output CovParms=CovP;
run;
/*隨機效應模型*/
/*根據固定效應模型的協方差參數，設置隨機效應模型的協方差*/
data covp;
set covp;
drop _LABEL_;
drop CovParm Group ;
rename Estimate=est;
run;
data start;
input est;
cards;
0.01
run;
data CovP;
set start CovP;
run;
proc print; run;
Proc mixed cl method=ml data=Example_data_CR_study;
class study;
model logodds= / s cl;
random int/ subject=study s;
repeated /group=study;
parms / parmsdata = CovP
eqcons = 2 to 23;
run;
SAS 代碼 3：單臂研究中基于 PROC NLMIXED 程序的 NNM 模型
proc nlmixed data=Example_data_CR qpoints=100;
parms theta= -2.5 thetai= -2.5 sigmasq=0.5;
logodds=log(td/(ttl-td));
sesq=1/td+1/(ttl-td);
model logodds ~ normal(thetai,sesq);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 4：單臂研究中基于 PROC NLMIXED 程序的 BNM 模型
proc nlmixed data=Example_data qpoints=100;
parms theta= -2.5 to -0.5 by 0.5 sigmasq=0.5;
pi=exp(logodds)/(1+exp(logodds));
model td ~ binomial(ttl,pi);
random logodds ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 5：雙臂研究中基于 PROC NLMIXED 程序的 NNM 模型
proc nlmixed data=Example_data_CR qpoints=100;
parms theta= -0.4 sigmasq=0.2;
logodds=log((td/(ttl-td))/(cd/(ctl-cd)));
sesq=1/td+1/(ttl-td)+1/cd+1/(ctl-cd);
model logodds ~ normal(thetai,sesq);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 6：雙臂研究中基于 PROC NLMIXED 程序的 BNM 模型
proc nlmixed data=Example_data qpoints=100;
parms theta= -2.5 to -0.5 by 0.5 sigmasq=0.5;
yi=td+cd;
pi=exp(log(ttl/ctl)+thetai)/(1+exp(log(ttl/ctl)+thetai));
model td ~ binomial(yi,pi);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
SAS 代碼 7：雙臂研究中基于 PROC NLMIXED 程序的 HNM 模型
proc nlmixed data=Example_data qpoints=100;
parms theta=-0.4 sigmasq=0 to 0.5 by 0.05;
y=td+cd;
num = lgamma(ttl+1) - lgamma(td+1) - lgamma(ttl-td+1)
+ lgamma(ctl+1) - lgamma(cd+1) - lgamma(ctl-cd+1)
+ td*thetai;
den = 0;
lo=max(0,y-ctl); hi=min(y,ttl);
do j=lo to hi;
z = lgamma(ttl+1) - lgamma(j+1) - lgamma(ttl-j+1)
+ lgamma(ctl+1) - lgamma(y-j+1) - lgamma(ctl-y+j+1)
+ j*thetai;
den = den + exp(z);
end;
ll = num - log(den);
model td ~ general(ll);
random thetai ~ normal(theta, sigmasq) subject=study;
run;
