單一率Meta分析的Matlab軟件實現(xiàn)*

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

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

實例分析

數(shù)據(jù)來自英國萊切斯特大學Eaden等[5]對大腸潰瘍患者發(fā)生結直腸癌風險的Meta分析，摘錄年齡組為11～12歲，數(shù)據(jù)見表1。

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

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

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

clear all;

clc;

s=[8 54

12 2295

9 429

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

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

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

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

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

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

w=1./se.^2;%各個研究的權重。

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

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

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

k=size(s,1);

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

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

If I2<0

I2=0;

楊琳說，以前朋友的老公。以前朋友是指歐陽橘紅。聽楊琳這一說，一種從來不曾有的溫暖，流入他的心中。老天爺不睜開眼睛看一看，這樣的好人，為什么偏偏得這種惡病？為什么好人就命不長呢？

end

If I2>0

I2=I2;

end

I2

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

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

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

HSE=0.5.*(log(Q)-log(k-1))./((2*Q)^0.5-(2*k-3)^0.5;%H檢驗統(tǒng)計量的標準誤。

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

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

If Qp<=0.1;%當Q檢驗p<=0.1時進行隨機效應模型。

If Q<=k-1

tau=0;

end

If Q>k-1

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

end

tau

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

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

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

End

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

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

Figure;%開始畫漏斗圖。

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

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))}]);%設置坐標軸。

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

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

box off

%漏斗圖開始檢驗。

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

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

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

a=mean(Y)-b.*mean(X);%計算截距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;%對a進行t檢驗，此處計算t檢驗統(tǒng)計量。

P=2*(1-tcdf(t,size(X,1)-2));%計算t檢驗統(tǒ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]);%定義坐標軸。

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ì)性檢驗：繪制Galbraith圖，結果見圖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-′);%在上兩個單位處畫等斜率直線。

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));%此處為給每個點加名稱。

end

注：以上程序為通用程序，可以根據(jù)異質(zhì)性檢驗結果自行選擇計算模型；程序中加冒號表示結果中隱藏不顯示，百分號后面的為注釋部分，在程序中不運行。

圖1 森林圖

圖2 漏斗圖

圖3 Egger′s檢驗圖

圖4 Galbraith圖

討 論

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

參 考 文 獻

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.俞慧強,鄭輝烈,李悅,等.Meta分析發(fā)表偏倚診斷方法研究.中國衛(wèi)生統(tǒng)計,2011,28(4):402-405.

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

4.王佩鑫,李宏田,劉建蒙.無對照二分類資料的Meta分析方法及Stata實現(xiàn).循證醫(yī)學,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檢驗的界面實現(xiàn).中國衛(wèi)生統(tǒng)計,2011,28(4):466-473.