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

最近我们被客户要求撰写关于贝叶斯非参数模型的研究报告，包括一些图形和统计输出。

概述

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

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

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

 

登录后复制

data(faithful)

• 1.

 

登录后复制

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)

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.

 

登录后复制

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

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.
• 8.
• 9.
• 10.
• 11.
• 12.
• 13.
• 14.
• 15.
• 16.
• 17.
• 18.
• 19.

 

 

登录后复制

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

• 1.
• 2.

 

 

登录后复制

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

• 1.
• 2.

在这个模型下，对于一个新的观察

，后验预测分布是最佳密度估计（在平方误差损失下）。这个估计的样本可以很容易地从我们的MCMC产生的样本中计算出来。

 

 

登录后复制

# 参数的后验样本
 samples[, "alpha"]
# 平均值的后验样本
 samples[, grep('muTilde', colnames(samples))] # 聚类平均数的后验样本。
# 集群方差的后验样本
samples[, grep('s2Tilde', colnames(samples))] # 聚类成员的后验样本。
# 聚类成员关系的后验样本
samples [, grep('xi', colnames(samples))] # 聚类成员的后验样本。

hist(y, freq = FALSE,
     xlab = "标准化对数尺度上的等待时间")
##对标准化对数网格的密度进行点式估计

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.
• 8.
• 9.
• 10.
• 11.
• 12.

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

 

 

登录后复制

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

hist(faithful$waiting, freq = FALSE

• 1.
• 2.
• 3.

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

生成混合分布的样本

 

登录后复制

truncG <- outputG$trunc # G的截断水平



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

hist(probY70 )

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.

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

 

登录后复制

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

• 1.
• 2.
• 3.

 

登录后复制

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

• 1.
• 2.

 

 

 

 

登录后复制

cModel <- compile

• 1.

 

 

 

 

登录后复制

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

• 1.

 

 

 

登录后复制

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

• 1.
• 2.

 

 

登录后复制

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

• 1.
• 2.

从混合分布中生成样本

 

 

登录后复制

outputG <- getSamplesDPmeasure(cmcmc)

• 1.

 

 

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

 

 

登录后复制

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


hist(waiting, freq = FALSE

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.

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

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

模型

 

登录后复制

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)

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.
• 8.

注意，截断水平

已被设置为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

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.

 

 

 

 

登录后复制

compile(Model)

• 1.

 

 

 

 

登录后复制

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

• 1.
• 2.

 

 

 

 

登录后复制

MCMC(cmcmc, niter = 24000)

• 1.

 

 

登录后复制

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

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

• 1.
• 2.
• 3.
• 4.

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

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

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

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

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

模型的制定

 

登录后复制

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)

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.
• 8.
• 9.

运行MCMC

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

 

 

登录后复制

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

• 1.
• 2.

 

 

 

 

登录后复制

par(mfrow = c(1, 4)

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

• 1.
• 2.
• 3.
• 4.

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

用于meta分析的基于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() # 药物的影响

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.
• 8.
• 9.
• 10.
• 11.
• 12.
• 13.
• 14.
• 15.

运行MCMC

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

 

 

登录后复制

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

MCMC(code = BNP, data = data)

• 1.
• 2.
• 3.

 

 

 

 

登录后复制

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

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

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.

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

参考文献

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.

---

原创作者: u_14293657 转载于: https://blog.51cto.com/u_14293657/11312494