SAS軟件是一款功能強大、國際認可度高的編程軟件，其可以實現包括網狀Meta分析在內的各種Meta分析。貝葉斯統計是一種重要的統計學方法，其可以通過馬爾科夫-蒙特卡洛（Markov Chain Monte Carlo，MCMC）算法進行各種統計推斷。在SAS軟件中運用PROC MCMC過程實現網狀Meta分析正是基于這一思想進行的，本文通過實例操作展示了該過程。
引用本文： 鄢金柱, 吳平勇, 張超, 劉俊, 曾憲濤. 應用SAS軟件PROC MCMC過程實現網狀Meta分析. 中國循證醫學雜志, 2014, 14(8): 1004-1008. doi: 10.7507/1672-2531.20140165 復制
網狀Meta分析是近年隨臨床需求出現的新的Meta分析方法[1, 2]，當前有較多軟件可實現，經典的如WinBUGS和OpenBUGS [3-5]，其他包括Stata [6]、R [7, 8]、GeMTC [9, 10]、JAGS [11]、ADDIS [12]及Microsoft Excel [13]等。此外，SAS軟件亦可以實現網狀Meta分析。
SAS軟件是諸多編程軟件中的一種，基于其強大的數據統計分析功能，其可以進行多種統計方法的網狀Meta分析，尤其是在大樣本數據分析中具有很強的優勢。目前已有文獻對SAS實現診斷性試驗的Meta分析[14, 15]和頭對頭比較的Meta分析[16-18]方法進行了介紹。本文將繼續以前文[5-13]介紹網狀Meta分析的示例為例，介紹網狀Meta分析在SAS中的實現。因SAS需要針對不同的情況進行編程且程序代碼較長，故為方便展示，我們選擇14種抗抑郁藥中的4種作為范例，展示如何使用SAS軟件PROC MCMC過程[19]進行網狀Meta分析，以期為讀者提供參考。
本文所用SAS為9.3版。
1 準備數據
在使用SAS進行統計分析前，我們需對數據進行整理。本文提取的4種藥物為placebo、fluoxetine、bupropion和sertaline [5-13]，整理后的數據見表 1。
表1
本文例題使用的數據
2 編寫程序
SAS是一款編程軟件，數據運算的每一步都需要通過代碼來實現。本文應用PROC MCMC過程及隨機效應模型來實現網狀Meta分析[19]，編寫的代碼如下，同時我們用“/* */”符號為每一步代碼做了注解。
Data betablock1;/*建立名稱為“betablock1”的數據集*/
;
run;
/*“Itrt”編入索引是為了在proc MCMC過程中使用(1=A, 2=B, 3=C, 4=D); */
/* “Recordis”是一個在MCMC隨機效應模型中使用的獨立、連續的標示符*/
data betablock2;
set betablock1 end=myend;
length Trt $ 8;
retain Nreff 0 NTrtlev 0 Nstudy 0;
drop ar an br bn cr cn dr dn Ntrtlev Record Nreff NStudy;
retain record 0;
Narm=(ar+an > 0)+(br+bn > 0)+(cr+cn > 0)+(dr+dn > 0);
if ar+an > 0 then do; r=ar; n=an; Trt="A"; ITrt=1; record=record+1; output; end;
if br+bn > 0 then do; r=br; n=bn; Trt="B"; ITrt=2; record=record+1; output; end;
if cr+cn > 0 then do; r=cr; n=cn; Trt="C"; ITrt=3; record=record+1; output; end;
if dr+dn > 0 then do; r=dr; n=dn; Trt="D"; ITrt=4; record=record+1; output; end;
Nstudy=Nstudy+1; Nreff=Nreff+6-1; NTrtlev=max(Ntrtlev, 6); /*標記1*/
if myend then do;
??call symput("Nrec", record);
??call symput("Nreff", Nreff);
??call symput("NTrtlev", NTrtlev);
??call symput("NStudy", NStudy);
end;
run;
??
/*去除各數據前端的空格*/
%let Ntrtlev= & Ntrtlev; %let Nstudy= & Nstudy; %let Nreff= & Nreff; %let Nrec= & Nrec;
/*展示各個數據*/
%put Nrec= & Nrec; %put Nreff= & Nreff; %put NTrtlev= & NTrtlev; %put NStudy= & NStudy;
???
/*產生虛擬數據保證程序順利計算*/
data a; run;
???
/*隨機效應模型計算*/
Title1 "隨機效應模型PROC MCMC過程";
/*下面的代碼將運行PROC MCMC過程并使所有數據讀入矩陣且保證不會重復讀取*/
proc mcmc data=a monitor=(beta2-beta & Ntrtlev LogOR OR sd v)STATISTICS(PERCENTAGE=(2.5 5 25 50 75 95 97.5))=ALL
NTU=10000 NBI=10000 NMC=1000000 thin=100 jointmodel;
array beta[ & Ntrtlev] beta1-beta & Ntrtlev;
array alpha[ & Nstudy] alpha1-alpha & Nstudy;
array LogOR[ & Ntrtlev];
array OR[ & Ntrtlev];
??
array randeff[ & NReff] Randeff1-Randeff & NReff;
array reffwt[ & Ntrtlev, %eval( & Ntrtlev-1)];
array reffind[ & Nrec];
array reffpt[ & Nrec];
??
/*內部存儲的數據*/
array Study[1]/nosymbols; array r[1]/nosymbols; array n[1]/nosymbols;
array itrt[1]/nosymbols; array Narm[1]/nosymbols;
?
/*固定效應模型約束*/
BEGINCNST;
/*讀取數據進入工作空間*/
??rc=read_array("betablock2", Study, "Study");
??rc=read_array("betablock2", R, "r");
??rc=read_array("betablock2", N, "n");
??rc=read_array("betablock2", itrt, "itrt");
??rc=read_array("betablock2", Narm, "Narm");
??Nrec=dim(Study, 1);
??/*k是指定當前的研究*/
??k=0;
??reffpoint=0;
??kk=1;
??do i=1 to & nrec;
????/*抵消該研究的第一個隨機效應; 檢查是否有重復納入研究*/
????if k=Study[i] then do; kk=kk+1;end;
????else do; k=Study[i]; reffpoint=reffpoint+kk-1;kk=1;end;
????reffind[i]=reffpoint;
????reffpt[i]=kk;
????end;
?
????/*此處定義約束條件; *該部分被設計得很靈活, 以至可用于任何不同的約束條件*/
????do i=1 to & Ntrtlev;
????do j=1 to i-2; reffwt[i, j]=sqrt(1/(2*j*(j+1))); end;
????if i > 1 then reffwt[i, i-
????do j=i to %eval( & Ntrtlev-1); reffwt[i, j]=0; end;
??end;
??******************;
??beta[1]=0;
?
ENDCNST;
?
/*指定先驗和初始值*/
parms alpha: 0;
parms beta2-beta & Ntrtlev 0;
parms sd 1;
parms Randeff: 0;
?
prior sd~uniform(0.001, 10);
v=(sd*sd);
?
/*保證先驗的合理性*/
prior Randeff:~normal(0, var=v);
prior alpha:~normal(0, var=10000);
prior beta2-beta & Ntrtlev~normal(0, var=10000);
?
ll=0;
do i=1 to & nrec;
/*二進制數據; 這里的隨機效應是該試驗不受約束的隨機效應的加權總和*/
??reff=0;
??do j=1 to Narm[i]-1;
????reff=reff+reffwt[reffpt[i], j]*randeff[reffind[i]+j];
??end;
l??inpred=alpha[study[i]]+beta[itrt[i]]+reff;
??p=1-(1/(1+exp(linpred)));
??ll=ll+r[i]*log(p)+(n[i]-r[i])*log(1-p);
end;
?
model general(ll);
/*標記2*/
LogOR[1]=(beta[2]-beta[1]); *B-A;
LogOR[2]=(beta[3]-beta[1]); *C-A;
LogOR[3]=(beta[4]-beta[1]); *D-A;
LogOR[4]=(beta[3]-beta[2]); *C-B;
LogOR[5]=(beta[4]-beta[2]); *D-B;
LogOR[6]=(beta[4]-beta[3]); *D-C;
?
OR[1]=exp((beta[2]-beta[1])); *B-A;
OR[2]=exp((beta[3]-beta[1])); *C-A;
OR[3]=exp((beta[4]-beta[1])); *D-A;
OR[4]=exp((beta[3]-beta[2])); *C-B;
OR[5]=exp((beta[4]-beta[2])); *D-B;
OR[6]=exp((beta[4]-beta[3])); *D-C;
ODS output postsummaries=library.post1BIG;
run;
ODS output clear;
?
以上為完整的代碼，點擊運行即可開始數據計算并得出結果。
SAS軟件的特點是根據數據進行代碼編寫，因此，應用此套代碼時要注意結合數據本身的計算要求做出修改：①建立數據集。這是進行數據計算的第一步，也是最重要的一步，本套代碼我們建立了9個變量，而讀者進行不同研究的計算時需根據自己的數據情況來建立數據集。②本文示例數據涉及4種藥物，共6種比較結果，在“標記1”處即編程為6，對于其他數目的比較結果需做相應更改；同時根據數據情況，在“標記2”處需作相應更改。
3 結果展示
編程好上述代碼后，提交運算，即可得出結果（表 2）。
表2
SAS軟件PROC MCMC過程計算的示例結果
此外，我們使用了WinBUGS軟件[5]對表 1中的數據進行計算，結果見表 3。對比表 2與表 3，其結果基本相同，這與迭代次數的設置及軟件間的差異有關。
表3
WinBUGS軟件計算的示例結果
4 結語
SAS軟件使用廣泛，使用不同的編碼可靈活處理各種數據[14-19]，本文介紹了通過PROC MCMC過程、隨機效應模型來計算網狀Meta分析。本文提供的代碼程序是在Jones等[19]提供的SAS軟件進行網狀Meta分析的PROC MCMC過程代碼上修改的，原文使用的是3種藥物[19]，本文使用的是4種藥物，讀者可以比較兩者的代碼異同。
因此，在處理3種藥物時，可以直接使用Jones等[19]提供的代碼，而在處理4種藥物時，可以使用本文提供的代碼；在處理其他數據時，我們只需根據具體情況稍作修改，即可完成網狀Meta分析。同時，從本文的代碼可以看出，代碼很長，這也是為何我們僅選擇了4種藥物作為演示，未使用全部14種藥物的原因；另一原因是為了可以讓讀者比較與3種藥物的代碼，便于學習。
貝葉斯理論有著廣泛的應用，除了用SAS的PROC MCMC過程實現網狀Meta分析外，還可以用Stata、WinBUGS和R等軟件來實現[
此外，可以看到，SAS的PROC MCMC過程并不能繪制網狀Meta分析相關的森林圖、網狀關系圖等圖形，需要另行編寫代碼來實現，亦可通過其他軟件[20]來實現。
網狀Meta分析是近年隨臨床需求出現的新的Meta分析方法[1, 2]，當前有較多軟件可實現，經典的如WinBUGS和OpenBUGS [3-5]，其他包括Stata [6]、R [7, 8]、GeMTC [9, 10]、JAGS [11]、ADDIS [12]及Microsoft Excel [13]等。此外，SAS軟件亦可以實現網狀Meta分析。
SAS軟件是諸多編程軟件中的一種，基于其強大的數據統計分析功能，其可以進行多種統計方法的網狀Meta分析，尤其是在大樣本數據分析中具有很強的優勢。目前已有文獻對SAS實現診斷性試驗的Meta分析[14, 15]和頭對頭比較的Meta分析[16-18]方法進行了介紹。本文將繼續以前文[5-13]介紹網狀Meta分析的示例為例，介紹網狀Meta分析在SAS中的實現。因SAS需要針對不同的情況進行編程且程序代碼較長，故為方便展示，我們選擇14種抗抑郁藥中的4種作為范例，展示如何使用SAS軟件PROC MCMC過程[19]進行網狀Meta分析，以期為讀者提供參考。
本文所用SAS為9.3版。
1 準備數據
在使用SAS進行統計分析前，我們需對數據進行整理。本文提取的4種藥物為placebo、fluoxetine、bupropion和sertaline [5-13]，整理后的數據見表 1。
表1
本文例題使用的數據
2 編寫程序
SAS是一款編程軟件，數據運算的每一步都需要通過代碼來實現。本文應用PROC MCMC過程及隨機效應模型來實現網狀Meta分析[19]，編寫的代碼如下，同時我們用“/* */”符號為每一步代碼做了注解。
Data betablock1;/*建立名稱為“betablock1”的數據集*/
;
run;
/*“Itrt”編入索引是為了在proc MCMC過程中使用(1=A, 2=B, 3=C, 4=D); */
/* “Recordis”是一個在MCMC隨機效應模型中使用的獨立、連續的標示符*/
data betablock2;
set betablock1 end=myend;
length Trt $ 8;
retain Nreff 0 NTrtlev 0 Nstudy 0;
drop ar an br bn cr cn dr dn Ntrtlev Record Nreff NStudy;
retain record 0;
Narm=(ar+an > 0)+(br+bn > 0)+(cr+cn > 0)+(dr+dn > 0);
if ar+an > 0 then do; r=ar; n=an; Trt="A"; ITrt=1; record=record+1; output; end;
if br+bn > 0 then do; r=br; n=bn; Trt="B"; ITrt=2; record=record+1; output; end;
if cr+cn > 0 then do; r=cr; n=cn; Trt="C"; ITrt=3; record=record+1; output; end;
if dr+dn > 0 then do; r=dr; n=dn; Trt="D"; ITrt=4; record=record+1; output; end;
Nstudy=Nstudy+1; Nreff=Nreff+6-1; NTrtlev=max(Ntrtlev, 6); /*標記1*/
if myend then do;
??call symput("Nrec", record);
??call symput("Nreff", Nreff);
??call symput("NTrtlev", NTrtlev);
??call symput("NStudy", NStudy);
end;
run;
??
/*去除各數據前端的空格*/
%let Ntrtlev= & Ntrtlev; %let Nstudy= & Nstudy; %let Nreff= & Nreff; %let Nrec= & Nrec;
/*展示各個數據*/
%put Nrec= & Nrec; %put Nreff= & Nreff; %put NTrtlev= & NTrtlev; %put NStudy= & NStudy;
???
/*產生虛擬數據保證程序順利計算*/
data a; run;
???
/*隨機效應模型計算*/
Title1 "隨機效應模型PROC MCMC過程";
/*下面的代碼將運行PROC MCMC過程并使所有數據讀入矩陣且保證不會重復讀取*/
proc mcmc data=a monitor=(beta2-beta & Ntrtlev LogOR OR sd v)STATISTICS(PERCENTAGE=(2.5 5 25 50 75 95 97.5))=ALL
NTU=10000 NBI=10000 NMC=1000000 thin=100 jointmodel;
array beta[ & Ntrtlev] beta1-beta & Ntrtlev;
array alpha[ & Nstudy] alpha1-alpha & Nstudy;
array LogOR[ & Ntrtlev];
array OR[ & Ntrtlev];
??
array randeff[ & NReff] Randeff1-Randeff & NReff;
array reffwt[ & Ntrtlev, %eval( & Ntrtlev-1)];
array reffind[ & Nrec];
array reffpt[ & Nrec];
??
/*內部存儲的數據*/
array Study[1]/nosymbols; array r[1]/nosymbols; array n[1]/nosymbols;
array itrt[1]/nosymbols; array Narm[1]/nosymbols;
?
/*固定效應模型約束*/
BEGINCNST;
/*讀取數據進入工作空間*/
??rc=read_array("betablock2", Study, "Study");
??rc=read_array("betablock2", R, "r");
??rc=read_array("betablock2", N, "n");
??rc=read_array("betablock2", itrt, "itrt");
??rc=read_array("betablock2", Narm, "Narm");
??Nrec=dim(Study, 1);
??/*k是指定當前的研究*/
??k=0;
??reffpoint=0;
??kk=1;
??do i=1 to & nrec;
????/*抵消該研究的第一個隨機效應; 檢查是否有重復納入研究*/
????if k=Study[i] then do; kk=kk+1;end;
????else do; k=Study[i]; reffpoint=reffpoint+kk-1;kk=1;end;
????reffind[i]=reffpoint;
????reffpt[i]=kk;
????end;
?
????/*此處定義約束條件; *該部分被設計得很靈活, 以至可用于任何不同的約束條件*/
????do i=1 to & Ntrtlev;
????do j=1 to i-2; reffwt[i, j]=sqrt(1/(2*j*(j+1))); end;
????if i > 1 then reffwt[i, i-
????do j=i to %eval( & Ntrtlev-1); reffwt[i, j]=0; end;
??end;
??******************;
??beta[1]=0;
?
ENDCNST;
?
/*指定先驗和初始值*/
parms alpha: 0;
parms beta2-beta & Ntrtlev 0;
parms sd 1;
parms Randeff: 0;
?
prior sd~uniform(0.001, 10);
v=(sd*sd);
?
/*保證先驗的合理性*/
prior Randeff:~normal(0, var=v);
prior alpha:~normal(0, var=10000);
prior beta2-beta & Ntrtlev~normal(0, var=10000);
?
ll=0;
do i=1 to & nrec;
/*二進制數據; 這里的隨機效應是該試驗不受約束的隨機效應的加權總和*/
??reff=0;
??do j=1 to Narm[i]-1;
????reff=reff+reffwt[reffpt[i], j]*randeff[reffind[i]+j];
??end;
l??inpred=alpha[study[i]]+beta[itrt[i]]+reff;
??p=1-(1/(1+exp(linpred)));
??ll=ll+r[i]*log(p)+(n[i]-r[i])*log(1-p);
end;
?
model general(ll);
/*標記2*/
LogOR[1]=(beta[2]-beta[1]); *B-A;
LogOR[2]=(beta[3]-beta[1]); *C-A;
LogOR[3]=(beta[4]-beta[1]); *D-A;
LogOR[4]=(beta[3]-beta[2]); *C-B;
LogOR[5]=(beta[4]-beta[2]); *D-B;
LogOR[6]=(beta[4]-beta[3]); *D-C;
?
OR[1]=exp((beta[2]-beta[1])); *B-A;
OR[2]=exp((beta[3]-beta[1])); *C-A;
OR[3]=exp((beta[4]-beta[1])); *D-A;
OR[4]=exp((beta[3]-beta[2])); *C-B;
OR[5]=exp((beta[4]-beta[2])); *D-B;
OR[6]=exp((beta[4]-beta[3])); *D-C;
ODS output postsummaries=library.post1BIG;
run;
ODS output clear;
?
以上為完整的代碼，點擊運行即可開始數據計算并得出結果。
SAS軟件的特點是根據數據進行代碼編寫，因此，應用此套代碼時要注意結合數據本身的計算要求做出修改：①建立數據集。這是進行數據計算的第一步，也是最重要的一步，本套代碼我們建立了9個變量，而讀者進行不同研究的計算時需根據自己的數據情況來建立數據集。②本文示例數據涉及4種藥物，共6種比較結果，在“標記1”處即編程為6，對于其他數目的比較結果需做相應更改；同時根據數據情況，在“標記2”處需作相應更改。
3 結果展示
編程好上述代碼后，提交運算，即可得出結果（表 2）。
表2
SAS軟件PROC MCMC過程計算的示例結果
此外，我們使用了WinBUGS軟件[5]對表 1中的數據進行計算，結果見表 3。對比表 2與表 3，其結果基本相同，這與迭代次數的設置及軟件間的差異有關。
表3
WinBUGS軟件計算的示例結果
4 結語
SAS軟件使用廣泛，使用不同的編碼可靈活處理各種數據[14-19]，本文介紹了通過PROC MCMC過程、隨機效應模型來計算網狀Meta分析。本文提供的代碼程序是在Jones等[19]提供的SAS軟件進行網狀Meta分析的PROC MCMC過程代碼上修改的，原文使用的是3種藥物[19]，本文使用的是4種藥物，讀者可以比較兩者的代碼異同。
因此，在處理3種藥物時，可以直接使用Jones等[19]提供的代碼，而在處理4種藥物時，可以使用本文提供的代碼；在處理其他數據時，我們只需根據具體情況稍作修改，即可完成網狀Meta分析。同時，從本文的代碼可以看出，代碼很長，這也是為何我們僅選擇了4種藥物作為演示，未使用全部14種藥物的原因；另一原因是為了可以讓讀者比較與3種藥物的代碼，便于學習。
貝葉斯理論有著廣泛的應用，除了用SAS的PROC MCMC過程實現網狀Meta分析外，還可以用Stata、WinBUGS和R等軟件來實現[
此外，可以看到，SAS的PROC MCMC過程并不能繪制網狀Meta分析相關的森林圖、網狀關系圖等圖形，需要另行編寫代碼來實現，亦可通過其他軟件[20]來實現。
