霍乱传染数据贝叶斯ZIP分析的后验预测评价*

王 婷 曾 平，2△

霍乱传染数据贝叶斯ZIP分析的后验预测评价*

王 婷1曾 平1，2△

1926年英国流行病学家麦克德里克(McKendrick)应用数学模型研究了1906年印度孟买一个村庄流行性霍乱传染的数据[1]。为了描述霍乱传播过程，麦克德里克首先考虑了简单的Poisson分布，在发现Poisson分布拟合效果很差后，麦克德里克又考虑了零截尾Poisson分布(zero-truncated Poisson)，Irwin采用EM算法改进了麦克德里克的方法[2]。但是无论哪种分布，麦克德里克都没能很好处理霍乱数据中存在的大量没有患者的家庭，即霍乱感染者个数为0。麦克德里克也注意到这个问题，并给出了一个十分合理的解释：霍乱病菌经水源传播，村庄中有许多水井，也许只有某些水井被病菌污染，而一些水井没有被污染，是干净清洁的。因此，饮用被病菌污染井水的家庭中有人患病，而饮用干净井水的村民自然不会患病，这就揭示了为什么数据中存在大量没有霍乱患者的家庭。根据麦克德里克的理由，Meng认为可采用了零膨胀Poisson模型(zero inflated Poisson，ZIP)来分析霍乱传染数据[3]。一个很自然的问题是，ZIP模型对该数据的拟合优度如何？如何做出评价？为此，本文主要研究基于后验预测分布的贝叶斯ZIP模型评价，并通过该数据展示相应的方法。

贝叶斯ZIP模型

霍乱传染数据描述了孟买一个村庄223户家庭流行性霍乱传染情况[2]，见表1。表中x表示被感染家庭中霍乱患者人数，nx表示患者人数为x的家庭个数。数据显示大约75%的家庭没有霍乱患者，大约14%的家庭有1个患者，7%的家庭有2个患者，患者人数超过3个的家庭数不足4%。霍乱传染数据的ZIP模型可表示如下[4-6]：

(1)

p表示家庭中成员不可能感染霍乱的概率，I为指示函数，当x=0取值为1，否则为0。ZIP分布实际是退化到0的分布(概率为p)和Poisson分布(概率为1-p)的混合分布。

(2)

∝表示两边相差一个常数因子。采用基于MCMC(Markov Chain Monte Carlo)的随机游走Metropolis模拟算法[7]，通过构造平稳分布为后验分布的Markov链来得到参数的后验样本(pt,μt)，Monte Carlo样本量为105。图1-2给出了参数p和μ的直方图和等高线图[8]，表2为后验统计量。随机游走Metropolis模拟在SAS9.2的MCMC程序下完成[7]，后续统计分析在R2.11.1下完成[8-9]。

*:数据来源于《现代数学手册-随机数学卷》。贝叶斯ZIP模型拟合时选择参数的中位数，分别为p=0.5930和μ=0.9470。

后验预测模型评价

模型后验评价是建立在后验预测分布基础上的。后验预测分布(posterior predictive distribution)是指在观察到数据X之后未来可能观察到数据X*的分布[10-11]。称后验是因为数据X*的分布是建立在已观察到数据X之上的条件分布，称预测是因为该分布是针对尚没有观察到的数据而言。X*的后验预测分布可表示为：

f(X*|X)=∫f(X*|X,p,μ)f(p,μ|X)dpdμ, =∫f(X*|p,u)f(p,μ|X)dpdμ,(3)

上式中第二个等式成立是基于f(X*|X,p,u)=f(X*|p,μ)的事实，也即是X*的后验预测分布是参数p,μ的条件分布，X只通过参数p,μ来影响X*的分布，和已经观察到的数据X无关，在已知p,μ时X*和X条件独立。从后验预测分布中产生的数据又称为重复数据(replicated data)，本文中X指223户家庭霍乱发病人数，X*指预测霍乱患者人数。后验预测分布量化了在已经观察到数据X后在未来能够再次观察到X的可能性，因此不但可以用来预测未知数据，还能进行模型检查评价。

从后验预测分布随机抽取样本X*，然后将这些样本和观察到的数据X比较，如果模型对数据拟合得好，那么重复数据应该看上去和数据X一致，两者之间的明显差别表明选择模型的偏离。可以通过定义一个检验量T来量化这种偏离程度[11-12]，T(X*)表示应用预测数据得到的检验量，为一个随机变量，T(X)表示应用观察数据得到的检验量，为一个确定的数值。在霍乱传染数据的ZIP模型中，合适的检验量是数据中0的家庭数T(X*)=n0。后验检验的策略是，如果ZIP模型拟合良好，那么从预测分布f(X*|X)产生的0的家庭数T(X*)=n0和T(X)=n0应该接近，T(X*)和T(X)之间系统的差别提示ZIP模型拟合不足。本文中T(X)=168，图3给出了霍乱传染数据和9个预测数据的频数分布图，其中预测数据由公式(3)随机产生。图3显示预测数据分布和实际观察数据图形相近。

霍乱传染数据ZIP模型的后验模型检验具体过程如下[11-12]：

(1)从后验分布f(p,μ|X)抽取(pt,μt),t=1，…,104,这一过程已经由随机游走Metropolis模拟完成；

(3)计算后验预测样本X*中患者为0的家庭数T(X*)；

(4)对所有104组参数(pt,μt)重复(1)～(3)步，则产生T(X*)t,t=1,…,104；

如果ZIP能拟合霍乱数据，则p不应该出现极端值，例如很大(大于0.95)或很小(小于0.05)。T(X*)的直方图见图4，参考线为168的位置。T(X*)≥T(X)和T(X*)≤T(X)的比例分别是0.507和0.537，两者之和不为1是因为不等式中都存在等号的原因。这表明如果霍乱患者传染数据服从ZIP分布，那么观察到高达75%的家庭没有患者的可能很大，这种情况不大可能是偶然发生的。最后的结果显示T(X*)在135～196之间，P25和P75分别为162和174，均数和中位数分别为167.5和168，标准差为9.04。特殊地，如果用p=0.593和μ=0.947预测霍乱发病人数，0的预测值为167，概率为75.01%，和实际168十分一致，除此之外，其他预测值和实际观察到的发病人数也很接近。这些数值表示ZIP模型很好地拟合了实际数据，能够有效地处理数据中存在的大量0记录的情况。

如果选择Poisson模型，先验选择共轭gamma分布，指定α=β=10-3，贝叶斯Poisson模型T(X*)的分布如图5。T(X*)介于112～184之间，P25和P75分别为145和158，均数和中位数分别为151.7和152，标准差为9.38。T(X*)≥T(X)和T(X*)≤T(X)的比例分别是0.044和0.966，这表示从Poisson模型中不大可能观察到高达75%的零记录，Poisson模型对数据拟合不够，至少不能处理数据中零过多的问题。

讨 论

Poisson分布和其他模型对霍乱传染数据拟合效果差的一个主要原因在于不能很好地处理那些没有患者的家庭。由于霍乱病毒只污染了村里部分水井，那么饮用健康井水的村民无疑是不会患病的，即便是饮用霍乱污染的井水，村民也不一定表现出明显的患病症状，正是这两种原因使得很多家庭没有霍乱患者，其中前者是Poisson分布不能解释的0的来源。我们采用ZIP模型对这个混合总体的不同人群分别建立模型，给前者指定为退化到0点的分布，后者指定为Poisson分布，ZIP模型就是这两个分布的混合，准确地捕捉到了数据结构所呈现的信息。贝叶斯ZIP模型结果显示，预测的家庭患病人数和实际观察数据十分接近，尤其是准确地预测了没有患者的家庭数。

为了评价麦克德里克问题ZIP模型的拟合优度，

我们采用后验预测分布对ZIP模型做出评价。后验预测检查需要定义一个合适的检验量T，T不但可以是数据的函数T(X)，如本文中所选择的检验量，这和经典统计假设检验统计量类似，还可以是数据和参数的函数T(X,θ)，后验预测分布重复数据的T(X*)和实际数据的T(X)之间的差异反映了模型和数据之间的吻合程度。后验预测检查的原理在于，如果实际数据违背了模型的重要假设，那么T(X*)相对T(X)就会表现出具有统计学意义的差别来。如果得到的后验预测p值接近0或1，暗示实际数据具有极端的检验量T，而选择的模型很可能是不合适的。在后验评价中采用重复数据的0的比例作为检验量，这也是选择ZIP模型的主要理由。因此，如果对ZIP模型是恰当的，那么由后验预测分布所产生的0就应该和实际观察到的0的比例一致。后验检查显示，ZIP模型能够很好地拟合实际数据，而Poisson模型不能很好地处理数据中存在的0过多问题。除此之外，差别检验量还可选择χ2值或变异系数。

[1]Http://en.wikipedia.org/wiki/Anderson_Gray_McKendrick.

[2]徐利治.现代数学手册-随机数学卷.武汉:华中科技大学出版社,1999：357-390.

[3]Meng XL.The EM Algorithm and Medical Studies:A Historical Link.Statistical Methods s in Medical Research,1997,6(1):3-23.

[4]曾平,刘桂芬,曹红艳.零膨胀模型在心肌缺血节段数影响因素研究中的应用.中国卫生统计,2008,25(5):464-466.

[5]曾平.零过多计数资料回归模型及其医学应用.硕士学位论文,太原:山西医科大学,2009.

[6]Lambert D.Zero-inflated Poisson Regression with an Application to Defects in Manufacturing.Technimetrics,1992,34(1):1-14.

[7]SAS Institute Inc.Preliminary Capabilities for Bayesian Analysis in SAS/STAT Software.Cary,NC,USA,2006.

[8]Albert J.Bayesian Computation with R.Second Edition.New York:Springer,2009.

[9]R Development Core Team.R:A language and environment for statistical computing.R Foundation for Statistical Computing,Vienna,Austria,2014.URL http://www.R-project.org.

[10]Ntzoufras I.Bayesian Modeling Using WinBUGS.New York:John Wiley & Sons,2009.

[11]Gelman A,Carlin JB,Stern HS,et al.Bayesian Data Analysis.Second Edition.London:Chapman & Hall,2004.

[12]Gelman A,Meng XL,Stern HS.Posterior Predictive Assessment of Model Fitness via Realized Discrepancies.Statistica Sinica,1996,6(4):733-807.

[13]Meng XL.Posterior predictive p-values.Annals of Statistics,1994,22(3):1142-1160.

(责任编辑：邓 妍)

*国家自然科学基金项目(81402765)；国家统计局全国统计科学研究项目(2014464)；江苏省教育厅高校哲学社会科学研究基金项目(2013SJB790059,2013SJD790032)

1.徐州医学院公共卫生学院流行病与卫生统计学教研室(221004)

2.南京医科大学公共卫生学院流行病与卫生统计学系

△ 通信作者：曾平，E-mail：zpstat@xzmc.edu.cn