鄭建清 黃碧芬 吳 敏 肖麗華

Meta分析被越來越廣泛地用于醫(yī)學(xué)研究，通過匯總來自多個(gè)獨(dú)立研究的資料來評估某種干預(yù)的有效性，也可以用于因果關(guān)系及率的綜合分析。當(dāng)研究人員所關(guān)心的結(jié)局變量是二分類變量時(shí)，這類數(shù)據(jù)可以被視為一個(gè)2×2的列聯(lián)表，每個(gè)單元格對應(yīng)于不同組中的事件數(shù)和無事件數(shù)(例如研究組和對照組中有疾病和無疾病的人數(shù))。然后，可以通過風(fēng)險(xiǎn)差(RD)、相對危險(xiǎn)度(RR)或者比值比(OR)等指標(biāo)來匯總這些研究中收集的信息。臨床中，另外一個(gè)常用的效應(yīng)指標(biāo)為發(fā)生率比率(IRR)，又稱為發(fā)病率密度，是基于每個(gè)研究組單獨(dú)記錄的事件計(jì)數(shù)隨時(shí)間的變化情況，例如每1 000人年。將這種類型的數(shù)據(jù)稱為發(fā)病率數(shù)據(jù)，也稱為發(fā)病率密度數(shù)據(jù)[1]。

目前已經(jīng)提出了兩種方法來進(jìn)行發(fā)病率數(shù)據(jù)的Meta分析，包括倒方差法和基于固定效應(yīng)模型的泊松回歸模型[2-4]。當(dāng)數(shù)據(jù)存在“0事件”(多個(gè)研究在研究的1或2個(gè)組中具有0事件)時(shí)，倒方差法是有問題的，因?yàn)楫?dāng)研究包含0事件時(shí)，研究特定的效應(yīng)值log IRR及其方差都無法定義。因此，最終分析的時(shí)候會省略含0事件研究。一種被廣泛推薦使用的補(bǔ)救措施是應(yīng)用連續(xù)性校正，但這種方法通常僅在只有少數(shù)研究包含0事件的情況下使用，否則會產(chǎn)生明顯的偏倚?；诠潭ㄐ?yīng)模型的泊松回歸模型可以解決暴露時(shí)間變化的情況，也可以用于解決數(shù)據(jù)中存在0事件的問題。另一種選擇是通過基于隨機(jī)效應(yīng)模型的泊松回歸模型來擴(kuò)展先前的模型，但這種方法很少被用于上述情況。其優(yōu)點(diǎn)在于，除了潛在地處理0事件和改變暴露時(shí)間之外，還可以解決干預(yù)效應(yīng)存在異質(zhì)性的問題[4]。

本研究的目的是介紹如何利用SAS軟件中的NLMIXED過程步實(shí)現(xiàn)基于廣義線性混合模型的泊松-正態(tài)模型(Poisson-Normal Model，PNM)和二項(xiàng)式-正態(tài)模型(Binomial-Normal Model，BNM)在稀有發(fā)病率數(shù)據(jù)中的Meta分析技術(shù)，尤其是包含0事件的數(shù)據(jù)。本文采用基于正態(tài)-正態(tài)模型(Normal -Normal Model，NNM)的方法(包括是否采用連續(xù)性校正)作為比較。本文介紹的基于廣義線性混合模型的PNM和BNM由Stijnen等人提出，并通過小型模擬試驗(yàn)證實(shí)是行之有效的方法[5]。本文提供了各種模型的SAS代碼，這些代碼簡潔明了，使用方便。

1 基于倒方差方法的正態(tài)-正態(tài)Meta分析模型簡介

其中假設(shè)標(biāo)準(zhǔn)誤Si是已知的，而均值θi是待估計(jì)的。由于異質(zhì)性的存在，真實(shí)值θi可能因研究而異，而θi也被假設(shè)服從正態(tài)分布，即：

θi?N(θ, σi2) (公式2)

參數(shù)θ是感興趣的主要參數(shù)，稱為整體效應(yīng)，而σ2描述了不同研究的真實(shí)效果的變異。公式1和2被稱為NNM。統(tǒng)計(jì)學(xué)推斷基于以下的似然性：

普通最大似然法(Maximum Likelihood，ML)或約束最大似然法(Restricted Maximum Likelihood，REML)均可用于估計(jì)參數(shù)。

當(dāng)原始效應(yīng)值是OR、風(fēng)險(xiǎn)比或相對風(fēng)險(xiǎn)度時(shí)，則應(yīng)該先將估計(jì)值進(jìn)行對數(shù)轉(zhuǎn)換，因?yàn)楹笳吒咏龖B(tài)分布的漸近理論。

1.2 倒方差方法理論及介紹 倒方差Meta方法通過對來自多個(gè)獨(dú)立研究的點(diǎn)估計(jì)的加權(quán)平均，計(jì)算干預(yù)效應(yīng)值的匯總估計(jì)來綜合來自多個(gè)來源的信息。這種方法目前被廣泛使用和推薦。對于NNM，倒方差的總體效應(yīng)估計(jì)可通過下列公式獲得

1.3 發(fā)病率密度數(shù)據(jù)的效應(yīng)參數(shù)及計(jì)算 發(fā)病率密度數(shù)據(jù)分為單臂研究和雙臂研究。對于前者，研究人員關(guān)心的效應(yīng)指標(biāo)是發(fā)病率(IR)，而后者是IRR。

1.3.1 單臂研究IR的Meta分析 對于單臂研究的IR密度數(shù)據(jù)背景下，每個(gè)研究中提供了事件Yi和總隨訪時(shí)間 Ti。每項(xiàng)研究的IR由Yi/Ti估算，而Meta分析的目的是估計(jì)總體IR。標(biāo)準(zhǔn)方法使用

如果某項(xiàng)研究的對照組或干預(yù)組中存在0事件時(shí)，通常是兩組分別增加0.5作“連續(xù)性校正”，盡管有學(xué)者提出可以使用其他數(shù)值進(jìn)行校正[6-8]。計(jì)算出單個(gè)研究的效應(yīng)值及其標(biāo)準(zhǔn)誤后，可以根據(jù)上述NNM模型進(jìn)行效應(yīng)值的合并。

1.3.2 雙臂研究IRR的Meta分析 對于雙臂研究，單項(xiàng)研究通常提供治療組和對照組。Meta分析的效應(yīng)值是IRR。令Yi0和Ti0為研究對照組中事件的數(shù)量和總的隨訪時(shí)間，類似地，治療組為Yi1和Ti1?，F(xiàn)在θi表示研究i的真實(shí)對數(shù)IRR，θ代表對應(yīng)的總體平均值。標(biāo)準(zhǔn)的Meta分析模型可以使用

對于含0事件研究，同樣需要上述校正方法進(jìn)行校正。

在單臂研究的發(fā)病率密度數(shù)據(jù)背景下，筆者主張用泊松分布代替正態(tài)分布：

Yi?Poisson(Tiexp(θi)) 其中θi?N(θ, σ2) (公式7)

PNM是具有偏移變量log(Ti)的隨機(jī)截距泊松回歸模型。

對于雙臂研究的發(fā)病率密度數(shù)據(jù)背景下，對于研究i，事件總數(shù)Yi= Yi0+ Yi1，而Yi1被假設(shè)服從二項(xiàng)分布：

其中θi?N(θ, σ2) (公式8)

BNM是一個(gè)隨機(jī)截距邏輯回歸模型，具有偏移量變量log(Ti1/Ti0)。

2 方法

2.1 待分析數(shù)據(jù) 來源于Schutz等發(fā)表的研究血管內(nèi)皮生長因子受體酪氨酸激酶抑制劑治療的癌癥患者發(fā)生致命不良事件(FAE)的風(fēng)險(xiǎn)[9]，“study”代表試驗(yàn)；“y1、n1、t1”分別代表治療組的事件數(shù)、組樣本量、中位生存時(shí)間；“y0、n0、t0”分別代表對照組的事件數(shù)、組樣本量、中位生存時(shí)間(表1)。假設(shè)關(guān)心的結(jié)局：①不同干預(yù)組單位人時(shí)的總體事件發(fā)生率以中位總生存為基準(zhǔn)，每1年總生存時(shí)間內(nèi)的總體事件發(fā)生率；②比較兩種干預(yù)對總體事件發(fā)生率的影響。

將表1的數(shù)據(jù)錄入SAS軟件，并存儲在名為Example_data的數(shù)據(jù)集中(因篇幅限制，數(shù)據(jù)步的SAS代碼簡要輸入)。在數(shù)據(jù)集中，存在2項(xiàng)單0試驗(yàn)，分別是study 3、study 10；3項(xiàng)雙0試驗(yàn)，分別是study 1、study 5、study 6，分別對含0試驗(yàn)進(jìn)行校正，校正后的數(shù)據(jù)保存在名為Example_data_CR的數(shù)據(jù)集。

表1 待分析數(shù)據(jù)集

2.2 SAS代碼介紹 本文提供的附錄代碼包括數(shù)據(jù)集的錄入，以及采用連續(xù)性校正的方法對數(shù)據(jù)集中的單0研究或雙0研究進(jìn)行校正，獲得數(shù)據(jù)集Example_data_CR。然后使用PROC NLMIXED程序，實(shí)現(xiàn)NNM模型、BNM模型、PNM模型的Meta分析。

3 結(jié)果

3.1 不同干預(yù)組總體事件發(fā)生率的Meta分析結(jié)果 表2顯示，采用0.5對0事件進(jìn)行校正后，NNM模型的Meta分析結(jié)果顯示，研究組每1人年發(fā)生事件數(shù)0.719，對照組為0.564；而PNM模型研究組為0.594，對照組為0.419。正如前述，經(jīng)過校正，研究組NNM模型每1人年事件發(fā)生數(shù)比PNM模型增加了17.4%，而對照組增加了25.7%，提示采用0.5校正的方法明顯增加偏倚。此外，采用0事件研究刪除的方法顯示，NNM模型的Meta分析總體變化不大，以研究組為例，刪除0事件研究后，每1人年發(fā)生事件數(shù)為0.721，與采用連續(xù)校正法相比，差異僅0.27%。但對于PNM模型，則變化較大，同樣以研究組為例，刪除0事件研究后，每1人年發(fā)生事件數(shù)為0.676，與采用連續(xù)校正法相比，差異仍達(dá)12.1%，提示刪除0事件研究對PNM模型影響較大，反過來提示PNM具有較好的靈敏性。

表2 不同干預(yù)組總體事件發(fā)生率的Meta分析結(jié)果

3.2 兩種干預(yù)對總體事件發(fā)生率的影響的Meta分析結(jié)果 表3顯示，運(yùn)行附錄代碼3后可以獲得NNM模型和BNM模型的Meta分析結(jié)果。

表3 不同干預(yù)組總體事件發(fā)生率的Meta分析結(jié)果

在表3中，NNM模型的Meta分析結(jié)果顯示，采用連續(xù)性校正和0事件研究刪除法，對效應(yīng)值OR影響率為7.7%，影響率計(jì)算方法為[(2.434-2.247)/2.434]×100%，而BNM模型則高達(dá)25.70%{[(3.171-2.356)/2.356]×100%}。在該示例中，BNM模型給出總發(fā)生率的處理效應(yīng)在幅度上稍微大于NNM模型，并且P值較小：NNM為0.0489，BNM方法為0.0118。BNM模型的效應(yīng)值更高，而P值則更小，因此結(jié)果更靈敏、準(zhǔn)確。根據(jù)表3中BNM模型的結(jié)果，可以認(rèn)為VEGFR-TKIs治療癌癥患者發(fā)生致命不良事件(FAE)的風(fēng)險(xiǎn)大于對照組，其效應(yīng)值OR=3.171，SE=1.502,95% CI：1.383～7.272。在這個(gè)例子中，BNM和NNM方法給出的跨研究標(biāo)準(zhǔn)差相似，但均<0.001，提示研究間差異很小(讀者可以運(yùn)行附錄代碼進(jìn)行查看，這里不贅述)。

4 討論

目前，關(guān)于發(fā)病率數(shù)據(jù)(隨訪資料中的事件數(shù)量)的Meta分析方法并未引起足夠的重視[10]。當(dāng)數(shù)據(jù)中存在0事件時(shí)，如何進(jìn)行Meta分析仍然是一個(gè)懸而未決的問題[11]。盡管如此，越來越多的臨床小樣本研究需要進(jìn)行Meta分析，以獲得高質(zhì)量的證據(jù)，尤其對某些新藥開發(fā)中少見或罕見不良反應(yīng)的發(fā)生率分析[11]。Spittal等[1]早已證明，當(dāng)納入分析的研究存在0事件時(shí)，采用倒方差方法獲得的匯總結(jié)果是有偏差的，而這種方法卻是最常用和最受推薦的方法之一。本文使用Cochrane手冊推薦的連續(xù)性校正進(jìn)行調(diào)整后，結(jié)果仍然存在偏倚風(fēng)險(xiǎn)[12]。

Guevara等提出使用泊松回歸對每個(gè)研究的效應(yīng)指標(biāo)進(jìn)行分析[3]。由于泊松分布包括0，對于含0事件研究而言，這種分布基礎(chǔ)表明可以基于泊松分布進(jìn)行發(fā)病率數(shù)據(jù)的Meta分析[1]。采用泊松分布進(jìn)行Meta分析的關(guān)鍵問題是，基線變異性和異質(zhì)性如何影響最終的總體IRR和其他固定效應(yīng)估計(jì)[4]。理論上，該方法可以通過使用每個(gè)研究效應(yīng)指標(biāo)來解釋基線變異性(即對照組中發(fā)病率的變化)。模擬數(shù)據(jù)表明，當(dāng)基線變異性較低時(shí)，基于泊松分布的Meta分析方法可以產(chǎn)生相對無偏的估計(jì)[1, 5]，證明這種方法是行之有效的。本文使用Proc NLMIXED來實(shí)現(xiàn)基于PNM的Meta分析技術(shù)，對示例數(shù)據(jù)進(jìn)行分析。在這個(gè)例子中，基于NNM模型和PNM模型擬合出的總體事件發(fā)生率存在比較大的的差異。正如預(yù)期的那樣，采用連續(xù)性校正的傳統(tǒng)方法的總體發(fā)病率略高，對照組和治療組分別高25.7%和17.4%，進(jìn)一步證明提示采用連續(xù)性校正的方法明顯增加偏倚。

本文同時(shí)給出刪除0事件研究的Meta分析結(jié)果。在實(shí)踐中，這種方法是不推薦使用的[13]。對于研究組和對照組同時(shí)包含0事件時(shí)，通常從Meta分析中排除(這是許多統(tǒng)計(jì)軟件的默認(rèn)選項(xiàng))。Cochrane手冊認(rèn)為，此類研究并未包含相關(guān)的OR/RR的信息，因此應(yīng)予以排除[12]。然而，一些研究人員指出，從倫理的角度來看，雙0研究中的患者應(yīng)該被納入分析[14]；而另一些人認(rèn)為，這些研究可能通過樣本量來提供相對治療效果的信息[15]。筆示例分析表明，在BNM模型背景下，通過刪除含0事件的研究，最終匯總效應(yīng)值差別可高達(dá)25.7%，顯然這是不可取的。

本文側(cè)重在介紹Stijnen等提出的模型方法。結(jié)合提供的示例和Stijnen等的示例，證明基于BNM模型、PNM模型的Meta分析模型在SAS軟件中可以很容易擬合。本文謹(jǐn)慎地推薦，對于發(fā)病率數(shù)據(jù)Meta分析，GLMM(即本文所述的BNM模型、PNM模型)應(yīng)該成為優(yōu)選技術(shù)。

附錄代碼1：資料錄入

data Example_data;

input studyy1n1t1y0n0t0;

t0=t0/12; t1=t1/12; /*代表1人年的結(jié)果*/

datalines;

1 0 149 6.5 0 75 4.2

……

;

/*連續(xù)性校正*/

data Example_data_CR;

set Example_data;

if (y1=0) or (y0=0) then do;y1 = y1+0.5;y0 = y0+0.5;end;

/* 研究組發(fā)病率的正態(tài)-正態(tài)模型 */

proc nlmixed data=Example_data_CR qpoints=100;

parms theta= 2.0 thetai= 2.0 sigmasq=0.3;

logodds=log(y1/t1);

sesq=1/y1;

model logodds ～ normal(thetai,sesq);

random thetai ～ normal(theta, sigmasq) subject=study;

run;

/* 研究組發(fā)病率的泊松-正態(tài)模型 */

proc nlmixed data=Example_data qpoints=100;

parms theta=1.5 sigmasq=0.3;

model y1 ～ poisson(t1*exp(thetai));

random thetai ～ normal(theta,sigmasq) subject=study;

run;

/* 對數(shù)事件發(fā)生率比值比的正態(tài)-正態(tài)模型*/

proc nlmixed data=Example_data_CR qpoints=100;

parms theta= 0.5 to 2.0 by 0.5 thetai= 1 sigmasq=0.3;

logodds=log((y1/t1)/(y0/t0));

sesq=1/y1+1/y0;

model logodds ～ normal(thetai,sesq);

random thetai ～ normal(theta, sigmasq) subject=study;

run;

/*對數(shù)事件發(fā)生率比值比的二項(xiàng)式-正態(tài)模型*/

proc nlmixed data=Example_data qpoints=100;

parms theta=0.5 to 2.0 by 0.5 sigmasq=0.3;

model y1 ～ binomial(y0+y1, 1/(1+t0/t1*exp(-thetai)));

random thetai ～ normal(theta,sigmasq) subject=study;

run;