贝叶斯方法计算后验分布函数常常是复杂的、高维的,甚至是非标准形式的,其积分计算十分困难[1,2]。譬如,在行为学、社会学、心理学和医学中,许多潜变量模型需要用到结构方程模型,而结构方程模型的应用,常常含有高纬复杂的模型,不能显示表达其分布[3],直接估计困难,可用贝叶斯估计方法,特别是GIBBS抽样方法和MCMC估计方法[4,5]。MCMC方法可解决此类问题[6,7],其实现由专业软件WINBUGS[8]或OpenBUGS[9]完成。GIBBS抽样是构建马尔科夫链最流行的MCMC方法[10,11],其计算工作量较大,对计算机软硬件有较高要求。近年来随着计算机软硬件的快速发展,贝叶斯统计方法受到研究者的重视,贝叶斯统计已经成为统计学中很热门的研究课题[12]。

但是WINBUGS或OpenBUGS软件数据管理功能比较局限。大型数据库管理,譬如数据挖掘常常需要从海量数据中抽样获取,SAS数据管理功能强大,数据管理在数据挖掘中有较大的优势,如果将两个软件合理有序的联合应用,便于计算机后台工作,可节省人机对话的时间和人力。本文探讨了应用OpenBUGS进行贝叶斯结构方程模型分析的SAS程序实现。

研究设计与方法

研究从简单实例介绍OpenBUGS贝叶斯分析原理与实现过程,然后构建复杂结构方程模型。应用SAS9.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,τβ)

由此,其模型设计的程序如下:

001

初始文件:

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,2,3,4,5,6,7,8,9,10,11,12]),结构方程参数9个(gam[1,2,3,4,5,6,7,8,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中抽样产生数据文件opb650obugfirst.txt和opb650obuglast.txt等文件*/

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

%macrosemmc;

%doB=1%to2;/*按比例1%,2%进行随机抽样*/

%if&B=1%then%do;%letrate=0.01;%letSMPLN=1;%end;

%if&B=2%then%do;%letrate=0.02;%letSMPLN=2;%end;

procsurveyselectdata=essay2out=opb&SMPLNmethod=sysseed=130032samprate=&rate;

run;

DATAopb&SMPLN.obug1;

SETopb&SMPLN;

DROP_;

format_all_e8.2;

RUN;

datatabvar.opb&SMPLN.obugfirst1tabvar.opb&SMPLN.obuglast;

setopb&SMPLN.obug1;

if_n_=round(65016/10*&rate,0.1)*10thenoutputtabvar.opb&SMPLN.obuglast;

elseoutputtabvar.opb&SMPLN.obugfirst1;

run;

datatabvar.opb&SMPLN.obugfirst;

settabvar.opb&SMPLN.obugfirst1;

com="";

run;

PROCEXPORTDATA=tabvar.opb&SMPLN.obugfirst

OUTFILE="F:＼yanjiu＼qinchronic＼dirichlet＼28u＼origion＼data＼opb&SMPLN.obugfirst.txt"

DBMS=TABREPLACE;

PUTNAMES=NO;

DELIMITER=",";

RUN;

PROCEXPORTDATA=tabvar.opb&SMPLN.obuglast

OUTFILE=“F:＼yanjiu＼qinchronic＼dirichlet＼28u＼origion＼data＼opb&SMPLN.obuglast.txt”

DBMS=TABREPLACE;

PUTNAMES=NO;

DELIMITER=",";

RUN;

%END;/*closethedo-loopforsamplesizecondition.*/

%mendsemmc;

%semmc;

run;

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

**SAS宏程序2

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

%macroOPENBUGSIO_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=);

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

/*读入模型文件,并指定逻辑文件名*/

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

filenamemodelin"&firstdir&dirmodelin＼model&dirmodelin＼&modelf.txt";run;

filenamemmodel"&datadir&modelf..txt";run;

/*procprinttofile=mmodelnew;*/

data_null_;

filemmodelnew;

run;

data_null_;

infilemodelin;

input;

filemmodelmod;put_infile_;run;

procprintto;run;

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

/*生成贝叶斯分析需要的数据文件*/

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

optionspagesize=32767nodatenocenternonumberformdlim="mprintsymbolgenmlogic;

title1";

filenamedat1"&datadir&data1..txt";run;

filenamedatafrst"&datadir&datafirst..txt";run;

filenamedatalast"&datadir&datalast..txt";run;

filenamedat2"&dir&data2..txt";run;

data_null_;

filedat1new;

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

data_null_;filedat2new;put@1"),.Dim=c(&n,&p)))";run;

data_null_;infiledatafrst;input;filedat1mod;put_infile_;run;

data_null_;infiledatalast;input;filedat1mod;put_infile_;run;

data_null_;infiledat2;input;filedat1mod;put_infile_;run;

procprintto;run;

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

/*产生贝叶斯分析需要的初始链文件*/

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

filenameinit1"&datadir&init1..txt";

filenameinit2"&dir&init2..txt";

filenameinit3"&dir&init3..txt";

filenameinit"&datadir&dirmodelin&init..txt";

data_null_;

fileinitnew;

run;

data_null_;infileinit1;input;fileinitmod;put_infile_;run;

data_null_;infileinit2;input;fileinitmod;put_infile_;run;

data_null_;infileinit3;input;fileinitmod;put_infile_;run;

procprintto;run;

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

/*产生贝叶斯分析需要的批处理文件*/

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

optionspagesize=32767nodatenocenternonumberformdlim="mprintsymbolgenmlogic;

title1";

filenamefileout2"C:＼ProgramFiles(x86)＼OpenBUGS＼OpenBUGS323＼&batchf..txt";

data_null_;

filefileout2;

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)";

procprintto;run;

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

/*通过SAS调用Openbugs进行贝叶斯分析*/

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

optionsxminnoxwait;

xcdC:＼ProgramFiles(x86)＼OpenBUGS＼OpenBUGS323;

xOpenBUGS.exe/PAR&batchf..txt;

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

/*从贝叶斯分析软件.ODC导出到文本文件,导出到电子表格*/

/*i=0,j=9+numberofparemetrics+6+1,rep>j*/

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

data_null_;

retaini0j200;

infile"&exceldir.log.txt"expandtabstruncover;

lengthtext$200;

inputtext$1-200;

ifscan(text,1)="node"theni=_n_;

callsymput("i",i+1);

ifupcase(scan(text,1))="SAVE"thenj=_n_;

callsymput("j",j-1);

run;

dataresults;

infile"&exceldir.log.txt"firstobs=&iobs=&jexpandtabstruncover;

lengthvar$20;

inputvar$meanseMCMCerlowCImedianuppCIstartsamp;

rep=&count;

run;

procappendbase=&dsoutdata=results;

run;

databayesian;/*dropredundancycaseandvariables*/

setresults;

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

|var="Auto-correlation"|var="Dynamic"|var="History"|var="Running"|var="CODA"thendelete;

run;

PROCEXPORTDATA=WORK.Bayesian

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

DBMS=CSVREPLACE;

PUTNAMES=YES;

RUN;

%MENDOPENBUGSIO_model_28u_att1;

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

**SAS宏程序3

/*批处理同参数,同MCMCGIBBS抽样数模型*/

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

%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:＼yanjiu＼qinchronic200＼dirichlet\,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,20,21]。模拟研究时,当研究者在复杂抽样,复杂模型贝叶斯应用分析时,可适当修正该宏程序后使用。

对于复杂模型贝叶斯分析,如果需要用SAS编写宏程序进行后台分析时,编写程序要点在于:一是根据OpenBUGS软件分析步骤,如何产生三个相关文件,即模型文件,数据文件和初始链文件;二是用SAS产生OpenBUGS软件需要的批处理文件“batchscript”;三是用SAS调用OpenBUGS软件执行,编写程序输出结果到电子表格中。关于文中的第三个宏程序,用于研究者批量调用宏时,可适当修改使用。

参考文献：

[22]董大均,张尔强,等.SAS语言指导.辽宁科学技术出版社,1991:359.

[23]胡良平,现代统计学与SAS应用,军事医学科学出版社,2000.

秦正积,严晓玲,沈毅,肖静,何书,任文龙.OpenBUGS处理结构方程模型贝叶斯分析在SAS宏程序中的实现[J].中国卫生统计,2020,37(03):475-480.

基金:江苏省卫生和计划委员会项目(Y2018099).