混合效应模型（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;