混合效應模型（mixed effect model，MIXED）是隨機效應模型和固定效應模型的推廣，即同時包含固定效應和隨機效應^[1]。多水平 Meta 分析模型是混合效應模型的應用推廣，將傳統模型中單一的隨機誤差項分解到與相對應的研究水平上，其擬合效應不僅能使 Meta 分析結果更為穩健與合理，而且能通過對協變量的解釋指導臨床具體問題。SAS PROC MIXED 是 SAS 的內置程序，是目前進行線性混合效應模型分析的最佳程序之一^[2]。非線性混合效應模型（non-linear mixed effect model，NLMIXED），也被稱為非線性隨機效應模型、多水平非線性模型或非線性分層模型。在 Meta 分析領域，MIXED 和 NLMIXED 不僅能識別和估計研究間和研究內的變異，而且也考慮了解釋變量（協變量）與反應變量參數間的統計學關系。不同的是，MIXED 允許固定效應和隨機效應進入模型的線性部分，而 NLMIXED 則進入非線性部分。此外，MIXED 要求反應變量服從正態分布，而 NLMIXED 適用的反應變量的分布更廣，可以服從正態分布、二項分布或泊松分布等。因此，NLMIXED 的應用范圍更廣。SAS 軟件從 8.0 版本引入 PROC NLMIXED 過程中，專門用來進行非線性混合效應模型的擬合與評價。近年來，基于 PROC NLMIXED 的 Meta 分析方法被越來越多系統評價員采用，帶來了廣闊的應用前景^[3]。本研究目的是通過實例分析，比較利用 SAS MIXED 和 SAS NLMIXED 實現線性或非線性多水平模型的 Meta 分析的差異，并提供相關的 SAS 代碼。

1   Meta 分析模型簡介

1.1   D-L 隨機效應模型理論及介紹

DerSimonian 和 Laird 提出的逆方差隨機效應模型是目前最受推薦和廣泛使用的 Meta 方法，該模型被認為是經典的隨機效應 Meta 分析模型^{[4, 5]}。它通過計算來自多個獨立研究的效應值的加權平均，最終獲得合并效應值，從而評價干預效果。其模型可以表示如下：

1.2   多水平 Meta 分析模型理論及介紹

盡管 D-L 模型深受歡迎，但這種分析技術仍然具有爭議。一個重要的原因是這種模型無法評價研究結果的異質性。系統評價員需要通過其他統計方法來理解和解釋異質性。多水平模型可以納入包含研究水平的協變量，從而相對容易地實現異質性分析^[6]。在觀察性研究中，協變量（covariates）可被用于調整非等效或非隨機對照試驗的組間差異，實現統計學上的校正，以保證組間可比性。在 Meta 分析領域，協變量同樣可用于調整 Meta 分析中的非等效研究結果，并獲得評估治療效果的實質性信息^[7]。根據模型理論的不同，多水平 Meta 分析模型可分為線性混合效應多水平模型和非線性混合效應多水平模型。

混合效應線性模型是一般線性回歸模型的擴展，其公式表示如下^{[8, 9]}。

對于納入分析的 k 個研究中，Y_i是i=1...k個研究的觀測效應值。X_i表示第 i 個研究中影響固定效應值的協變量（covariates）。β 是固定效應參數的向量集。Z_i表示影響隨機效應值的協變量。b 是隨機效應參數的向量，或者是基于研究間水平（between-study level）的殘差項，而 e 是反映研究內水平殘差的矩陣。

在 Meta 分析領域，非線性混合效應多水平模型建立的統計學模型就是確定研究間變異和研究內變異的過程。對于納入分析的 k 個研究（i=1...k）中，通常有 p 個待分析變量（j=1...p），從而得到干預資料（X_ij，Y_ij）。欲分析Y_ij與X_ij之間的非線性關系，可建立兩個水平正態誤差非線性混合效應模型：

水平 1：研究內模型

式中Y_ij：第 i 個研究中干預 j 的觀測值；：第 i 個研究內r×1維回歸參數的向量集；X_ij：研究內p_i×t維的回歸設計矩陣，代表 t 個獨立的待分析研究內變量。代表維的X_ij和之間的非線性函數。：代表正態隨機誤差向量。

水平 2：研究間模型

式中A_i：是r×s已知設計矩陣。已知設計矩陣。β：表示未知的s×1維固定效應總體參數向量。b_i：未知的獨立同分布正態隨機效應參數向量。b_i服從。反映的是研究間方差-協方差矩陣，通常設置為無結構協方差矩陣，也可以根據研究目的，選擇自回歸矩陣或復合對稱矩陣等。公式中，b_i和是相互獨立的。

為便于書寫及理解，可將公式（3）寫出向量模式：

式中

因此，服從，其中矩陣的維數取決于 i。

對于 k 個研究，則滿足：

和 。

此時，完整的非線性混合效應模型可以寫成：

式中，

與線性混合效應模型一樣，非線性混合效應模型假設研究內的協方差結構為常數（即每個研究內的和γ是常數）。當研究內變異存在異方差性時，可以通過選擇條件獨立矩陣、自回歸矩陣、復合對稱矩陣等來解決。在 SAS 軟件中，混合效應線性模型是通過 PROC MIXED 過程步擬合，而非線性混合效應模型通過 PROC NLMIXED 過程步擬合。

2   資料與方法

2.1   示例介紹

為便于說明上述模型在 SAS 軟件中的具體操作方法，同時比較 3 種模型的 Meta 分析結果是否存在差異，我們結合實例數據進行模型擬合。本文的示例數據來源于 Shim 等^[10]發表的輔助手術降低宮頸癌根治性同步放化療患者局部復發風險的系統評價。該系統評價旨在探討輔助子宮切除術對同期放化療局部晚期宮頸癌患者預后的影響，共納入 8 個研究。本文只選取其中局部復發和淋巴結轉移比例數據作為示例。收集到的數據如下：“studyName”代表納入研究；“n11、n12、tot1、LNM1”分別代表治療組的復發例數、無復發例數、組樣本量、淋巴結轉移比例；“n21、n22、tot2、LNM2”分別代表對照組的復發例數、無復發例數、組樣本量、淋巴結轉移比例。假設我們關心的結局如下：① 不同治療模式（干預）對局部復發率的影響，主要的效應指標是對數比值比（Ln_OR）；② 淋巴結轉移比例作為協變量對局部復發的影響。該數據集的詳細信息如表 1 所示。

將表 1 的數據錄入 SAS 軟件，并存儲在名為 Sample_data_original 的數據集中。在單變量模型背景下，我們引入一個研究內水平的協變量 covariate_LNM，代表研究水平的淋巴結轉移率，其計算公式為（研究組淋巴結轉移例數+對照組轉移例數）/總樣本量；同時對 covariate_LNM 進行中心化，獲得協變量 covariate_LNM_center。

2.2   SAS 代碼介紹

數據集錄入及對連續變量 covariate_LNM 進行中心化處理的代碼見附錄。使用 PROC MIXED、PROC NLMIXED 程序，分別擬合線性或非線性多水平模型的 Meta 分析。同時給出基于 D-L 隨機效應模型的 Meta 分析結果，作為擬合的對照標準。受篇幅限制，計算 D-L 隨機效應模型的 SAS 代碼，可通過聯系本文通訊作者獲取。

3   結果

3.1   沒有協變量的不同模型 Meta 分析結果

在沒有協變量的情況下，基于雙變量隨機效應模型的 PROC MIXED 和非線性混合效應模型的 PROC NLMIXED 的 OR 合并效應值分別為［0.63，95%CI（0.46，0.87），P=0.005 7］和［0.60，95%CI（0.39，0.81），P=0.000 3］。非線性混合效應模型的擬合結果最接近 D-L 模型；而基于 PROC MIXED 過程步的單變量隨機效應模型給出的 OR 估計值與 D-L 模型相似，但標準誤偏大（表 2）。

3.2   引入協變量的不同模型 Meta 分析結果

在帶有協變量情況下，單變量隨機效應模型給出的標準誤最大，雙變量隨機效應模型和非線性混合效應模型 OR 效應值為［0.65，95%CI（0.47，0.91），P=0.011］和［0.59，95%CI（0.38，0.80），P=0.000 3］。協變量 OR 效應值為［2.70，95%CI（0.16，45.23），P>0.05］和［1.86，95%CI（?0.07，3.79），P=0.06］（表 3）。

4   討論

Meta 分析方法在醫學研究中越來越受重視，因具有相似治療方案的臨床研究可通過 Meta 分析獲得有關治療效果的合并效應估計從而起到擴大樣本量的作用。由 DerSimonian 和 Laird 提出的隨機效應模型（通常稱為 DerSimonian-Laird 模型，或 D-L 模型）是最常用的 Meta 分析模型^[4]。但該模型無法實現協變量分析，因此無法評價研究結果的異質性大小。多水平 Meta 分析模型可納入具有研究水平的協變量，從而用于調整 Meta 分析研究中的協變量對總體治療效應的影響信息。在許多統計學領域中，變量之間的關系是比較復雜的。最常見的處理方法是假設變量之間存在簡單的線性關系，從而使統計過程變得簡單。目前基于線性假設的多水平效應模型是多水平 Meta 分析最常用的模型^[2]。然而在實際情況中，這種線性假設可能并不適合，故變量間的非線性關系可能更具有一般性^[11]。線性混合效應模型并不適用于具有非線性關系的數據。非線性關系會導致數據建模存在一定的局限性。近年來，在線性混合效應模型基礎上，統計學家對線性模型進行擴展，根據變量間的非線性關系，建立不同的非線性混合效應模型，并根據因變量的分布特征擬合混合分布模型，即非線性混合效應模型^[11]。在 SAS 軟件中，分別采用 PROC MIXED 和 PROC NLMIXED 擬合線性混合效應模型和非線性混合效應模型。一般來說，PROC MIXED 適用于連續型變量，并且假設數據服從正態分布；而 PROC NLMIXED 適用于離散型變量，如多水平 logistic 回歸模型、多水平泊松模型等。而在實際應用中，由于 PROC NLMIXED 同樣可擬合正態模型數據，因此應用范圍更廣^[12]。多水平 Meta 分析模型可在 STATA、R 語言等常用統計軟件上實現，但 SAS 軟件更方便。SAS 可直接調用 SAS MIXED 和 NLMIXED 程序實現統計分析，而且擁有更強大的編程能力，因此是實現多水平 Meta 分析模型的優選軟件。

PROC NLMIXED 允許數據的參數服從一種條件分布（隨機效應），這種條件既可以是標準分布（正態分布、二項式分布、泊松分布等），也可以是 SAS 使用者自定義開發的數據分布。PROC NLMIXED 通過最大化隨機效應的聯合似然值來擬合非線性混合模型，并采用不同的積分逼近方法獲得聯合似然值，主要包括自適應高斯積分和一階泰勒級數積分。PROC NLMIXED 可通過多種技術來計算最大似然值，默認的是對偶擬牛頓算法（dual quasi-Newton algorithm）。使用 PROC NLMIXED 方法擬合的模型可以看做是使用 PROC MIXED 過程擬合的隨機參數模型的一般化。這就允許隨機參數使得模型非線性化，而在 MIXED 過程中是線性的。在 PROC MIXED 中最大似然值的計算方法可以是最大似然法（maximum likelihood，ML）和受約束最大似然法（restricted maximum likelihood，REML），而在 PRCO NLMIXED 中只能采用 ML，這是因為 NLMIXED 中 REML 的模擬包括所有固定參數的高維積分，這種積分在解析形式中一般不可用^[13]。在 PROC NLMIXED 過程步出現之前，很難通過多水平隨機效應模型擬合二項式數據。鑒于 PROC NLMIXED 功能的強大，越來越多系統評價員采用該模型擬合 Meta 分析，但遺憾的是，目前沒有學者比較二種統計模型的實施效果差別。本文的目的是比較基于多水平隨機效應模型的 PROC MIXED 和 PROC NLMIXED 擬合二項式結局 Meta 分析的差別。在 Meta 分析研究中，比值比（odds ratio，OR）是最常用的統計學指標，雖然 OR 經常不服從正態分布，但其對數轉換值 ln(OR) 近似服從正態分布^[14]。盡管本文以 OR 值作為示例，但本文提出的擬合模型同樣適用于其他效應指標，如相對危險度（relative risk，RR）和比例之間的差異。

本文介紹了沒有協變量的 Meta 分析模型和帶有協變量的 Meta 分析模型。在沒有協變量的情況下，比較了目前常用的 4 種 Meta 分析模型的結果差異。由于 D-L 模型是最常用的模型，因此 D-L 模型的計算結果被作為對照。D-L 模型計算可通過多種統計軟件實現，包括 RevMan 軟件、STATA、R 軟件等。在沒有協變量的情況下，非線性混合效應模型的擬合結果最接近 D-L 模型；而基于 PROC MIXED 過程步的單變量隨機效應模型給出的 OR 估計值與 D-L 模型相似，但標準誤偏大，這是因為 PROC MIXED 采用的是最大似然法。

本文同時給出了雙變量隨機效應模型的 Meta 分析結果。雙變量隨機效應模型將研究組干預和對照組干預分別擬合成 exp 和 con 兩個變量，經 PROC MIXED 擬合后，可給出 Meta 分析的效應值 OR，并給出相應的協方差結構^[15]。因此，雙變量隨機效應模型自然地解釋了研究中治療（或暴露）與對照組之間的潛在相關性。但需要注意的是，采用雙變量隨機效應模型擬合的協變量效應值 OR 具有很大的標準誤，說明使用該模型擬合協變量并不準確，在實踐中應少用。在 4 種模型中，單變量隨機效應模型給出的標準誤最大，而且 P 值效應最小，相應的可信區間也更寬，因此，我們認為該模型擬合效率較差，不推薦使用。

在考慮協變量的情況下，筆者給出了三種模型的統計結果。正如前述，單變量隨機效應模型給出的標準誤最大，從而給出了很寬的置信區間，結果 P 值大于0.05，與其他模型統計結果相反。在非線性混合效應模型構建中，采用對數正態模型（Logit-Normal Model）作為鏈接函數^[16]。這兩種模型均能給出參數估計的協方差矩陣，以非線性混合效應模型為例，給出了 PROC NLMIXED 的協方差結果。基于兩種模型的分析結果，表明輔助手術可以降低宮頸癌根治性放化療患者的局部復發風險，而淋巴結是否轉移對局部復發未見顯著影響。

需要注意的是：① 不同的方差-協方差結構，反映了數據之間不同的相關性模式。使用雙變量隨機效應模型擬合線性混合效應模型時，有時需要指定或修改協方差的結構形式，以本文示例數據為例，如果將協方差結構指定為“不規則結構（UN）”則模型不收斂，而指定為“因子分析結構 fa0(2)”時，模型得到了很好地收斂。收斂性問題是 PROC MIXED 經常遇見的問題，需要使用者根據實際調整。對多種協方差結構均能收斂時，使用者可以通過 Akaike 信息準則（Akaike information criterion，AIC）、Schwarz 信息準則（Schwarz information criterion，BIC）、-2 倍的對數似然函數（-2 res log likelihood）來確定。AIC 和 BIC 取值越小，說明對應的方差和協方差結構對數據的擬合效果越好。② 當事件數較少時，基于線性模型的 PROC MIXED 經常會失效，此時非線性混合效應模型表現出對稀疏數據的可靠處理，在事件數較少時，PROC NLMIXED 仍可以表現出強大而靈活的建模能力。③ PROC NLMIXED 具有強大的編程能力，可以根據研究者的需要進行代碼修改，甚至自定義數據的分布，這是 PROC NLMIXED 最大的優點。④ 在本文示例數據中，PROC NLMIXED 將協變量信息處理為研究間水平，即協變量以研究為單位，而 PROC MIXED 將協變量信息處理為研究內水平，即協變量以組別為單位。兩者處理結果稍有不同。

綜上所述，鑒于 PROC NLMIXED 具有強大的編程能力以及非線性混合效應模型對稀疏數據的靈活的建模能力，PROC NLMIXED 在 Meta 分析領域將發揮越來越重要的作用。

附錄代碼：

/*數據集建立及資料錄入*/

title "數據集建立及資料錄入";

data sample_data_original;

　input studyName$ n11 tot1 n21 tot2 LNM1 LNM2;

　　　n12=tot1-n11;n22=tot2-n21;

　　　or = (n11*n22)/(n21*n12);Ln_OR = log(or);

　　　v_Ln_OR = 1/n11 + 1/n12 + 1/n21 + 1/n22;

　　　w_Ln_OR = 1/v_Ln_OR;

　　　row = _n_;col =_n_;value = v_Ln_OR;study = _n_;

　　　covariate_LNM=round(LNM1/LNM2,0.00001);

　datalines;

　　1　11　119　21　121　0.06　0.04

……

run;

/*對數值變量 covariate_LNM 進行中心化*/

proc sql noprint;

select mean(covariate_LNM) into: mean_covariate_LNM from sample_data_original;

quit;

　%put &mean_covariate_LNM.;

data sample_data_original;

　set sample_data_original;

　mean_covariate_LNM=&mean_covariate_LNM.;

　covariate_LNM_center=mean_covariate_LNM-covariate_LNM;

run;

proc print data=sample_data_original;run;

title "模型 1：基于線性模型的沒有協變量的隨機效應模型";

proc mixed data=sample_data_original cl;

　class study;

　model Ln_OR = / ddfm=sat s cl;

　random study / gdata=sample_data_original;

　ods output SolutionF=solF_Liner_REM;

run;

Data OR_solF_Liner_REM;

　set solF_Liner_REM;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_solF_Liner_REM;run;

title "模型 2：基于線性模型的帶有協變量的隨機效應模型";

proc mixed data=sample_data_original cl;

　class study;

　model Ln_OR =covariate_LNM_center / ddfm=sat s cl;

　random study / gdata=sample_data_original;

　ods output SolutionF=solF_Liner_REM_co;

run;

Data OR_solF_Liner_REM_co;

　set solF_Liner_REM_co;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_solF_Liner_REM_co;run;

/*構造二項式觀測數據集*/

data sample_data_binomial;

　set sample_data_original;

　　event = n11;noevent = n12;

　　tot = tot1;covariate_LNM=LNM1;

　　　ln_OR=log(n11/n12);est=1/n11+1/n12;

　　treat = 1;

　　output;

　　event = n21;tot = tot2;

　　noevent = n22;covariate_LNM=LNM2;

　　　ln_OR=log(n21/n22);est=1/n21+1/n22;

　　treat = 0;

　　output;

　　keep study treat event noevent tot covariate_LNM ln_OR est;

run;

/*對數值變量 covariate_LNM 進行中心化*/

proc sql noprint;

　　select mean(covariate_LNM) into: mean_covariate_LNM from sample_data_binomial;

quit;

　　%put &mean_covariate_LNM.;

data sample_data_binomial;

　　set sample_data_binomial;

　　　mean_covariate_LNM=&mean_covariate_LNM.;

　　covariate_LNM_center=covariate_LNM-mean_covariate_LNM;

　　　arm=_n_;exp=(treat=1); con=(treat=0);

run;

proc print data=sample_data_binomial;run;

data covvarsa;

input est;

　　datalines;

　　0.0002

　　0.0001

　　0.0003

;

data covvars2;

set covvarsa sample_data_binomial;

run;

Proc print data= covvars2;run;

title "模型 3：基于線性模型的沒有協變量的雙變量模型，采用 fa0(2) 協方差結構";

Proc mixed cl method=ml data=sample_data_binomial asycov;

　class study arm;

　model ln_OR= exp con / noint s cl covb ddf=1000,1000;

　random exp con/ subject=study type=fa0(2) s;

　repeated /group=arm;

　estimate 'difference' exp 1 con -1/cl df=1000;

　parms /parmsdata=covvars2 eqcons=4 to 19;

　ods output CovParms=CovParms_Bivariate;

　ods output SolutionF=solF_Bivariate;

　ods output Estimates=Estimates_Bivariate;

run;

Data OR_Estimates_Bivariate;

　set Estimates_Bivariate;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_Estimates_Bivariate;run;

title "模型 4：基于線性模型的含協變量的雙變量模型，采用 fa0(2) 協方差結構";

Proc mixed cl method=ml data=sample_data_binomial asycov;

　class study arm;

　model ln_OR= exp con covariate_LNM_center/ noint s cl covb ddf=1000, 1000;

　random exp con/ subject=study type=fa0(2) s;

　repeated /group=arm;

　estimate 'difference' exp 1 con -1/cl df=1000;

　parms /parmsdata=covvars2 eqcons=4 to 19;

　ods output CovParms=CovParms_Bivariate;

　ods output SolutionF=solF_Bivariate;

　ods output Estimates=Estimates_Bivariate;

run;

Data OR_Estimates_Bivariate;

　set Estimates_Bivariate;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_Estimates_Bivariate;run;

data sample_data_bi_NLMIXED;

　set sample_data_original;

　event = n11;tot =tot1;treat = 1;

　a_fem = covariate_LNM_center;

　b_fem = covariate_LNM_center;

　output;

　event = n21;tot = tot2;treat = 0;

　a_fem = 0;

　b_fem = covariate_LNM_center;;

　output;

　keep study treat event tot a_fem b_fem;

run;

title "模型 5：基于非線性模型的無協變量的對數正態模型";

proc nlmixed data=sample_data_bi_NLMIXED cov corr;

　parms beta0=1 beta1=-2 vc=1;

　bounds vc >= 0;

　eta = beta0 + beta1*treat + u;

　expeta = exp(eta);

　p = expeta/(1+expeta);

　model event ～ binomial(tot,p);

　random u ～ normal(0,vc) subject=study;

　estimate 'OR' exp(beta1);

run;

title "模型 6：基于非線性模型的具有協變量的對數正態模型";

proc nlmixed data=sample_data_bi_NLMIXED cov corr;

　parms beta0=1 beta1=-2 beta2=1 vc=1;

　bounds vc >= 0;

　eta = beta0 + beta1*treat + beta2*a_fem + u;

　expeta = exp(eta);

　p = expeta/(1+expeta);

　model event ～ binomial(tot,p);

　random u ～ normal(0,vc) subject=study;

　estimate 'OR' exp(beta1);

　estimate 'OR Prop LNM' exp(beta2);

run;

混合效應模型（mixed effect model，MIXED）是隨機效應模型和固定效應模型的推廣，即同時包含固定效應和隨機效應^[1]。多水平 Meta 分析模型是混合效應模型的應用推廣，將傳統模型中單一的隨機誤差項分解到與相對應的研究水平上，其擬合效應不僅能使 Meta 分析結果更為穩健與合理，而且能通過對協變量的解釋指導臨床具體問題。SAS PROC MIXED 是 SAS 的內置程序，是目前進行線性混合效應模型分析的最佳程序之一^[2]。非線性混合效應模型（non-linear mixed effect model，NLMIXED），也被稱為非線性隨機效應模型、多水平非線性模型或非線性分層模型。在 Meta 分析領域，MIXED 和 NLMIXED 不僅能識別和估計研究間和研究內的變異，而且也考慮了解釋變量（協變量）與反應變量參數間的統計學關系。不同的是，MIXED 允許固定效應和隨機效應進入模型的線性部分，而 NLMIXED 則進入非線性部分。此外，MIXED 要求反應變量服從正態分布，而 NLMIXED 適用的反應變量的分布更廣，可以服從正態分布、二項分布或泊松分布等。因此，NLMIXED 的應用范圍更廣。SAS 軟件從 8.0 版本引入 PROC NLMIXED 過程中，專門用來進行非線性混合效應模型的擬合與評價。近年來，基于 PROC NLMIXED 的 Meta 分析方法被越來越多系統評價員采用，帶來了廣闊的應用前景^[3]。本研究目的是通過實例分析，比較利用 SAS MIXED 和 SAS NLMIXED 實現線性或非線性多水平模型的 Meta 分析的差異，并提供相關的 SAS 代碼。

1   Meta 分析模型簡介

1.1   D-L 隨機效應模型理論及介紹

DerSimonian 和 Laird 提出的逆方差隨機效應模型是目前最受推薦和廣泛使用的 Meta 方法，該模型被認為是經典的隨機效應 Meta 分析模型^{[4, 5]}。它通過計算來自多個獨立研究的效應值的加權平均，最終獲得合并效應值，從而評價干預效果。其模型可以表示如下：

1.2   多水平 Meta 分析模型理論及介紹

盡管 D-L 模型深受歡迎，但這種分析技術仍然具有爭議。一個重要的原因是這種模型無法評價研究結果的異質性。系統評價員需要通過其他統計方法來理解和解釋異質性。多水平模型可以納入包含研究水平的協變量，從而相對容易地實現異質性分析^[6]。在觀察性研究中，協變量（covariates）可被用于調整非等效或非隨機對照試驗的組間差異，實現統計學上的校正，以保證組間可比性。在 Meta 分析領域，協變量同樣可用于調整 Meta 分析中的非等效研究結果，并獲得評估治療效果的實質性信息^[7]。根據模型理論的不同，多水平 Meta 分析模型可分為線性混合效應多水平模型和非線性混合效應多水平模型。

混合效應線性模型是一般線性回歸模型的擴展，其公式表示如下^{[8, 9]}。

對于納入分析的 k 個研究中，Y_i是i=1...k個研究的觀測效應值。X_i表示第 i 個研究中影響固定效應值的協變量（covariates）。β 是固定效應參數的向量集。Z_i表示影響隨機效應值的協變量。b 是隨機效應參數的向量，或者是基于研究間水平（between-study level）的殘差項，而 e 是反映研究內水平殘差的矩陣。

在 Meta 分析領域，非線性混合效應多水平模型建立的統計學模型就是確定研究間變異和研究內變異的過程。對于納入分析的 k 個研究（i=1...k）中，通常有 p 個待分析變量（j=1...p），從而得到干預資料（X_ij，Y_ij）。欲分析Y_ij與X_ij之間的非線性關系，可建立兩個水平正態誤差非線性混合效應模型：

水平 1：研究內模型

式中Y_ij：第 i 個研究中干預 j 的觀測值；：第 i 個研究內r×1維回歸參數的向量集；X_ij：研究內p_i×t維的回歸設計矩陣，代表 t 個獨立的待分析研究內變量。代表維的X_ij和之間的非線性函數。：代表正態隨機誤差向量。

水平 2：研究間模型

式中A_i：是r×s已知設計矩陣。已知設計矩陣。β：表示未知的s×1維固定效應總體參數向量。b_i：未知的獨立同分布正態隨機效應參數向量。b_i服從。反映的是研究間方差-協方差矩陣，通常設置為無結構協方差矩陣，也可以根據研究目的，選擇自回歸矩陣或復合對稱矩陣等。公式中，b_i和是相互獨立的。

為便于書寫及理解，可將公式（3）寫出向量模式：

式中

因此，服從，其中矩陣的維數取決于 i。

對于 k 個研究，則滿足：

和 。

此時，完整的非線性混合效應模型可以寫成：

式中，

與線性混合效應模型一樣，非線性混合效應模型假設研究內的協方差結構為常數（即每個研究內的和γ是常數）。當研究內變異存在異方差性時，可以通過選擇條件獨立矩陣、自回歸矩陣、復合對稱矩陣等來解決。在 SAS 軟件中，混合效應線性模型是通過 PROC MIXED 過程步擬合，而非線性混合效應模型通過 PROC NLMIXED 過程步擬合。

2   資料與方法

2.1   示例介紹

為便于說明上述模型在 SAS 軟件中的具體操作方法，同時比較 3 種模型的 Meta 分析結果是否存在差異，我們結合實例數據進行模型擬合。本文的示例數據來源于 Shim 等^[10]發表的輔助手術降低宮頸癌根治性同步放化療患者局部復發風險的系統評價。該系統評價旨在探討輔助子宮切除術對同期放化療局部晚期宮頸癌患者預后的影響，共納入 8 個研究。本文只選取其中局部復發和淋巴結轉移比例數據作為示例。收集到的數據如下：“studyName”代表納入研究；“n11、n12、tot1、LNM1”分別代表治療組的復發例數、無復發例數、組樣本量、淋巴結轉移比例；“n21、n22、tot2、LNM2”分別代表對照組的復發例數、無復發例數、組樣本量、淋巴結轉移比例。假設我們關心的結局如下：① 不同治療模式（干預）對局部復發率的影響，主要的效應指標是對數比值比（Ln_OR）；② 淋巴結轉移比例作為協變量對局部復發的影響。該數據集的詳細信息如表 1 所示。

將表 1 的數據錄入 SAS 軟件，并存儲在名為 Sample_data_original 的數據集中。在單變量模型背景下，我們引入一個研究內水平的協變量 covariate_LNM，代表研究水平的淋巴結轉移率，其計算公式為（研究組淋巴結轉移例數+對照組轉移例數）/總樣本量；同時對 covariate_LNM 進行中心化，獲得協變量 covariate_LNM_center。

2.2   SAS 代碼介紹

數據集錄入及對連續變量 covariate_LNM 進行中心化處理的代碼見附錄。使用 PROC MIXED、PROC NLMIXED 程序，分別擬合線性或非線性多水平模型的 Meta 分析。同時給出基于 D-L 隨機效應模型的 Meta 分析結果，作為擬合的對照標準。受篇幅限制，計算 D-L 隨機效應模型的 SAS 代碼，可通過聯系本文通訊作者獲取。

3   結果

3.1   沒有協變量的不同模型 Meta 分析結果

在沒有協變量的情況下，基于雙變量隨機效應模型的 PROC MIXED 和非線性混合效應模型的 PROC NLMIXED 的 OR 合并效應值分別為［0.63，95%CI（0.46，0.87），P=0.005 7］和［0.60，95%CI（0.39，0.81），P=0.000 3］。非線性混合效應模型的擬合結果最接近 D-L 模型；而基于 PROC MIXED 過程步的單變量隨機效應模型給出的 OR 估計值與 D-L 模型相似，但標準誤偏大（表 2）。

3.2   引入協變量的不同模型 Meta 分析結果

在帶有協變量情況下，單變量隨機效應模型給出的標準誤最大，雙變量隨機效應模型和非線性混合效應模型 OR 效應值為［0.65，95%CI（0.47，0.91），P=0.011］和［0.59，95%CI（0.38，0.80），P=0.000 3］。協變量 OR 效應值為［2.70，95%CI（0.16，45.23），P>0.05］和［1.86，95%CI（?0.07，3.79），P=0.06］（表 3）。

4   討論

Meta 分析方法在醫學研究中越來越受重視，因具有相似治療方案的臨床研究可通過 Meta 分析獲得有關治療效果的合并效應估計從而起到擴大樣本量的作用。由 DerSimonian 和 Laird 提出的隨機效應模型（通常稱為 DerSimonian-Laird 模型，或 D-L 模型）是最常用的 Meta 分析模型^[4]。但該模型無法實現協變量分析，因此無法評價研究結果的異質性大小。多水平 Meta 分析模型可納入具有研究水平的協變量，從而用于調整 Meta 分析研究中的協變量對總體治療效應的影響信息。在許多統計學領域中，變量之間的關系是比較復雜的。最常見的處理方法是假設變量之間存在簡單的線性關系，從而使統計過程變得簡單。目前基于線性假設的多水平效應模型是多水平 Meta 分析最常用的模型^[2]。然而在實際情況中，這種線性假設可能并不適合，故變量間的非線性關系可能更具有一般性^[11]。線性混合效應模型并不適用于具有非線性關系的數據。非線性關系會導致數據建模存在一定的局限性。近年來，在線性混合效應模型基礎上，統計學家對線性模型進行擴展，根據變量間的非線性關系，建立不同的非線性混合效應模型，并根據因變量的分布特征擬合混合分布模型，即非線性混合效應模型^[11]。在 SAS 軟件中，分別采用 PROC MIXED 和 PROC NLMIXED 擬合線性混合效應模型和非線性混合效應模型。一般來說，PROC MIXED 適用于連續型變量，并且假設數據服從正態分布；而 PROC NLMIXED 適用于離散型變量，如多水平 logistic 回歸模型、多水平泊松模型等。而在實際應用中，由于 PROC NLMIXED 同樣可擬合正態模型數據，因此應用范圍更廣^[12]。多水平 Meta 分析模型可在 STATA、R 語言等常用統計軟件上實現，但 SAS 軟件更方便。SAS 可直接調用 SAS MIXED 和 NLMIXED 程序實現統計分析，而且擁有更強大的編程能力，因此是實現多水平 Meta 分析模型的優選軟件。

PROC NLMIXED 允許數據的參數服從一種條件分布（隨機效應），這種條件既可以是標準分布（正態分布、二項式分布、泊松分布等），也可以是 SAS 使用者自定義開發的數據分布。PROC NLMIXED 通過最大化隨機效應的聯合似然值來擬合非線性混合模型，并采用不同的積分逼近方法獲得聯合似然值，主要包括自適應高斯積分和一階泰勒級數積分。PROC NLMIXED 可通過多種技術來計算最大似然值，默認的是對偶擬牛頓算法（dual quasi-Newton algorithm）。使用 PROC NLMIXED 方法擬合的模型可以看做是使用 PROC MIXED 過程擬合的隨機參數模型的一般化。這就允許隨機參數使得模型非線性化，而在 MIXED 過程中是線性的。在 PROC MIXED 中最大似然值的計算方法可以是最大似然法（maximum likelihood，ML）和受約束最大似然法（restricted maximum likelihood，REML），而在 PRCO NLMIXED 中只能采用 ML，這是因為 NLMIXED 中 REML 的模擬包括所有固定參數的高維積分，這種積分在解析形式中一般不可用^[13]。在 PROC NLMIXED 過程步出現之前，很難通過多水平隨機效應模型擬合二項式數據。鑒于 PROC NLMIXED 功能的強大，越來越多系統評價員采用該模型擬合 Meta 分析，但遺憾的是，目前沒有學者比較二種統計模型的實施效果差別。本文的目的是比較基于多水平隨機效應模型的 PROC MIXED 和 PROC NLMIXED 擬合二項式結局 Meta 分析的差別。在 Meta 分析研究中，比值比（odds ratio，OR）是最常用的統計學指標，雖然 OR 經常不服從正態分布，但其對數轉換值 ln(OR) 近似服從正態分布^[14]。盡管本文以 OR 值作為示例，但本文提出的擬合模型同樣適用于其他效應指標，如相對危險度（relative risk，RR）和比例之間的差異。

本文介紹了沒有協變量的 Meta 分析模型和帶有協變量的 Meta 分析模型。在沒有協變量的情況下，比較了目前常用的 4 種 Meta 分析模型的結果差異。由于 D-L 模型是最常用的模型，因此 D-L 模型的計算結果被作為對照。D-L 模型計算可通過多種統計軟件實現，包括 RevMan 軟件、STATA、R 軟件等。在沒有協變量的情況下，非線性混合效應模型的擬合結果最接近 D-L 模型；而基于 PROC MIXED 過程步的單變量隨機效應模型給出的 OR 估計值與 D-L 模型相似，但標準誤偏大，這是因為 PROC MIXED 采用的是最大似然法。

本文同時給出了雙變量隨機效應模型的 Meta 分析結果。雙變量隨機效應模型將研究組干預和對照組干預分別擬合成 exp 和 con 兩個變量，經 PROC MIXED 擬合后，可給出 Meta 分析的效應值 OR，并給出相應的協方差結構^[15]。因此，雙變量隨機效應模型自然地解釋了研究中治療（或暴露）與對照組之間的潛在相關性。但需要注意的是，采用雙變量隨機效應模型擬合的協變量效應值 OR 具有很大的標準誤，說明使用該模型擬合協變量并不準確，在實踐中應少用。在 4 種模型中，單變量隨機效應模型給出的標準誤最大，而且 P 值效應最小，相應的可信區間也更寬，因此，我們認為該模型擬合效率較差，不推薦使用。

在考慮協變量的情況下，筆者給出了三種模型的統計結果。正如前述，單變量隨機效應模型給出的標準誤最大，從而給出了很寬的置信區間，結果 P 值大于0.05，與其他模型統計結果相反。在非線性混合效應模型構建中，采用對數正態模型（Logit-Normal Model）作為鏈接函數^[16]。這兩種模型均能給出參數估計的協方差矩陣，以非線性混合效應模型為例，給出了 PROC NLMIXED 的協方差結果。基于兩種模型的分析結果，表明輔助手術可以降低宮頸癌根治性放化療患者的局部復發風險，而淋巴結是否轉移對局部復發未見顯著影響。

需要注意的是：① 不同的方差-協方差結構，反映了數據之間不同的相關性模式。使用雙變量隨機效應模型擬合線性混合效應模型時，有時需要指定或修改協方差的結構形式，以本文示例數據為例，如果將協方差結構指定為“不規則結構（UN）”則模型不收斂，而指定為“因子分析結構 fa0(2)”時，模型得到了很好地收斂。收斂性問題是 PROC MIXED 經常遇見的問題，需要使用者根據實際調整。對多種協方差結構均能收斂時，使用者可以通過 Akaike 信息準則（Akaike information criterion，AIC）、Schwarz 信息準則（Schwarz information criterion，BIC）、-2 倍的對數似然函數（-2 res log likelihood）來確定。AIC 和 BIC 取值越小，說明對應的方差和協方差結構對數據的擬合效果越好。② 當事件數較少時，基于線性模型的 PROC MIXED 經常會失效，此時非線性混合效應模型表現出對稀疏數據的可靠處理，在事件數較少時，PROC NLMIXED 仍可以表現出強大而靈活的建模能力。③ PROC NLMIXED 具有強大的編程能力，可以根據研究者的需要進行代碼修改，甚至自定義數據的分布，這是 PROC NLMIXED 最大的優點。④ 在本文示例數據中，PROC NLMIXED 將協變量信息處理為研究間水平，即協變量以研究為單位，而 PROC MIXED 將協變量信息處理為研究內水平，即協變量以組別為單位。兩者處理結果稍有不同。

綜上所述，鑒于 PROC NLMIXED 具有強大的編程能力以及非線性混合效應模型對稀疏數據的靈活的建模能力，PROC NLMIXED 在 Meta 分析領域將發揮越來越重要的作用。

附錄代碼：

/*數據集建立及資料錄入*/

title "數據集建立及資料錄入";

data sample_data_original;

　input studyName$ n11 tot1 n21 tot2 LNM1 LNM2;

　　　n12=tot1-n11;n22=tot2-n21;

　　　or = (n11*n22)/(n21*n12);Ln_OR = log(or);

　　　v_Ln_OR = 1/n11 + 1/n12 + 1/n21 + 1/n22;

　　　w_Ln_OR = 1/v_Ln_OR;

　　　row = _n_;col =_n_;value = v_Ln_OR;study = _n_;

　　　covariate_LNM=round(LNM1/LNM2,0.00001);

　datalines;

　　1　11　119　21　121　0.06　0.04

……

run;

/*對數值變量 covariate_LNM 進行中心化*/

proc sql noprint;

select mean(covariate_LNM) into: mean_covariate_LNM from sample_data_original;

quit;

　%put &mean_covariate_LNM.;

data sample_data_original;

　set sample_data_original;

　mean_covariate_LNM=&mean_covariate_LNM.;

　covariate_LNM_center=mean_covariate_LNM-covariate_LNM;

run;

proc print data=sample_data_original;run;

title "模型 1：基于線性模型的沒有協變量的隨機效應模型";

proc mixed data=sample_data_original cl;

　class study;

　model Ln_OR = / ddfm=sat s cl;

　random study / gdata=sample_data_original;

　ods output SolutionF=solF_Liner_REM;

run;

Data OR_solF_Liner_REM;

　set solF_Liner_REM;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_solF_Liner_REM;run;

title "模型 2：基于線性模型的帶有協變量的隨機效應模型";

proc mixed data=sample_data_original cl;

　class study;

　model Ln_OR =covariate_LNM_center / ddfm=sat s cl;

　random study / gdata=sample_data_original;

　ods output SolutionF=solF_Liner_REM_co;

run;

Data OR_solF_Liner_REM_co;

　set solF_Liner_REM_co;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_solF_Liner_REM_co;run;

/*構造二項式觀測數據集*/

data sample_data_binomial;

　set sample_data_original;

　　event = n11;noevent = n12;

　　tot = tot1;covariate_LNM=LNM1;

　　　ln_OR=log(n11/n12);est=1/n11+1/n12;

　　treat = 1;

　　output;

　　event = n21;tot = tot2;

　　noevent = n22;covariate_LNM=LNM2;

　　　ln_OR=log(n21/n22);est=1/n21+1/n22;

　　treat = 0;

　　output;

　　keep study treat event noevent tot covariate_LNM ln_OR est;

run;

/*對數值變量 covariate_LNM 進行中心化*/

proc sql noprint;

　　select mean(covariate_LNM) into: mean_covariate_LNM from sample_data_binomial;

quit;

　　%put &mean_covariate_LNM.;

data sample_data_binomial;

　　set sample_data_binomial;

　　　mean_covariate_LNM=&mean_covariate_LNM.;

　　covariate_LNM_center=covariate_LNM-mean_covariate_LNM;

　　　arm=_n_;exp=(treat=1); con=(treat=0);

run;

proc print data=sample_data_binomial;run;

data covvarsa;

input est;

　　datalines;

　　0.0002

　　0.0001

　　0.0003

;

data covvars2;

set covvarsa sample_data_binomial;

run;

Proc print data= covvars2;run;

title "模型 3：基于線性模型的沒有協變量的雙變量模型，采用 fa0(2) 協方差結構";

Proc mixed cl method=ml data=sample_data_binomial asycov;

　class study arm;

　model ln_OR= exp con / noint s cl covb ddf=1000,1000;

　random exp con/ subject=study type=fa0(2) s;

　repeated /group=arm;

　estimate 'difference' exp 1 con -1/cl df=1000;

　parms /parmsdata=covvars2 eqcons=4 to 19;

　ods output CovParms=CovParms_Bivariate;

　ods output SolutionF=solF_Bivariate;

　ods output Estimates=Estimates_Bivariate;

run;

Data OR_Estimates_Bivariate;

　set Estimates_Bivariate;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_Estimates_Bivariate;run;

title "模型 4：基于線性模型的含協變量的雙變量模型，采用 fa0(2) 協方差結構";

Proc mixed cl method=ml data=sample_data_binomial asycov;

　class study arm;

　model ln_OR= exp con covariate_LNM_center/ noint s cl covb ddf=1000, 1000;

　random exp con/ subject=study type=fa0(2) s;

　repeated /group=arm;

　estimate 'difference' exp 1 con -1/cl df=1000;

　parms /parmsdata=covvars2 eqcons=4 to 19;

　ods output CovParms=CovParms_Bivariate;

　ods output SolutionF=solF_Bivariate;

　ods output Estimates=Estimates_Bivariate;

run;

Data OR_Estimates_Bivariate;

　set Estimates_Bivariate;

　OR=exp(Estimate);

　OR_SE=(exp(Upper)-exp(Lower))/2*1.96;

　OR_Low=exp(Lower);

　OR_Up=exp(Upper);

　P_Value=Probt;

　keep OR OR_SE OR_Low OR_Up P_Value;

run;

proc print data=OR_Estimates_Bivariate;run;

data sample_data_bi_NLMIXED;

　set sample_data_original;

　event = n11;tot =tot1;treat = 1;

　a_fem = covariate_LNM_center;

　b_fem = covariate_LNM_center;

　output;

　event = n21;tot = tot2;treat = 0;

　a_fem = 0;

　b_fem = covariate_LNM_center;;

　output;

　keep study treat event tot a_fem b_fem;

run;

title "模型 5：基于非線性模型的無協變量的對數正態模型";

proc nlmixed data=sample_data_bi_NLMIXED cov corr;

　parms beta0=1 beta1=-2 vc=1;

　bounds vc >= 0;

　eta = beta0 + beta1*treat + u;

　expeta = exp(eta);

　p = expeta/(1+expeta);

　model event ～ binomial(tot,p);

　random u ～ normal(0,vc) subject=study;

　estimate 'OR' exp(beta1);

run;

title "模型 6：基于非線性模型的具有協變量的對數正態模型";

proc nlmixed data=sample_data_bi_NLMIXED cov corr;

　parms beta0=1 beta1=-2 beta2=1 vc=1;

　bounds vc >= 0;

　eta = beta0 + beta1*treat + beta2*a_fem + u;

　expeta = exp(eta);

　p = expeta/(1+expeta);

　model event ～ binomial(tot,p);

　random u ～ normal(0,vc) subject=study;

　estimate 'OR' exp(beta1);

　estimate 'OR Prop LNM' exp(beta2);

run;