秦正積 嚴(yán)曉玲 沈 毅 肖 靜 何 書 任文龍

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

但是WINBUGS或OpenBUGS軟件數(shù)據(jù)管理功能比較局限。大型數(shù)據(jù)庫管理，譬如數(shù)據(jù)挖掘常常需要從海量數(shù)據(jù)中抽樣獲取，SAS數(shù)據(jù)管理功能強大，數(shù)據(jù)管理在數(shù)據(jù)挖掘中有較大的優(yōu)勢，如果將兩個軟件合理有序的聯(lián)合應(yīng)用，便于計算機后臺工作，可節(jié)省人機對話的時間和人力。本文探討了應(yīng)用OpenBUGS進(jìn)行貝葉斯結(jié)構(gòu)方程模型分析的SAS程序?qū)崿F(xiàn)。

研究設(shè)計與方法

研究從簡單實例介紹OpenBUGS貝葉斯分析原理與實現(xiàn)過程，然后構(gòu)建復(fù)雜結(jié)構(gòu)方程模型。應(yīng)用SAS 9.4和OpenBUGS軟件的GIBBS抽樣方法進(jìn)行貝葉斯分析，并編寫有關(guān)的SAS宏程序從實例展示其后臺實現(xiàn)。

貝葉斯分析原理與實例

1.貝葉斯分析原理與OpenBUGS分析程序基本結(jié)構(gòu)

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

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

為了簡明的闡述軟件執(zhí)行的過程，此處利用軟件的例子進(jìn)行說明。

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

表1 不同年齡小鼠體重測量值統(tǒng)計表

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

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

αi～Normal(αc,τα)

βi～Normal(βc,τβ)

由此，其模型設(shè)計的程序如下：

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)。

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

表2 貝葉斯分析參數(shù)估計

OpenBUGS處理結(jié)構(gòu)方程模型貝葉斯分析的SAS宏程序

1.結(jié)構(gòu)方程模型實例

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

數(shù)據(jù)文件是有關(guān)慢性病的大量人群的調(diào)查研究資料，研究從該數(shù)據(jù)庫進(jìn)行了不同比例的抽樣，得到相應(yīng)文件，用于復(fù)雜結(jié)構(gòu)方程模型的貝葉斯統(tǒng)計分析。

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

結(jié)構(gòu)方程模型見圖1。

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

圖1 結(jié)構(gòu)方程模型示意圖

2.SAS宏程序?qū)崿F(xiàn)

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

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

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

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

%macro semmc；

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

%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；

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

/*生成貝葉斯分析需要的數(shù)據(jù)文件*/

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

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；

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

/*產(chǎn)生貝葉斯分析需要的初始鏈文件*/

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

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；

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

/*產(chǎn)生貝葉斯分析需要的批處理文件*/

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

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調(diào)用Openbugs進(jìn)行貝葉斯分析*/

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

options xmin noxwait；

x cd C：Program Files(x86)OpenBUGSOpenBUGS323；

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

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

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

/*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

/*批處理同參數(shù)，同MCMC GIBBS 抽樣數(shù)模型*/

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

%macro

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

%OPENBUGSIO_model_28u_att1(…)；

…

%mendopenbugswuxi；

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

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

/*嵌套之前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是一個專用貝葉斯分析軟件，可應(yīng)用于簡單或復(fù)雜分布的貝葉斯參數(shù)估計，特別是對非顯式表達(dá)模型的參數(shù)估計有極強的優(yōu)勢[3]。其經(jīng)歷了多個版本變更，近年來發(fā)展很快，仍存在一些局限，如數(shù)據(jù)管理欠完美。本文討論了利用SAS進(jìn)行數(shù)據(jù)文件的管理，然后用SAS宏調(diào)用數(shù)據(jù)文件后臺運行，節(jié)省時間，避免人力資源浪費。另外近年發(fā)展的免費開源軟件R也提供了一些MCMC貝葉斯分析的程序[17-18]。鑒于SAS軟件的強大數(shù)據(jù)管理功能，如果有條件建議采用SAS；可是SAS軟件為商業(yè)軟件，并非開源軟件，如果從免費開源角度，可以采用R軟件。

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

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