OpenBUGS處理結構方程模型貝葉斯分析在SAS宏程序中的實現*

秦正積 嚴曉玲 沈 毅 肖 靜 何 書 任文龍

貝葉斯方法計算后驗分布函數常常是復雜的、高維的，甚至是非標準形式的，其積分計算十分困難[1-2]。譬如，在行為學、社會學、心理學和醫學中，許多潛變量模型需要用到結構方程模型，而結構方程模型的應用，常常含有高緯復雜的模型，不能顯示表達其分布[3]，直接估計困難，可用貝葉斯估計方法，特別是GIBBS抽樣方法和MCMC(Markov Chain Monte Carlo)估計方法[4-5]。MCMC方法可解決此類問題[6-7]，其實現由專業軟件WINBUGS[8]或OpenBUGS[9]完成。GIBBS 抽樣是構建馬爾科夫鏈最流行的MCMC方法[10-11]，其計算工作量較大，對計算機軟硬件有較高要求。近年來隨著計算機軟硬件的快速發展，貝葉斯統計方法受到研究者的重視，貝葉斯統計已經成為統計學中很熱門的研究課題[12]。

但是WINBUGS或OpenBUGS軟件數據管理功能比較局限。大型數據庫管理，譬如數據挖掘常常需要從海量數據中抽樣獲取，SAS數據管理功能強大，數據管理在數據挖掘中有較大的優勢，如果將兩個軟件合理有序的聯合應用，便于計算機后臺工作，可節省人機對話的時間和人力。本文探討了應用OpenBUGS進行貝葉斯結構方程模型分析的SAS程序實現。

研究設計與方法

研究從簡單實例介紹OpenBUGS貝葉斯分析原理與實現過程，然后構建復雜結構方程模型。應用SAS 9.4和OpenBUGS軟件的GIBBS抽樣方法進行貝葉斯分析，并編寫有關的SAS宏程序從實例展示其后臺實現。

貝葉斯分析原理與實例

1.貝葉斯分析原理與OpenBUGS分析程序基本結構

貝葉斯參數估計可用軟件Winbugs或OpenBUGS實現。OpenBUGS軟件是在WinBUGS軟件基礎上研制的一款開源的貝葉斯統計推斷軟件。OpenBUGS進行MCMC貝葉斯分析需要4個文件，其分別完成不同功能。Model文件提供了OpenBUGS軟件需要的貝葉斯邏輯函數和經驗分布。Data提供了分析用的數據文件，用于構建數據提供的似然函數。Inits提供初始值，對于MCMC抽樣模擬的鏈可以采用一條，也可有多條。如果有多條，每條給予不同的初始值，多條鏈的后驗參數最終會收斂到接近一致。另外還有一個語法文件，提供了分析過程語法和結果的自動保存。OpenBUGS分析基本結構可參閱文獻[8]和[13].

2.軟件OpenBUGS示例程序及其解釋

為了簡明的闡述軟件執行的過程，此處利用軟件的例子進行說明。

示例1取自Gelf等[11]的第6部分，研究者測得30只小鼠5周體重。用Yij表示年齡在Xj時的體重，見表1。

表1 不同年齡小鼠體重測量值統計表

研究假定模型是線性增長曲線的隨機效應，模型假定如下：

yij～Normal(αi+βi(xj-xbar),τc)

αi～Normal(αc,τα)

βi～Normal(βc,τβ)

由此，其模型設計的程序如下：

model

{

for(i in 1∶N){

for(j in 1∶T){

Y[i，j]～dnorm(mu[i，j]，tau.c)

mu[i，j]<-alpha[i]+ beta[i]*(x[j]-xbar)

}

alpha[i]～dnorm(alpha.c，alpha.tau)

beta[i]～dnorm(beta.c，beta.tau)

}

tau.c～dgamma(0.001，0.001)

sigma<-1/sqrt(tau.c)

alpha.c～dnorm(0.0，1.0E-6)

alpha.tau～dgamma(0.001，0.001)

beta.c～dnorm(0.0，1.0E-6)

beta.tau～dgamma(0.001，0.001)

alpha0<-alpha.c-xbar*beta.c

}

初始文件：

list(alpha=c(250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250，250)，

beta =c(6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6，6)，

alpha.c=150，beta.c=10，tau.c=1，alpha.tau=1，beta.tau=1)。

數據文件可參考Gelf等[11]。模型MCMC抽樣10000次后去掉初始1000次后，計算其統計量見表2。

表2 貝葉斯分析參數估計

OpenBUGS處理結構方程模型貝葉斯分析的SAS宏程序

1.結構方程模型實例

本研究可采用不同樣本含量數據，樣本數據文件是從記錄數為60000余條的較大型數據文件中抽樣產生的，然后運用OpenBUGS軟件采用MCMC方法進行貝葉斯統計分析，其計算工作量較大，計算時間較長[2]，研究人員為了統籌安排，常常需要在計算機后臺進行，因此有必要編輯SAS宏程序。

數據文件是有關慢性病的大量人群的調查研究資料，研究從該數據庫進行了不同比例的抽樣，得到相應文件，用于復雜結構方程模型的貝葉斯統計分析。

先驗采用含dirichlet的無信息先驗，其中貝葉斯dirichlet先驗可參閱有關文獻[14-15]。

結構方程模型見圖1。

模型含測量方程載荷參數21個(lam[1-12])，結構方程參數9個(gam[1-9])，應用貝葉斯分析法進行估計。由于本文只關心程序是否能運行，所以不對模型合理性進行探討。

圖1 結構方程模型示意圖

2.SAS宏程序實現

結構方程模型的貝葉斯統計分析在SAS宏中實現，程序包括三個SAS宏程序。第一個宏程序用于從大數據中抽樣產生子數據文件，分別為數據的第一行構成的文件opb1obugfirst，其余部分構成文件opb1obuglast，供下一步SAS宏程序產生貝葉斯分析用的數據文件。第二個SAS宏程序命名為“OPENBUGSIO_model_28u_att1”，通過SAS程序產生貝葉斯分析需要的四個文件，包括模型文件、樣本數據文件、初始鏈數據文件、OpenBUGS批處理文件。最后通過SAS調用OpenBUGS批處理文件，然后運用SAS軟件導出參數估計到excel文件。第三個宏程序將第一個宏或幾個宏嵌套其中，以期能組合不同宏，通過改變宏參數，如改變MCMC迭代抽樣的次數及其樣本含量等參數，實現后臺SAS程序批量分析。其中，程序參考了部分SAS參考書[22-23]與網址[24]和部分人員的程序[5，16]，宏程序及解釋如下：

**********************；

**SAS宏程序1/*從已有的含65016條記錄的大型數據文件essay2中抽樣產生數據文件opb650 obugfirst.txt和opb650 obuglast.txt等文件*/

**********************；

%macro semmc；

%do B=1%to2；/*按比例1%，2%進行隨機抽樣*/

%if &B=1 %then %do；%let rate=0.01；%let SMPLN=1；%end；

%if &B=2 %then %do；%let rate=0.02；%let SMPLN=2；%end；

proc surveyselect data=essay2 out=opb&SMPLN method=sys seed=130032 samprate=&rate ；

run；

DATA opb&SMPLN.obug1；

SET opb&SMPLN；

DROP _；

format _all_ e8.2；

RUN；

data tabvar.opb&SMPLN.obugfirst1 tabvar.opb&SMPLN.obuglast；

set opb&SMPLN.obug1；

if _n_=round(65016/10*&rate，0.1)*10 then output tabvar.opb&SMPLN.obuglast；

else output tabvar.opb&SMPLN.obugfirst1；

run；

data tabvar.opb&SMPLN.obugfirst；

set tabvar.opb&SMPLN.obugfirst1；

com=""；

run；

PROC EXPORT DATA=tabvar.opb&SMPLN.obugfirst

OUTFILE= "F：yanjiuqinchronicdirichlet28uorigiondataopb&SMPLN.obugfirst.txt"

DBMS=TAB REPLACE；

PUTNAMES=NO；

DELIMITER="，"；

RUN；

PROC EXPORT DATA=tabvar.opb&SMPLN.obuglast

OUTFILE= “F：yanjiuqinchronicdirichlet28uorigiondataopb&SMPLN.obuglast.txt”

DBMS=TAB REPLACE；

PUTNAMES=NO；

DELIMITER="，"；

RUN；

%END；/*close the do-loop for sample size condition.*/

%mend semmc；

%semmc；

run；

**********************；

**SAS宏程序2

**********************；

%macro OPENBUGSIO_model_28u_att1(dsout=，count=，exceldir=，outf=，numit=，modelcompile=，n=，p=，dir=，datadir=，dirmodelin=，modelf=，dirmodelopenbug=，data1=，datafirst=，datalast=，data2=，init1=，init2=，init3=，suffix=，batchf=，firstdir=，dirout=)；

**********************；

/*讀入模型文件，并指定邏輯文件名*/

**********************；

filename modelin "&firstdir&dirmodelinmodel&dirmodelin&modelf.txt"；run；

filename mmodel "&datadir&modelf..txt"；run；

/*proc printto file=mmodel new；*/

data _null_；

file mmodel new；

run；

data _null_；

infile modelin；

input；

file mmodel mod；put _infile_；run；

proc printto；run；

**********************；

/*生成貝葉斯分析需要的數據文件*/

**********************；

options pagesize=32767 nodate nocenter nonumber formdlim="mprint symbolgen mlogic；

title1"；

filename dat1 "&datadir&data1..txt"；run；

filename datafrst "&datadir&datafirst..txt"；run；

filename datalast "&datadir&datalast..txt"；run；

filename dat2 "&dir&data2..txt"；run；

data _null_；

file dat1 new；

put@1"list(N=650，P=28，alpha=c(1，1，1，1)，m=4，y=structure(.Data=c("；

data _null_；file dat2 new；put@1")，.Dim=c(&n，&p)))"；run；

data _null_；infile datafrst；input；file dat1 mod；put _infile_；run；

data _null_；infile datalast；input；file dat1 mod；put _infile_；run；

data _null_；infile dat2；input；file dat1 mod；put _infile_；run；

proc printto；run；

**********************；

/*產生貝葉斯分析需要的初始鏈文件*/

**********************；

filename init1 "&datadir&init1..txt"；

filename init2 "&dir&init2..txt"；

filename init3 "&dir&init3..txt"；

filename init "&datadir&dirmodelin&init..txt"；

data _null_；

file init new；

run；

data _null_；infile init1；input；file init mod；put _infile_；run；

data _null_；infile init2；input；file init mod；put _infile_；run；

data _null_；infile init3；input；file init mod；put _infile_；run；

proc printto；run；

**********************；

/*產生貝葉斯分析需要的批處理文件*/

**********************；

options pagesize=32767 nodate nocenter nonumber formdlim="mprint symbolgen mlogic；

title1"；

filename fileout2 "C：Program Files(x86)OpenBUGSOpenBUGS323&batchf..txt"；

data _null_；

file fileout2；

put@1"modelDisplay('log')"；

put@1"modelCheck('&dirmodelopenbug&modelf..txt')"；

put@1"modelData('&dirmodelopenbug&data1..txt')"；

put@1"modelCompile(&modelcompile)"；

put@1"modelInits('&dirmodelopenbug&dirmodelin&init..txt'，3)"；

put@1"modelGenInits()"；

put@1"modelUpdate(&dirout.updates'，'thin'，1，T)"；

put@1"modelSaveState('&dirout.outcomexample.odc')"；

put@1"modelSetRN(5)"；

put@1"modelPrecision(4)"；

put@1"samplesSet('lam')"；

put@1"samplesSet('tlam')"；

put@1"samplesSet('attlam')"；

put@1"samplesSet('u')"；

put@1"samplesSet('p')"；

put@1"samplesBeg(1)"；

put@1"samplesEnd(10000000)"；

put@1"samplesSet('gam')"；

put@1"modelUpdate(&numit)"；

put@1"samplesBeg(10)"；

put@1"samplesEnd(1000000)"；

put@1"samplesFirstChain(1)"；

put@1"samplesLastChain(3)"；

put@1"samplesThin(1)"；

put@1"samplesStats('*')"；

put@1"samplesDensity('*')"；

put@1"samplesAutoC('*')"；

put@1"samplesTrace('*')"；

put@1"samplesHistory('*')"；

put@1"samplesQuantiles('*')"；

put@1"samplesBgr('*')"；

put@1"samplesCoda('lam'，'gam')"；

put@1"dicset('*')"；

put@1"dicstats('*')"；

put@1"modelSaveLog('&dirout.log.&suffix')"；

put@1"modelSaveLog('&dirout.log.txt')"；

put@1"modelQuit(Y)"；

proc printto；run；

**********************;

/*通過SAS調用Openbugs進行貝葉斯分析*/

**********************;

options xmin noxwait；

x cd C：Program Files(x86)OpenBUGSOpenBUGS323；

x OpenBUGS.exe /PAR &batchf..txt；

**********************;

/*從貝葉斯分析軟件.ODC導出到文本文件，導出到電子表格*/

/*i=0，j=9+number of paremetrics+6+1，rep>j*/

**********************;

data _null_；

retain i 0 j 200；

infile "&exceldir.log.txt" expandtabs truncover；

length text $200；

input text $ 1-200；

if scan(text，1)="node" then i=_n_；

call symput("i"，i+1)；

if upcase(scan(text，1))="SAVE" then j=_n_；

call symput("j"，j-1)；

run；

data results；

infile "&exceldir.log.txt" firstobs=&i obs=&j expandtabs truncover；

length var $20；

input var $ mean se MCMCer lowCI median uppCI start samp；

rep=&count；

run；

proc append base=&dsout data=results；

run；

data bayesian；/*drop redundancy case and variables*/

set results；

if var="model"|var="data"|var="initial"|var="100"|var="Node"|var="mean"|var="Posterior"

|var="Auto-correlation"|var="Dynamic"|var="History"|var="Running"|var="CODA" then delete；

run；

PROC EXPORT DATA= WORK.Bayesian

OUTFILE="&exceldir&modelf.&datafirst..csv"

DBMS=CSV REPLACE；

PUTNAMES=YES；

RUN；

%MENDOPENBUGSIO_model_28u_att1；

**********************；

**SAS宏程序3

/*批處理同參數，同MCMC GIBBS 抽樣數模型*/

**********************;

%macro

openbugswuxi(firstdir=，numit=，n=，dirmodelin=，p=，SMPLN=，modelf=，init1=，init2=，init3=，count=，dirout=)；

%OPENBUGSIO_model_28u_att1(…)；

…

%mendopenbugswuxi；

運行該宏程序前，需要提供貝葉斯統計軟件OpenBUGS分析用的文本文件，包括一個模型文件、數據文件、初始鏈文件。此處的模型文件為model_28u_att1。數據文件為兩個opb1obugfirst和opb1obuglast，為SAS抽樣產生，分別為數據的第一行構成文件opb1obugfirst，其余部分構成文件opb1obuglast。如果不通過SAS產生數據文件，可以適當修改SAS宏程序，直接調用數據文件。初始鏈文件包括3個初始鏈，init1、init2和init3。在此處的宏程序中，應用的文件類型為文本文件，如果是其他類型文件，可作適當修改。

**********************;

/*嵌套之前3宏程序，后臺運行宏*/

**********************;

%openbugswuxi(firstdir=F：yanjiuqinchronic200dirichlet，numit=1000，n=650，dirmodelin=28u，p=28，SMPLN=1，modelf=model_28u_att1，init1=init1，init2=init2，init3=init3，count=200)；

run；

討 論

OpenBUGS是一個專用貝葉斯分析軟件，可應用于簡單或復雜分布的貝葉斯參數估計，特別是對非顯式表達模型的參數估計有極強的優勢[3]。其經歷了多個版本變更，近年來發展很快，仍存在一些局限，如數據管理欠完美。本文討論了利用SAS進行數據文件的管理，然后用SAS宏調用數據文件后臺運行，節省時間，避免人力資源浪費。另外近年發展的免費開源軟件R也提供了一些MCMC貝葉斯分析的程序[17-18]。鑒于SAS軟件的強大數據管理功能，如果有條件建議采用SAS；可是SAS軟件為商業軟件，并非開源軟件，如果從免費開源角度，可以采用R軟件。

本文先用簡單例子簡要介紹了貝葉斯分析的基本步驟和要素，以便理解OpenBUGS的貝葉斯分析基本過程。在此基礎上，為了充分發揮SAS宏數據管理和后臺批處理優勢，以含dirichlet先驗分布的結構方程模型編寫了復雜SAS宏程序實現。dirichlet先驗主要是用于多分類變量先驗，具有自動聚類功能，是近些年無信息先驗的一個重要應用，本模型中應用它是考慮到先驗信息不充分。由于本文的研究重點是探討SAS軟件用于復雜結構方程模型的貝葉斯分析實現，所以對dirichlet先驗不作較多介紹，可參考[19-21]。模擬研究時，當研究者在復雜抽樣，復雜模型貝葉斯應用分析時，可適當修正該宏程序后使用。

對于復雜模型貝葉斯分析，如果需要用SAS編寫宏程序進行后臺分析時，編寫程序要點在于：一是根據OpenBUGS軟件分析步驟，如何產生三個相關文件，即模型文件，數據文件和初始鏈文件；二是用SAS產生OpenBUGS軟件需要的批處理文件“batch script”；三是用SAS調用OpenBUGS軟件執行，編寫程序輸出結果到電子表格中。關于文中的第三個宏程序，用于研究者批量調用宏時，可適當修改使用。