肖 靜 吳桂云 高月霞 沈 毅 李奕辰 陳德喜 蔡 輝△

Meta分析是對(duì)以往的研究結(jié)果進(jìn)行系統(tǒng)定量綜合的統(tǒng)計(jì)學(xué)方法[1]。近年來(lái)，隨著循證醫(yī)學(xué)的迅速興起，Meta分析研究的結(jié)論在有關(guān)疾病或醫(yī)學(xué)健康問(wèn)題的防治、治療等方面越來(lái)越受到重視[2]。其中單一率資料的Meta分析是一種只提供人群總?cè)藬?shù)和事件發(fā)生人數(shù)，不同于其它有對(duì)照人群的Meta分析，大多為患病率、檢出率、病死率和感染率等[3]，屬于橫斷面研究，是揭示暴露與疾病關(guān)系不可或缺的，因此基于計(jì)算機(jī)的實(shí)現(xiàn)是很有必要的[4]。本文運(yùn)用實(shí)例數(shù)據(jù)基于Matlab軟件編寫了Meta分析以及發(fā)表偏倚和異質(zhì)性檢驗(yàn)實(shí)現(xiàn)的相關(guān)程序，供讀者參考。

實(shí)例分析

數(shù)據(jù)來(lái)自英國(guó)萊切斯特大學(xué)Eaden等[5]對(duì)大腸潰瘍患者發(fā)生結(jié)直腸癌風(fēng)險(xiǎn)的Meta分析，摘錄年齡組為11～12歲，數(shù)據(jù)見(jiàn)表1。

表1 患大腸潰瘍?nèi)巳喊l(fā)生結(jié)直腸癌風(fēng)險(xiǎn)Meta分析數(shù)據(jù)(部分)

1.Matlab程序?qū)崿F(xiàn)過(guò)程：

(1)建立數(shù)據(jù)集并計(jì)算各研究效應(yīng)量：

clear all;

clc;

s=[8 54

12 2295

9 429

……];%錄入數(shù)據(jù)，第一列為各個(gè)研究事件發(fā)生例數(shù)，第二列為各個(gè)研究總個(gè)體數(shù)。

r=s(:,1)./s(:,2);%取出s矩陣中的第一列除以第二列，計(jì)算各個(gè)研究事件的發(fā)生率ri。

a=s(:,1);%定義s矩陣中第一列為ai。

n=s(:,2);%定義s矩陣中第二列數(shù)據(jù)為ni。

Dr=exp(log(r)-1.96*se);%計(jì)算各個(gè)研究發(fā)生率95%的下限，注意Matlab軟件中l(wèi)og表示ln。

Ur=exp(log(r)+1.96*se);%計(jì)算各個(gè)研究發(fā)生率95%的上限。

w=1./se.^2;%各個(gè)研究的權(quán)重。

Dhr=exp(log(hr)-1.96*sehr);

Uhr=exp(log(hr)+1.96*sehr);

Q=sum(w.*(log(r)-log(hr)).^2);%Q檢驗(yàn)。

k=size(s,1);

Qp=1-chi2cdf(Q,k-1);%計(jì)算Q檢驗(yàn)的p值。

I2=(Q-(k-1))/Q*100;%計(jì)算I2統(tǒng)計(jì)量。

If I2<0

I2=0;

楊琳說(shuō)，以前朋友的老公。以前朋友是指歐陽(yáng)橘紅。聽楊琳這一說(shuō)，一種從來(lái)不曾有的溫暖，流入他的心中。老天爺不睜開眼睛看一看，這樣的好人，為什么偏偏得這種惡?。繛槭裁春萌司兔婚L(zhǎng)呢？

end

If I2>0

I2=I2;

end

I2

Z=abs(log(hr)/sehr);%I2檢驗(yàn)統(tǒng)計(jì)量值的計(jì)算。

Zp=2*(1-normcdf(Z,0,1));%I2檢驗(yàn)統(tǒng)計(jì)量對(duì)應(yīng)的P值。

H=(Q/(k-1))^0.5;%H檢驗(yàn)統(tǒng)計(jì)量的計(jì)算。

HSE=0.5.*(log(Q)-log(k-1))./((2*Q)^0.5-(2*k-3)^0.5;%H檢驗(yàn)統(tǒng)計(jì)量的標(biāo)準(zhǔn)誤。

Hs=exp(log(H)+1.96*HSE);

Hx=exp(log(H)-1.96*HSE);

If Qp<=0.1;%當(dāng)Q檢驗(yàn)p<=0.1時(shí)進(jìn)行隨機(jī)效應(yīng)模型。

If Q<=k-1

tau=0;

end

If Q>k-1

tau=(Q-(k-1))./(sum(w)-(sum(w.^2)./(sum(w))));%隨機(jī)效應(yīng)模型中τ2的計(jì)算。

end

tau

wj=1./((se).^2+tau);%w′的計(jì)算。

Zj=abs(log(hrj)/sehrj);%z′的計(jì)算。

Zpj=2*(1-normcdf(Zj,0,1));%z′對(duì)應(yīng)的P值的計(jì)算。

End

(3)森林圖：此處由于森林圖的程序繁雜略去，感興趣的讀者可以聯(lián)系本文作者索取，結(jié)果見(jiàn)圖1所示。

(4)發(fā)表偏倚的檢驗(yàn)：繪制漏斗圖以及對(duì)漏斗圖進(jìn)行檢驗(yàn)(Egger′s檢驗(yàn))，相關(guān)結(jié)果見(jiàn)圖2、圖3。

Figure;%開始畫漏斗圖。

y=se;%以精度的倒數(shù)為縱坐標(biāo)。

plot(x,y,′o′);%繪制圖。

hold on

plot([log(hrj)log(hrj)],[0 1.2],′k-′)

plot([log(hrj)log(hrj)-1.96*1.2],[0 1.2],′k′)

plot([log(hrj)log(hrj)+1.96*1.2],[0,1.2],′k:′)

set(gca,′tickdir′,′out′,′ydir′,′reverse′,′fontsize′,15,′linewidth′,1.5,′xlim′,[-8-2],′ylim′,[0,1.2],′xtick′,[-8,-6,-4,-2],′xticklabel′,[{num2str(exp(-8))},{num2str(exp(-6))},{num2str(exp(-4))},{num2str(exp(-2))}]);%設(shè)置坐標(biāo)軸。

xlabel(′log(r)′,′fontsize′,16);%給X軸加標(biāo)簽。

ylabel(′se′,′fontsize′,16);%給Y軸加標(biāo)簽。

box off

%漏斗圖開始檢驗(yàn)。

Y=log(r)./se;%定義Y數(shù)值。

X=1./se;%定義X數(shù)值。

b=sum((X-mean(X)).*(Y-mean(Y)))./sum((X-mean(X)).^2);%計(jì)算斜率b。

a=mean(Y)-b.*mean(X);%計(jì)算截距a。

figure

plot(X,Y,′O′,′markersize′,5);

hold on

plot([0,max(X)],[a,a+b.*max(X)],′k-′);

Yi=a+b.*X;%建立回歸方程。

Syx=(sum((Y-Yi).^2)./(size(X,1)-2)).^0.5;

Sa=Syx.*(sum((X).^2)./(size(X,1).*sum((X-mean(X)).^2))).^0.5;

t=abs((a-0)/Sa;%對(duì)a進(jìn)行t檢驗(yàn)，此處計(jì)算t檢驗(yàn)統(tǒng)計(jì)量。

P=2*(1-tcdf(t,size(X,1)-2));%計(jì)算t檢驗(yàn)統(tǒng)計(jì)量對(duì)應(yīng)的自由度k-2下的P值。

tinv(0.975,size(X,1)-2);

ax=a-tinv(0.975,size(X,1)-2)*Sa;

as=a+tinv(0.975,size(X,1)-2)*Sa;

set(gca,′tickdir′,′out′,′fontsize′,15,′linewidth′,1.5,′xlim′,[-0.18],′ylim′,[-328],′xtick′,[1:2:8],′ytick′,[-32:4:8]);%定義坐標(biāo)軸。

xlabel(′1/se′,′fontsize′,16);

ylabel(′log(r)/se′,′fontsize′,16);

box off

plot([0 0],[axas],′b-′,′linewidth′,1);%畫靠近y的線。

plot(0,a,′rs′,′markersize′,6,′markerfacecolor′,′w′);%畫區(qū)間線中間的正方形框，為了直觀表現(xiàn)。

plot([0-0.05 0+0.05],[axax],′k-′,′linewidth′,1);

plot([0-0.05 0+0.05],[asas],′k-′,′linewidth′,1);

(5)異質(zhì)性檢驗(yàn)：繪制Galbraith圖，結(jié)果見(jiàn)圖4。

figure

plot(X,Y,′.′,′markersize′,8);

hold on

plot([0,max(X)],[0,log(hr).*max(X)],′k-′);

hold on

plot([0,max(X)],[2,log(hr).*max(X)+2],′k-′);%在上兩個(gè)單位處畫等斜率直線。

plot([0,max(X)],[-2,log(hr).*max(X)-2],′k-′);

set(gca,′tickdir′,′out′,′fontsize′,15,′linewidth′,1.5,′xlim′,[0 6.5],′ylim′,[-35 2]);

xlabel(′1/se′,′fontsize′,16);

ylabel(′log(r)/se′,′fontsize′,16);

box off

for i=1:size(X,1)

text(X(i)+0.02,Y(i),num2str(i));%此處為給每個(gè)點(diǎn)加名稱。

end

注：以上程序?yàn)橥ㄓ贸绦?，可以根?jù)異質(zhì)性檢驗(yàn)結(jié)果自行選擇計(jì)算模型；程序中加冒號(hào)表示結(jié)果中隱藏不顯示，百分號(hào)后面的為注釋部分，在程序中不運(yùn)行。

圖1 森林圖

圖2 漏斗圖

圖3 Egger′s檢驗(yàn)圖

圖4 Galbraith圖

討 論

隨著循證醫(yī)學(xué)的發(fā)展，Meta分析已被公認(rèn)為客觀評(píng)價(jià)和合成針對(duì)某一特定問(wèn)題研究證據(jù)的最佳手段[3]。對(duì)于目前國(guó)內(nèi)通用的Meta分析軟件RevMan不能進(jìn)行Meta回歸分析、累積Meta分析和Egger′s檢驗(yàn)等，Stata軟件雖為目前Meta分析最受推崇的軟件，但無(wú)單一率Meta分析的固定模塊。本文基于Matlab軟件編寫了單一率資料Meta分析的實(shí)現(xiàn)程序，非專業(yè)人員只需替換數(shù)據(jù)即可方便實(shí)現(xiàn)。Matlab人際交互性強(qiáng)，具有強(qiáng)大的數(shù)據(jù)處理功能[6]，且分析結(jié)果的計(jì)算精度可以通過(guò)程序自由控制，同時(shí)Matlab具有強(qiáng)大的作圖功能，所作圖形比同類軟件更為美觀。本文通過(guò)實(shí)例分析驗(yàn)證了利用Matlab程序?qū)崿F(xiàn)單一率Meta分析結(jié)果的可行性和有效性，并完成了異質(zhì)性和發(fā)表偏倚的檢驗(yàn)。

參 考 文 獻(xiàn)

1.Egger M,Smith G D,Altman D G.Systematic reviews in health care.Meta-analysis in context.2nd edition.London:BMJ Publishing Group,2001.

2.俞慧強(qiáng),鄭輝烈,李悅,等.Meta分析發(fā)表偏倚診斷方法研究.中國(guó)衛(wèi)生統(tǒng)計(jì),2011,28(4):402-405.

3.曾憲濤,冷衛(wèi)東,李勝,等.如何正確理解及使用GRADE系統(tǒng).中國(guó)循證醫(yī)學(xué)雜志,2011,11(9):985-990.

4.王佩鑫,李宏田,劉建蒙.無(wú)對(duì)照二分類資料的Meta分析方法及Stata實(shí)現(xiàn).循證醫(yī)學(xué),2012,12(1):52-64.

5.Eaden J A,Abrams K R,Mayberry J F.The risk of colorectal cancer in ulcerative colitis:a meta-analysis.Gut,2001,48(4):526-535.

6.胡小剛,陳劍鴻,孫鳳軍,等.基于Matlab的Kruskalk-Wallis和Nemenyi檢驗(yàn)的界面實(shí)現(xiàn).中國(guó)衛(wèi)生統(tǒng)計(jì),2011,28(4):466-473.