R语言贝叶斯非参数模型:密度估计、非参数化随机效应meta分析心肌梗死数据

最近,我们使用贝叶斯非参数(BNP)混合模型进行马尔科夫链蒙特卡洛(MCMC)推断。

由Kaizong Ye,Weilong Zhang撰写

在这篇文章中,我们通过展示如何使用具有不同内核的非参数混合模型进行密度估计。在后面的文章中,我们将采用参数化的广义线性混合模型,并展示如何切换到非参数化的随机效应表示,避免了正态分布的随机效应假设。


使用Dirichlet Process Mixture模型进行基本密度估计

提供了通过Dirichlet过程混合(DPM)模型进行非参数密度估计的机制(Ferguson, 1974; Lo, 1984; Escobar, 1994; Escobar and West, 1995)。对于一个独立和相同分布的样本 ,该模型的形式为


这个模型实现是灵活的,运行任意核的混合。, 可以是共轭的,也可以是不共轭的(也是任意的)基度量 . 在共轭核/基数测量对的情况下,能够检测共轭的存在,并利用它来提高采样器的性能。

×

简单说贝叶斯非参数模型就是参数是无限维的贝叶斯模型。

PS:贝叶斯非参和非参(比如K近邻、决策树之类的)有共同点,就是模型可以自动适应数据。不同的是贝叶斯非参是在贝叶斯的框架下提出的。

模型是用来描述数据的。概率模型就是从概率的角度描述数据:对于给定的数据集[公式],我们认为这个数据集里面的元素是从某个概率分布[公式]独立采样得到的,即[公式]

模型估计问题就是给定[公式]找到这样一个[公式]。之后我们会找一个参数控制的模型的集合[公式],把在函数空间中寻找最佳[公式]的问题转换成某参数空间[公式]中寻找最佳参数[公式]的问题。例如我们可以定义某目标函数[公式],例如似然函数[公式],这里[公式]为分布[公式]的概率密度函数,然后在[公式]上优化[公式]就可以得到最佳参数[公式]和对应的模型[公式]

上面这个方法叫频率方法。贝叶斯是另一回事儿,不过[公式]还是有的,但是我的目标并不是找一个特定的[公式]作为最优解,而是在[公式]上定义一个概率分布[公式]作为[公式]的先验分布,然后再根据数据[公式]计算[公式]的后验分布,即[公式],这里[公式][公式]对应的密度函数。

举个例子,给你个硬币,让你估计正面的概率(伯努利分布的参数),你扔了100次,有60次朝上,频率方法估计的结果就是正面向上的概率为0.6,贝叶斯方法的结果可能是这么个样子:正面向上概率为0.6的概率是0.9,概率为0.4的概率为0.1

而贝叶斯参数和非参数的区别主要就在这个[公式]上面。简单说如果[公式]是有限维的,比如[公式],那这个模型就叫参数模型,如果[公式]是无限维的就叫非参数模型。非参数模型中的[公式]一般要借助随机过程来描述。



为了说明这些能力,我们考虑对R中提供的Faithful火山数据集的喷发间隔时间的概率密度函数进行估计。

观测值  对应于数据框架的第二列,而 .

使用CRP表示法拟合高斯_location-scale_ 分布混合分布

模型说明

我们首先考虑用混合正态分布的_location-scale_Dirichlet过程s来拟合转换后的数据


视频

R语言中RStan贝叶斯层次模型分析示例

探索见解

其中 对应的是正态-逆伽马分布。这个模型可以解释为提供一个贝叶斯版本的核密度估计 用于使用高斯核和自适应带宽。在数据的原始尺度上,这可以转化为一个自适应的对数高斯核密度估计。

引入辅助变量,表明混合的哪个成分产生了每个观测值,并对随机量进行积分,我们得到模型的CRP表示(Blackwell and MacQueen, 1973)。

其中

是向量中唯一值的数量,是第个唯一值在中出现的次数。这个说明清楚地表明,每个观测值都属于最多正态分布聚类中的任何一个,并且CRP分布与分区结构的先验分布相对应。

 

这个模型的说明是这样的

请注意,在模型代码中,参数向量muTilde和s2Tilde的长度被设置为.我们这样做是因为目前的实现要求提前设置参数向量的长度,并且不允许它们的数量在迭代之间变化。因此,如果我们要确保算法总是按预期执行,我们需要在最坏的情况下工作,即有多少个成分就有多少个观测值的情况。但它的效率也有点低,无论是在内存需求方面(当 规模大时,需要维护大量未占用的成分)还是在计算负担方面(每次迭代都需要更新大量不需要后验推理的参数)。当我们在下面使用伽马分布的混合时,我们将展示一个能提高效率的计算捷径。

还需要注意的是,的值控制着我们先验预期的成分数量,的值越大,对应于数据占据的成分数量越多。因此,通过指定一个先验值,我们为模型增加了灵活性。对Gamma先验的特殊选择允许使用数据增强方案从相应的全条件分布中有效取样。也可以选择其他的先验,在这种情况下,这个参数的默认采样是一个自适应的随机游走Metropolis-Hastings算法。

运行MCMC算法

下面的代码设置了数据和常数,初始化了参数,定义了模型对象,并建立和运行了MCMC算法。默认采样器是一个折叠的吉布斯采样器(Neal, 2000)。

我们可以从参数的后验分布中提取样本,并创建痕迹图、直方图和任何其他感兴趣的总结。例如,对于参数,我们有。


WINBUGS对随机波动率模型进行贝叶斯估计与比较

阅读文章



随时关注您喜欢的主题


在这个模型下,对于一个新的观察,后验预测分布是最佳密度估计(在平方误差损失下)。这个估计的样本可以很容易地从我们的MCMC产生的样本中计算出来。

然而,回顾一下,这是对等待时间的对数的密度估计。为了获得原始尺度上的密度,我们需要对内核进行适当的转换。

无论是哪种情况,都有明显的证据表明,数据中的等待时间有两个组成部分。

生成混合分布的样本

虽然混合分布的线性函数的后验分布的样本(比如上面的预测分布)可以直接从折叠采样器的实现中计算出来,但是对于非线性函数的推断需要我们首先从混合分布中生成样本。我们可以从随机度量中获得后验样本。需要注意的是,为了从 ,得到后验样本,我们需要监控所有参与其计算的随机变量,即成员变量xi,聚类参数muTilde和s2Tilde,以及浓度参数alpha。

下面的代码从随机测量中生成后验样本。cMCMC对象包括模型和参数的后验样本。

函数估计了一个截断水平,即truncG。后验样本是一个带列的矩阵,其中参数分布向量的维度(在本例中为)。

下面的代码使用随机测量的后验样本来计算的后验样本。请注意,这些样本是基于转换后的模型计算的,大于70的值对应于上述定义的网格上大于0.035的值。

使用CRP表示法拟合伽马混合分布

不限于在DPM模型中使用高斯核。就Old Faithful数据而言,除了我们在上一节中介绍的对数尺度上的高斯核的混合分布外,还有一种选择是数据原始尺度上的伽马混合分布。

模型

在这种情况下,模型的形式为

其中对应于两个独立Gamma分布的乘积。下面的代码提供了该模型。

请注意,在这种情况下,向量beta和lambda的长度为 。这样做是为了减少与采样算法有关的计算和存储负担。你可以把这种方法看作是对过程的截断,只不过它可以被认为是*精确的截断。事实上,在CRP表示法下,只要采样器的成分数严格低于采样器每次迭代的参数向量的长度,使用长度短于样本中观察值的参数向量就会生成一个合适的算法。

运行MCMC算法

下面的代码设置了模型数据和常数,初始化了参数,定义了模型对象,并建立和运行了Gamma混合分布的MCMC算法。请注意,在构建MCMC时,会产生一个关于聚类参数数量的警告信息。这是因为betaTilde和lambdaTilde的长度小于。另外,请注意,在执行过程中没有产生错误信息,这表明所需的集群数量未超过50个的上限。

在这种情况下,我们使用参数的后验样本来构建一个轨迹图,并估计的后验分布。

从混合分布中生成样本

和以前一样,我们从后验分布中获得样本。

我们使用这些样本来创建一个数据密度的估计值,以及一个95%置信带。

我们再次看到,数据的密度是双峰的,看起来与我们之前得到的数据非常相似。

使用stick-breaking 表示法拟合伽马DP混合分布

模型

Dirichlet过程混合物的另一种表示方法是使用随机分布的stick-breaking表示(Sethuraman, 1994)。

引入辅助变量,表明哪个成分产生了每个观测值,上一节讨论的Gamma密度的混合物的相应模型的形式为

其中 是两个独立Gamma分布的乘积。

. 下面的代码提供了该模型说明。

注意,截断水平已被设置为Trunc值,该值将在函数的常数参数中定义。

运行MCMC算法

下面的代码设置了模型数据和常量,初始化了参数,定义了模型对象,并建立和运行了Gamma混合分布的MCMC算法。当使用stick-breaking表示时,会指定一个分块Gibbs抽样器(Ishwaran, 2001; Ishwaran and James, 2002)。

使用stick-breaking近似法会自动提供随机分布的近似值,即 。下面的代码使用来自样本对象的后验样本计算后验样本,并从中计算出数据的密度估计。

正如预期的那样,这个估计值看起来与我们通过CRP表示的过程获得的估计值相同。

贝叶斯非参数化:非参数化随机效应

我们将采用一个参数化的广义线性混合模型,并展示如何切换到非参数化的随机效应表示,避免了正态分布的随机效应假设。

心肌梗死(MIs)的参数化meta分析

我们将在对以前非常流行的糖尿病药物 “Avandia “的副作用进行meta分析的背景下,说明使用非参数混合模型对随机效应分布进行建模。我们分析的数据在引起对这种药物的安全性的严重质疑方面发挥了作用。问题是使用”Avandia “是否会增加心肌梗死(心脏病发作)的风险。,每项研究都有治疗和对照组。

模型的制定

我们首先进行基于标准的广义线性混合模型(GLMM)的meta分析。向量n和x分别包含对照组的患者总数和每项研究中对照组的心肌梗死患者人数。同样,向量m和y包含接受药物的病人的类似信息。该模型的形式为

其中,随机效应、遵循共同的正态分布被合理地赋予非信息性先验。参数量化了对照组和治疗组之间的风险差异,而参数则量化了研究的具体变化。

这个模型可以用以下代码指定。

运行MCMC

让我们来运行一个基本的MCMC。

结果表明,对照组和治疗组之间存在着整体的风险差异。但是正态性假设呢?我们的结论对该假设是否稳健?也许随机效应的分布是偏斜的。

用于meta分析的基于DP的随机效应模型

模型

现在,我们对使用非参数分布。更具体地说,我们假设每个都是由位置尺度的正态混合分布产生的。

这种模型引起了随机效应之间的聚类。与密度估计问题的情况一样,DP先验允许数据决定分量的数量,从最少的一个分量(即简化为参数模型)到最多的分量,即每个观测值有一个分量。如果数据支持这种行为,这允许随机效应的分布是多模态的,大大增加了其灵活性。这个模型可以用以下代码指定。

运行MCMC

以下代码对模型进行了编译,并对模型运行了一个压缩Gibbs抽样

主要推论似乎对原始的参数化假设很稳健。这可能是由于没有太多证据表明随机效应分布中缺乏正态性。

参考文献

Blackwell, D. and MacQueen, J. 1973. Ferguson distributions via Polya urn schemes. The Annals of Statistics 1:353-355.

Ferguson, T.S. 1974. Prior distribution on the spaces of probability measures. Annals of Statistics 2:615-629.

Lo, A.Y. 1984. On a class of Bayesian nonparametric estimates I: Density estimates. The Annals of Statistics 12:351-357.


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498