关于数据挖掘:R语言贝叶斯非参数模型密度估计非参数化随机效应meta分析心肌梗死数据附代码数据

原文链接：http://tecdat.cn/?p=23785

原文出处：拓端数据部落公众号

最近咱们被客户要求撰写对于贝叶斯非参数模型的钻研报告，包含一些图形和统计输入。

概述

最近，咱们应用贝叶斯非参数（BNP）混合模型进行马尔科夫链蒙特卡洛（MCMC）推断。

在这篇文章中，咱们通过展现如何应用具备不同内核的非参数混合模型进行密度估计。在前面的文章中，咱们将采纳参数化的狭义线性混合模型，并展现如何切换到非参数化的随机效应示意，防止了正态分布的随机效应假如。

应用 Dirichlet Process Mixture 模型进行根本密度估计

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

这个模型实现是灵便的，运行任意核的混合。, 能够是共轭的，也能够是不共轭的（也是任意的）基度量 . 在共轭核 / 基数测量对的状况下，可能检测共轭的存在，并利用它来进步采样器的性能。

为了阐明这些能力，咱们思考对 R 中提供的 Faithful 火山数据集的喷发间隔时间的概率密度函数进行预计。

data(faithful)

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

应用 CRP 表示法拟合高斯 location-scale 散布混合散布

模型阐明

咱们首先思考用混合正态分布的 location-scaleDirichlet 过程 s 来拟合转换后的数据

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

引入辅助变量，表明混合的哪个成分产生了每个观测值，并对随机量进行积分，咱们失去模型的 CRP 示意（Blackwell and MacQueen, 1973）。

其中

是向量中惟一值的数量，是第个惟一值在中呈现的次数。这个阐明分明地表明，每个观测值都属于最多正态分布聚类中的任何一个，并且 CRP 散布与分区构造的先验散布绝对应。

这个模型的阐明是这样的

    y[i] ~ dnorm(mu[i], var = s2[i])
    mu[i] <- muTilde[xi[i]]
    s2[i] <- s2Tilde[xi[i]]
  xi[1:n] ~ dCRP(alpha, size = n)
    muTilde[i] ~ dnorm(0, var = s2Tilde[i])
    s2Tilde[i] ~ dinvgamma(2, 1)
  alpha ~ dgamma(1, 1)

请留神，在模型代码中，参数向量 muTilde 和 s2Tilde 的长度被设置为. 咱们这样做是因为目前的实现要求提前设置参数向量的长度，并且不容许它们的数量在迭代之间变动。因而，如果咱们要确保算法总是按预期执行，咱们须要在最坏的状况下工作，即有多少个成分就有多少个观测值的状况。但它的效率也有点低，无论是在内存需要方面（当 规模大时，须要保护大量未占用的成分）还是在计算累赘方面（每次迭代都须要更新大量不须要后验推理的参数）。当咱们在上面应用伽马散布的混合时，咱们将展现一个能提高效率的计算捷径。

还须要留神的是，的值管制着咱们先验预期的成分数量，的值越大，对应于数据占据的成分数量越多。因而，通过指定一个先验值，咱们为模型减少了灵活性。对 Gamma 先验的非凡抉择容许应用数据加强计划从相应的全条件散布中无效取样。也能够抉择其余的先验，在这种状况下，这个参数的默认采样是一个自适应的随机游走 Metropolis-Hastings 算法。

运行 MCMC 算法

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

# 模型数据
y = standlFaithful
# 模型常量
n = length(standlFaithful))
# 参数初始化
list(xi = sample(1:10, size=n, replace=TRUE),
# 创立和编译模型
Model(code,  data,  inits,  consts)
## 定义模型...
## 建设模型...
## 设置数据和初始值...
## 在模型上运行计算（随后的任何错误报告可能只是反映了模型变量的缺失值）... 
## 查看模型的大小和尺寸......。## 模型构建实现。## 编译实现。#MCMC 的配置、创立和编译
MCMC(conf)
## 编译...... 这可能须要一分钟
## 编译实现。

咱们能够从参数的后验散布中提取样本，并创立痕迹图、直方图和任何其余感兴趣的总结。例如，对于参数，咱们有。\
 

# 参数的痕迹图
ts.plot(samples[ , "alpha"], xlab = "iteration", ylab = expression(alpha))

# 后验直方图
hist(samples[ , "alpha"], xlab = expression(alpha), main = "", ylab ="Frequency")

在这个模型下，对于一个新的察看，后验预测散布是最佳密度估计（在平方误差损失下）。这个预计的样本能够很容易地从咱们的 MCMC 产生的样本中计算出来。

# 参数的后验样本
 samples[, "alpha"]
# 平均值的后验样本
 samples[, grep('muTilde', colnames(samples))] # 聚类平均数的后验样本。# 集群方差的后验样本
samples[, grep('s2Tilde', colnames(samples))] # 聚类成员的后验样本。# 聚类成员关系的后验样本
samples [, grep('xi', colnames(samples))] # 聚类成员的后验样本。hist(y, freq = FALSE,
     xlab = "标准化对数尺度上的等待时间")
## 对标准化对数网格的密度进行点式预计 

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

standlGrid*sd(lFaithful) + mean(lFaithful) # 对数尺度上的网格

hist(faithful$waiting, freq = FALSE
  

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

生成混合散布的样本

尽管混合散布的线性函数的后验散布的样本（比方下面的预测散布）能够间接从折叠采样器的实现中计算出来，然而对于非线性函数的推断须要咱们首先从混合散布中生成样本。咱们能够从随机度量中取得后验样本。须要留神的是，为了从，失去后验样本，咱们须要监控所有参加其计算的随机变量，即成员变量 xi，聚类参数 muTilde 和 s2Tilde，以及浓度参数 alpha。

上面的代码从随机测量中生成后验样本。cMCMC 对象包含模型和参数的后验样本。函数预计了一个截断程度，即 truncG。后验样本是一个带列的矩阵，其中参数散布向量的维度（在本例中为）。

outputG <- getSamplesDPmeasure(cmcmc)

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

     truncG <- outputG$trunc # G 的截断程度



probY70 <- rep(0, nrow(samples))  # P(y.tilde>70) 的后验样本

hist(probY70)

应用 CRP 表示法拟合伽马混合散布

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

模型

在这种状况下，模型的模式为

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

    y[i] ~ dgamma(shape = beta[i], scale = lambda[i])
    beta[i] <- betaTilde[xi[i]]
    lambda[i] <- lambdaTilde[xi[i]]

请留神，在这种状况下，向量 beta 和 lambda 的长度为。这样做是为了缩小与采样算法无关的计算和存储累赘。你能够把这种办法看作是对过程的截断，只不过它能够被认为是 * 准确的截断。事实上，在 CRP 表示法下，只有采样器的成分数严格低于采样器每次迭代的参数向量的长度，应用长度短于样本中察看值的参数向量就会生成一个适合的算法。

运行 MCMC 算法

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

data <- list(y = waiting)
Model(code, data = data)

cModel <- compile

samples <- runMCMC(cmcmc, niter = 7000, nburnin = 2000, setSeed = TRUE)

在这种状况下，咱们应用参数的后验样本来构建一个轨迹图，并预计的后验散布。\
 

# 参数的后验样本的跟踪图
ts.plot(samples[ , 'alpha'], xlab = "iteration", ylab = expression(alpha))

# 参数的后验样本的直方图 
hist(samples[, 'alpha'])

从混合散布中生成样本

和以前一样，咱们从后验散布中取得样本。\
 

outputG <- getSamplesDPmeasure(cmcmc)

咱们应用这些样本来创立一个数据密度的估计值，以及一个 95% 相信带。

for(iter in seq_len)) {density[iter,] <- sapply(grid, function(x)
    sum(weightSamples[iter,] * dgamma)))
}


hist(waiting, freq = FALSE

咱们再次看到，数据的密度是双峰的，看起来与咱们之前失去的数据十分类似。

应用 stick-breaking 表示法拟合伽马 DP 混合散布

模型

Dirichlet 过程混合物的另一种示意办法是应用随机散布的 stick-breaking 示意（Sethuraman, 1994）。

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

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

与. 上面的代码提供了该模型阐明。

      y[i] ~ dgamma(shape = beta[i], scale = lambda[i])
      beta[i] <- betaStar[z[i]]
      lambda[i] <- lambdaStar[z[i]]
      z[i] ~ dcat(w[1:Trunc])
     # stick-breaking 
      v[i] ~ dbeta(1, alpha)
    w[1:Trunc] <- stick_breaking(v[1:(Trunc-1)]) # stick-breaking 权重
      betaStar[i] ~ dgamma(shape = 71, scale = 2)

留神，截断程度已被设置为 Trunc 值，该值将在函数的常数参数中定义。

运行 MCMC 算法

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

data <- list(y = waiting)
consts <-length(waiting)
betaStar = rgamma
              lambdaStar = rgamma
              v = rbeta
              z = sample
              alpha = 1

compile(Model)

MCMC(rModel, c("w", "betaStar", "lambdaStar", 'z', 'alpha'))
comp(mcmc)

MCMC(cmcmc, niter = 24000)

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

  densitySamples[i,] <- sapply(grid, function(x) sum(weightSamples  * dgamma(x, shape ,
                                    scale )))

hist(waiting ylim=c(0,0.05),
 

正如预期的那样，这个估计值看起来与咱们通过 CRP 示意的过程取得的估计值雷同。

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

咱们将采纳一个参数化的狭义线性混合模型，并展现如何切换到非参数化的随机效应示意，防止了正态分布的随机效应假如。

心肌梗死（MIs）的参数化 meta 剖析

咱们将在对以前十分风行的糖尿病药物 “Avandia “ 的副作用进行 meta 剖析的背景下，阐明应用非参数混合模型对随机效应散布进行建模。咱们剖析的数据在引起对这种药物的安全性的重大质疑方面施展了作用。问题是应用 ”Avandia “ 是否会减少心肌梗死（心脏病发生）的危险。，每项钻研都有医治和对照组。

模型的制订

咱们首先进行基于规范的狭义线性混合模型（GLMM）的 meta 剖析。向量 n 和 x 别离蕴含对照组的患者总数和每项钻研中对照组的心肌梗死患者人数。同样，向量 m 和 y 蕴含承受药物的病人的相似信息。该模型的模式为

\
其中，随机效应、遵循独特的正态分布、和被正当地赋予非信息性先验。参数量化了对照组和医治组之间的危险差别，而参数则量化了钻研的具体变动。

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

        y[i] ~ dbin(size = m[i], prob = q[i]) # 药物 MIs
        x[i] ~ dbin(size = n[i], prob = p[i]) # 管制 MIs
        q[i] <- expit(theta + gamma[i]) # 药物的对数指数
        p[i] <- expit(gamma[i]) #对照组对数
        gamma[i] ~ dnorm(mu, var = tau2) # 钻研成果
    theta ~ dflat() # 药物的影响
    # 随机效应超参数
    mu ~ dnorm(0, 10)
    tau2 ~ dinvgamma(2, 1)

运行 MCMC

让咱们来运行一个根本的 MCMC。

MCMC(codeParam,  data, inits,
                      constants, monitors = c("mu", "tau2", "theta", "gamma")

par(mfrow = c(1, 4)

hist(gammaMn)
hist(samples[1000, gammaCols)

结果表明，对照组和医治组之间存在着整体的危险差别。然而正态性假如呢？咱们的论断对该假如是否持重？兴许随机效应的散布是偏斜的。

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

模型

当初，咱们对应用非参数散布。更具体地说，咱们假如每个都是由地位尺度的正态混合散布产生的。

\
这种模型引起了随机效应之间的聚类。与密度估计问题的状况一样，DP 先验容许数据决定重量的数量，从起码的一个重量（即简化为参数模型）到最多的重量，即每个观测值有一个重量。如果数据反对这种行为，这容许随机效应的散布是多模态的，大大增加了其灵活性。这个模型能够用以下代码指定。

        y[i] ~ dbin(size = m[i], prob = q[i]) # 药物 MIs
        x[i] ~ dbin(size = n[i], prob = p[i]) # MIs
        q[i] <- expit(theta + gamma[i]) # 药物的对数指数
        p[i] <- expit(gamma[i]) # 对照组对数值
        gamma[i] ~ dnorm(mu[i], var = tau2[i]) # 来自混合物的随机效应。mu[i]<- muTilde[xi[i]]                 # 来自聚类的随机效应的平均值 xi[i]
        tau2[i] <- tau2Tilde[xi[i]]           # 来自群组 xi[i] 的随机效应变量
    # 从根底测量中提取混合成分参数
        muTilde[i] ~ dnorm(mu0, var = var0)
        tau2Tilde[i] ~ dinvgamma(a0, b0)
    # 用于将钻研报告聚类为混合成分的 CRP
    xi[1:nStudies] ~ dCRP(alpha, size = nStudies)
    # 超参数
  
    theta ~ dflat() # 药物的影响 

运行 MCMC

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

inits <- list(gamma = rnorm(nStudies))

MCMC(code = BNP, data = data)

hist(samplesBNP[, 'theta'], xlab = expression(theta), main = 'avandia 的影响')
              main = "随机效应散布")
                   main = "随机效应散布")

# 推断出了多少个混合成分？xiRes <- samplesBNP[, xiCols].

次要推论仿佛对原始的参数化假如很持重。这可能是因为没有太多证据表明随机效应散布中不足正态性。

参考文献

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.

最受欢迎的见解

1.matlab 应用贝叶斯优化的深度学习

2.matlab 贝叶斯隐马尔可夫 hmm 模型实现

3.R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真

4.R 语言中的 block Gibbs 吉布斯采样贝叶斯多元线性回归

5.R 语言中的 Stan 概率编程 MCMC 采样的贝叶斯模型

6.Python 用 PyMC3 实现贝叶斯线性回归模型

7.R 语言应用贝叶斯 层次模型进行空间数据分析

8.R 语言随机搜寻变量抉择 SSVS 预计贝叶斯向量自回归（BVAR）模型

9.matlab 贝叶斯隐马尔可夫 hmm 模型实现