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软件执行，编写程序输出结果到电子表格中。关于文中的第三个宏程序，用于研究者批量调用宏时，可适当修改使用。