INLA教程入门
周瑜
2024-10-24
INLA
INLA 是一种快速的替代方法，用于拟合贝叶斯模型，相较于 MCMC，各有优劣。虽然 MCMC 是一种比较直观的方法，甚至在简单的场景下自己动手实现也不难，但 INLA 的算法对我来说一开始有些难理解。
现在我觉得好一些了，但当我在做博士论文时第一次学习 INLA 时，几乎没有太多友好的入门材料。我花了好几个月的时间，反复研究 INLA 的经典论文中的证明，试图弄清其中的道理。
为了帮助自己，我整理了一份面向应用统计学家的总结——就像我自己一样，熟悉正态分布、泰勒级数展开和贝叶斯定理等理论，但可能最近没怎么碰过 Hessian 矩阵。
这是一个非常简化的解释；经典论文非常详尽且清晰，值得一读（甚至读五十遍），如果你想更深入地理解其中的细节。这些细节既有趣又重要。但这是我试图给出的“精髓”，也是我刚开始学习时希望能有的东西。
简介：为什么选择INLA
层次模型广泛应用于各个学科中，用来表示数据中复杂的依赖结构。在贝叶斯框架下，可以通过合适的先验分布对模型参数和潜在变量或过程中的不确定性进行编码。
尽管后验分布几乎从未有显示解，但已经通过马尔可夫链蒙特卡洛（MCMC）方法开发出通用解决方案，比如像 WinBUGS 这样的软件。然而，MCMC 方法速度较慢，扩展性差，对于一些复杂模型，可能会失败（模型无法收敛）。一些新软件（如 JAGS、Stan）尝试解决这些问题。
MCMC 是一种渐近精确的方法，而 INLA 是一种近似方法。有时使用 INLA 会引来质疑或问题，例如“那近似误差怎么办？你对结果有多大的信心？” 的确，INLA 有时会失效，但要记住我们生活在现实世界中，不是在理想化的渐近环境中。经验上，许多模拟研究和我个人的多次观察都表明，MCMC 误差和 INLA 误差往往非常相似。
INLA 是马尔可夫链蒙特卡洛方法的一种快速替代方法，适用于广泛的潜在高斯模型（LGMs）。许多常见的模型可以重新构造成 LGMs 的形式，比如广义线性模型（GLM(M)）、广义加性模型（GAM(M)）、时间序列、空间模型、测量误差模型等。
要理解 INLA 的核心工作原理，我们需要熟悉以下概念： 1. 贝叶斯推断 2. 潜在高斯模型（LGMs） 3. 高斯马尔科夫随机场（GMRFs） 4. 拉普拉斯近似
如果你有兴趣学习 INLA，我假设你已经熟悉贝叶斯推断，但我们还是会先进行一个非常简短的回顾。
贝叶斯推断
在贝叶斯模型中，我们通常希望得到模型的后验分布（例如，参数在给定数据的条件下的分布），或者预测后验分布（用于预测/预测——即给定观测值的情况下新值的分布）。让我们专注于前者，后者可以逻辑推导。
后验分布等于数据似然函数乘以前验分布，再除以归一化常数（使得后验分布积分为1）：
\[ p(\theta \mid \mathbf{y})=\frac{p(\mathbf{y} \mid \theta) p(\theta)}{\int p(\mathbf{y} \mid \theta) p(\theta) d \theta} \propto p(\mathbf{y} \mid \theta) p(\theta) \]
在频率派框架中，我们通常使用数值方法（如牛顿-拉夫森法）最大化数据的似然函数，来获得给定参数的点估计（我们将参数视为固定但未知的）。然后我们使用理论重抽样的思想来估计围绕该参数估计的相应置信区间。实际上，我们通常使用显式重抽样的形式，如自助法，但这是另一个话题。在贝叶斯分析中，我们得到参数的后验分布（我们将其视为随机变量），可以为其提供总结统计量（如均值、中位数或众数），并直接通过分位数获得可信区间。
如果我们对每个参数使用共轭先验（即数学上简便的典型分布），那么我们可能可以以闭式形式确定后验分布。通常情况下，分母中的积分是不可解的，因此我们使用诸如 MCMC 的数值方法，从条件分布中采样，并估计每个感兴趣参数的边际分布。
潜在高斯模型
一般形式的潜在高斯模型（LGM）可以表示为：
\[ \mathbf{y} \mid \mathbf{x}, \theta_2 \sim \prod_i p\left(y_i \mid \eta_i, \theta_2\right) \quad \text{"似然函数"} \]
\[ \mathbf{x} \mid \theta_1 \sim p\left(\mathbf{x} \mid \theta_1\right)=\mathcal{N}(0, \mathbf{\Sigma}) \quad \text{"潜在场"} \]
\[ \theta=\left[\theta_1, \theta_2\right]^T \sim p(\theta) \quad \text{"超先验"} \]
其中，\(\mathbf{y}\) 是观测数据集，\(\mathbf{x}\) 不是协变量，而是线性预测器中所有参数的联合分布（包括其自身），\(\theta\) 是潜在场中的非高斯超参数。
那么，我们熟悉的模型，如我之前提到的广义线性（混合）模型（GLM(M)），如何融入这个框架呢？我们可以考虑广义（线性/加性）（混合）模型的一个非常通用的形式：
\[ \begin{gathered} \mathbf{y} \sim \prod_i^N p\left(y_i \mid \mu_i\right) \\ g\left(\mu_i\right)=\eta_i=\alpha+\sum_{k=1}^{n_\beta} \beta_k z_{k i}+\sum_{j=1}^{n_f} f^{(j)}\left(w_{j i}\right)+\varepsilon_i \end{gathered} \]
其中 \(g(\cdot)\) 是连接函数，\(\alpha\) 是截距，\(\beta\) 是（线性）协变量 \(\mathbf{z}\) 的熟悉回归参数，\(\mathbf{f}=f_1(\cdot), f_2(\cdot), \ldots, f_{n_\beta}(\cdot)\) 是一些协变量 \(\mathbf{w}\) 上的函数集合，它们可能是非线性的，比如随机（独立同分布）效应、空间或时间相关效应、平滑样条等。
线性回归是一个特殊情况，当结果变量为高斯分布，连接函数为恒等函数，并且 \(\mathbf{f}() = 0\) 时。在经典的BYM空间模型中，\(f_1(\cdot) \sim CAR\)，而 \(f_2(\cdot) \sim N(0, \sigma_{f_2}^2)\)，连接函数通常是对数函数或逻辑函数，因为你通常有二项或泊松分布的结果变量。许多熟悉的模型都是这个最通用形式的特例（有些特例更为明显）。
当我们假设这些参数具有联合高斯分布时，这个模型就是LGM。这可以通过给每个参数设定正态先验来实现：
\[ \mathbf{x}=[\eta, \alpha, \beta, \mathbf{f}(\cdot)] \sim \mathcal{N}(0, \mathbf{\Sigma}) \]
如果我们可以假设 \(\mathbf{x}\) 中的条件独立性，那么这个潜在场 \(\mathbf{x}\) 就是一个高斯马尔可夫随机场（GMRF）。在许多情况下，这个假设是合理的。我将在下一个部分介绍GMRF时对此做进一步讨论。需要注意的是，\(\mathbf{x}\) 的维度通常很大，因为它通过线性预测器与数据规模成正比，但非高斯超参数 \(\theta\) 的维度通常较小。
高斯马尔科夫随机场(GMRF)
GMRF 是具有马尔可夫性质的多元正态分布随机向量：例如，当 \(i \neq j\) 时，\(x_i \perp x_j \mid \mathbf{x}_{ij}\)。这里的符号 \(-ij\) 指的是“除了 \(i\) 和 \(j\) 之外的所有元素”。空间或时间自回归过程就是这种马尔可夫性质的一个自然例子。
回忆一下，精度是方差的倒数。同样地，精度矩阵 \((\mathbf{Q})\) 是协方差矩阵 \((\mathbf{\Sigma})\) 的逆矩阵，例如 \(\mathbf{x} \sim \mathcal{N}(0, \mathbf{\Sigma})\)，其中 \(\mathbf{Q}=\mathbf{\Sigma}^{-1}\)。此外，矩阵求逆在非平凡情况下是一个复杂且计算密集的操作，计算复杂度为 \(\mathcal{O}(n^3)\)。稀疏矩阵（含有大量零元素）则更容易求逆。例如，对角矩阵是稀疏的，且其求逆实际上就是求元素的倒数，单位矩阵的逆矩阵本身就是单位矩阵：\(I^{-1}=I\)。即使某些非对角元素非零，零元素的存在也大大简化了计算。
我们已经非常熟悉如下结论：
\[ x_i \perp x_j \Leftrightarrow \boldsymbol{\Sigma}_{ij}=0 \]
因此，如果我们希望得到稀疏的协方差矩阵，必须假设参数之间有很强的边缘独立性，而这通常是不合理的假设。然而，条件独立性（通过马尔可夫性质）则往往是合理的。GMRF 方法的强大之处在于 Havard Rue 等人展示了如何将条件独立性编码进精度矩阵，并且如何利用这一点来提高计算效率。经过大量推导，可以得出如下结果：
\[ x_i \perp x_j \mid \mathbf{x}_{-ij} \Leftrightarrow \mathbf{Q}_{ij}=0 \]
因此，在 GMRF 中的马尔可夫假设导致了一个稀疏的精度矩阵。这种稀疏性也传递到 Cholesky 分解矩阵中 \(\left(\mathbf{Q}=\mathbf{L} \mathbf{L}^{T}\right)\)，从而大大加快了计算速度。由于这是大矩阵，因此它在计算上是最重要的部分。
拉普拉斯近似(Laplace approximations)
回忆基本的泰勒级数展开，其中一个函数在某个点 \(a\) 处可以展开为一个（有时是无穷的）项的和，使用有限数量的项可以作为一个近似（通常使用前三项，直到二阶导数）。
\[ f(x)=f(a)+\frac{f^{\prime}(a)}{1!}(x-a)+\frac{f^{\prime\prime}(a)}{2!}(x-a)^2+\frac{f^{\prime\prime\prime}(a)}{3!}(x-a)^3+\ldots \]
假设我们有基本的抛物线 \(y=x^2\)。在 \(a=2\) 处展开（这里可以选择任意值）：
\[ f(x)=x^2, \quad f^{\prime}(x)=2x, \quad f^{\prime\prime}(x)=2, \quad f^{\prime\prime\prime}(x)=0 \]
因此：
\[ f(x)=x^2=2^2+2(2)(x-2)+\frac{2}{2}(x-2)^2 \]
# 定义函数 f(x) = x^2
f <- function(x) {
x^2
}
# 定义泰勒级数
taylor_series <- function(x) {
2^2 + 2 * 2 * (x - 2) + (2/2) * (x - 2)^2
}
# 定义 x 的范围
x <- seq(-200, 200, length.out = 400)
# 计算 f(x) 和泰勒级数的值
y_f <- f(x)
y_taylor <- taylor_series(x)
# 绘制图形
plot(x, y_f, type = 'l', col = 'blue', lwd = 4,
xlab = 'x', ylab = 'y',
xlim = c(-200, 200), ylim = c(0, 40000),
main = 'Plot of f(x) and its Taylor Series')
lines(x, y_taylor, col = 'red', lty = 2, lwd = 2)
# 添加图例
legend("topright", legend = c("f(x) = x^2", "Taylor Series"),
col = c("blue", "red"), lty = c(1, 2), lwd = 2)
# 添加网格
grid()
拉普拉斯近似用正态分布估计任何分布（PDF）。它使用泰勒级数展开式的前三项（二次函数）在一个函数的模态点 \(\hat{x}\) 处对 \(\log g(x)\) 进行近似（取对数简化了微分过程）：
\[ \log g(x) \approx \log g(\hat{x}) + \frac{\partial \log g(\hat{x})}{\partial x}(x - \hat{x}) + \frac{\partial^2 \log g(\hat{x})}{2 \partial x^2}(x - \hat{x})^2 \]
由于模态处的导数为零，上述表达式简化了，如果我们说基于曲率估计的方差是 \(\hat{\sigma}^2 = -\frac{1}{\frac{\partial^2 \log g(\hat{x})}{2 \partial x^2}}\)，那么我们可以将这些表达式重新写为：
\[ \log g(x) \approx \log g(\hat{x}) - \frac{1}{2 \sigma^2}(x - \hat{x})^2 \]
这应该变得更加熟悉了！我们可以使用这个结果进行正态近似。两边取指数和积分并移出常数项：
\[ \int g(x) dx = \int \exp [\log g(x)] dx \approx \text { 常数 } \int \exp \left[-\frac{(x - \hat{x})^2}{2 \hat{\sigma^2}}\right] dx \]
使用拉普拉斯近似，分布 \(f(x)\) 近似为正态分布，均值为 \(\hat{x}\)，可以通过解方程 \(f^{\prime}(x) = 0\) 找到模态，方差为 \(\hat{\sigma}^2 = -\frac{1}{f^{\prime \prime}(x)}\) （在模态处）。还不完全信服？这里用 \(\chi^2\) 分布做一个快速的例子（因为它容易手动求导；分母中的伽马函数项只是常数）。
\[ \begin{gathered} f(x; k) = \frac{x^{(k / 2 - 1)} e^{-x / 2}}{2^{k / 2} \Gamma\left(\frac{k}{2}\right)}, \quad x \geq 0 \\ \log f(x) = \left(\frac{k}{2} - 1\right) \log x - \frac{x}{2} \\ \log f^{\prime}(x) = \frac{\frac{k}{2} - 1}{x} - \frac{1}{2} = 0 \\ \log f^{\prime \prime}(x) = -\frac{\frac{k}{2} - 1}{x^2} \\ \therefore \chi_k^2 \stackrel{L A}{\sim} N\left(\hat{x} = k - 2, \hat{\sigma}^2 = 2(k - 2)\right) \end{gathered} \]
当自由度较高时（即分布更接近正态分布时），这个近似显然效果更好。这是一个非常简单的单变量例子，目的是让你相信，当分布大致接近正态分布时，这种方法往往效果很好。正如你所预料的那样，这个方法可以扩展到多变量情况，在这种情况下使用的是用于描述曲率的Hessian矩阵，接下来就涉及大量的数学推导了……
# 加载必要的库
if (!require("ggplot2")) install.packages("ggplot2")## 载入需要的程序包：ggplot2## Warning: 程序包'ggplot2'是用R版本4.4.1 来建造的if (!require("gridExtra")) install.packages("gridExtra")## 载入需要的程序包：gridExtralibrary(ggplot2)
library(gridExtra)
# 定义卡方分布的密度函数和拉普拉斯近似的正态分布参数
chi2_pdf <- function(x, k) {
return(dchisq(x, df = k))
}
laplace_approx_pdf <- function(x, k) {
mu <- k - 2
sigma2 <- 2 * (k - 2)
return(dnorm(x, mean = mu, sd = sqrt(sigma2)))
}
# 自由度列表
degrees_of_freedom <- c(6, 15, 25, 50)
# 定义 x 的范围
x <- seq(0, 100, length.out = 1000)
# 创建绘图函数
plot_density <- function(k, x_limits, y_limits) {
data <- data.frame(x = x,
chi2_density = chi2_pdf(x, k),
laplace_density = laplace_approx_pdf(x, k))
ggplot(data, aes(x = x)) +
geom_line(aes(y = chi2_density, color = "Chi^2 PDF"), linewidth = 1) +
geom_line(aes(y = laplace_density, color = "Laplace Approximation"), linetype = "dashed", linewidth = 1) +
labs(title = paste("k =", k),
x = "x", y = "Density") +
scale_color_manual(values = c("Chi^2 PDF" = "blue", "Laplace Approximation" = "red")) +
theme_minimal() +
theme(legend.title = element_blank(), legend.position = "top") +
coord_cartesian(xlim = x_limits, ylim = y_limits)
}
# 自由度列表和对应的坐标轴范围
degrees_of_freedom <- c(6, 15, 25, 50)
x_limits_list <- list(c(0, 30), c(0, 30), c(0, 50), c(0, 100))
y_limits_list <- list(c(0, 0.14), c(0, 0.08), c(0, 0.06), c(0, 0.04))
# 创建每个自由度的图
plots <- lapply(1:4, function(i) plot_density(degrees_of_freedom[i], x_limits_list[[i]], y_limits_list[[i]]))
# 将4个图排列成2行2列的布局
grid.arrange(grobs = plots, nrow = 2, ncol = 2)
# 设置图形的布局为两行两列
par(mfrow = c(2, 2))
# 定义卡方分布和拉普拉斯近似的函数
f <- function(x, k) {
return(dchisq(x, df = k))
}
laplace_approx <- function(x, k) {
mu <- k - 2
sigma2 <- 2 * (k - 2)
return(dnorm(x, mean = mu, sd = sqrt(sigma2)))
}
# 定义要绘制的自由度列表和相应的坐标范围
degrees_of_freedom <- c(6, 15, 25, 50)
x_limits <- list(c(0, 30), c(0, 30), c(0, 50), c(0, 100))
y_limits <- list(c(0, 0.14), c(0, 0.08), c(0, 0.06), c(0, 0.04))
# 绘制每个自由度对应的图
for (i in 1:4) {
k <- degrees_of_freedom[i]
x <- seq(x_limits[[i]][1], x_limits[[i]][2], length.out = 200)
# 计算 f(x) 和拉普拉斯近似的值
y_f <- f(x, k)
y_laplace <- laplace_approx(x, k)
# 绘制卡方分布和拉普拉斯近似曲线
plot(x, y_f, type = 'l', col = 'blue', lwd = 2,
xlab = 'x', ylab = 'Density',
xlim = x_limits[[i]], ylim = y_limits[[i]],
main = paste("Degrees of Freedom =", k))
lines(x, y_laplace, col = 'red', lty = 2, lwd = 2)
# 添加图例
# legend("topright", legend = c("Chi^2 PDF", "Laplace Approximation"),
# col = c("blue", "red"), lty = c(1, 2), lwd = 1)
# 添加网格
grid()
}
上图展示了自由度分别为 k=6,15,25,50
时的卡方分布密度函数（蓝线）以及通过拉普拉斯近似得到的正态分布近似（红色虚线）。随着自由度的增加，卡方分布逐渐向正态分布靠近，说明拉普拉斯近似在高自由度情况下效果更好
整合：集成嵌套拉普拉斯近似法
我们已经讨论了许多不同的部分，但对于我们的拉普拉斯高斯模型 (LGMs)，我们希望获得潜在场（例如回归参数）元素的边际分布，以及可能来自超先验分布中的一些元素的边际分布（例如自回归模型中的相关参数或随机效应的方差）——我们希望拟合一个分层模型！
\[ \begin{gathered} p\left(x_j \mid \mathbf{y}\right)=\int p\left(x_j, \theta \mid \mathbf{y}\right) d \theta=\int p\left(x_j \mid \theta, \mathbf{y}\right) p(\theta \mid \mathbf{y}) d \theta \\ p\left(\theta_k \mid \mathbf{y}\right)=\int p(\theta \mid \mathbf{y}) d \theta_{-k} \end{gathered} \]
因此，我们需要近似的项为：\(p(\theta \mid \mathbf{y})\) 和 \(p\left(x_j \mid \theta, \mathbf{y}\right)\)。第一个项可以用来估计我们对 \(\theta\) 感兴趣的所有边际分布，同时也需要用来估计潜在场的边际分布。
回顾概率论的基本规则：
\[ \begin{aligned} p(b) & =\frac{p(a, b)}{p(a \mid b)} \\ p(b \mid c) & =\frac{p(a, b \mid c)}{p(a \mid b, c)} \end{aligned} \]
我们可以利用此规则重新写出边际分布：
\[ p(\theta \mid \mathbf{y})=\frac{p(\mathbf{x}, \theta \mid \mathbf{y})}{p(\mathbf{x} \mid \theta, \mathbf{y})} \]
我们可以展开分子部分（细节省略）并对未知的分母使用拉普拉斯近似：
\[ \left.p(\theta \mid \mathbf{y}) \approx \frac{p(\mathbf{y} \mid \mathbf{x}, \theta) p(\mathbf{x} \mid \theta) p(\theta)}{\tilde{p}(\mathbf{x} \mid \theta, \mathbf{y})}\right|_{x=x^*(\theta)}=\tilde{p}(\theta \mid \mathbf{y}) \]
其中 \(x^*(\theta)\) 是在给定 \(\theta\) 配置下 \(\mathbf{x}\) 的模态值。获取 \(p\left(x_j \mid \theta, \mathbf{y}\right)\) 也使用相同的逻辑，但更加复杂，因为 \(\mathbf{x}\) 的维数要大得多。
最明显的选择是直接对 \(p\left(x_j \mid \theta, \mathbf{y}\right)\) 使用正态近似，因为我们已经在上面的分母中计算了该项。此方法非常快速，但假设过于强烈，结果通常较差（这些条件分布往往呈现轻微偏态或重尾特性）。
或者，我们可以将 \(\mathbf{x}\) 分解为 \(\mathbf{x}=\left[x_j, \mathbf{x}_{-j}\right]\)，并对潜在场 \(\mathbf{x}\) 的每个元素使用拉普拉斯近似：
\[ p\left(x_j \mid \theta, \mathbf{y}\right) \propto \frac{p(\mathbf{x}, \theta \mid \mathbf{y})}{p\left(\mathbf{x}_{-j} \mid x_j, \theta, \mathbf{y}\right)} \propto \frac{p(\theta) p(\mathbf{x} \mid \theta) p(\mathbf{y} \mid \mathbf{x})}{p\left(\mathbf{x}_{-j} \mid x_j, \theta, \mathbf{y}\right)} \]
这种方法效果很好，因为条件分布 \(p\left(\mathbf{x}_{-j} \mid x_j, \theta, \mathbf{y}\right)\) 通常接近正态，但计算量更大。
第三种方法（INLA软件中的默认方法）使用一种简化的拉普拉斯近似（可以理解为前两种方法的折中），这种方法计算效率极高且近似效果很好。它使用泰勒展开到三阶来估计 \(p\left(x_j \mid \theta, \mathbf{y}\right)\) 的分子和分母。
在 INLA 中可以使用其他两种方法。我个人在开始时使用近似拟合模型，因为它计算速度非常快。近似结果和MCMC方法（或真实值，在做模拟时）大体一致，但经常有显著差异。所以我在建模和验证时使用近似，然后在最终模型中使用简化或完整的拉普拉斯近似。我个人从未遇到过能够在简化和完整拉普拉斯近似之间显著区分的情况。
INLA 算法使用类似牛顿法的方式来探索超参数的联合后验分布 \(\tilde{p}(\theta \mid \mathbf{y})\)，从而找到合适的数值积分点：
\[ p\left(x_j \mid \mathbf{y}\right) \approx \sum_{h=1}^H \tilde{p}\left(x_j \mid \theta_h^*, \mathbf{y}\right) \hat{p}\left(\theta_h^* \mid \mathbf{y}\right) \Delta_h \]
现在我们可以理解 INLA 的核心步骤了！
集成 (Integrated)：使用数值积分
嵌套 (Nested)：因为需要 \(p(\theta \mid \mathbf{y})\) 才能得到 \(p\left(x_j \mid \mathbf{y}\right)\)
拉普拉斯近似 (Laplace Approximations)：用于获得正态近似的参数的方法
2011年，Finn Lindgren等人的研究提出了一种基于随机偏微分方程 (SPDE) 的INLA方法，这种方法允许在具有曲面的点数据上进行更真实的空间建模，使得大范围地理区域能够使用连续的协方差函数来进行建模。今天我不会讨论SPDE，但未来可能会专门发帖介绍。如果你刚开始学习 INLA，除非你只对 SPDE部分感兴趣，否则我建议先掌握基础知识。毕竟，随机微积分并不是我们大多数应用研究人员每天会用到的东西。
我知道刚才的信息量很大——但实际上我只是触及了皮毛！我推荐阅读 2009 年的原始 INLA 论文（参考文献中可见）。现在，我们切换到 R-INLA 代码的展示，并在一些简单的示例中对比其输出与 MCMC 的结果。
使用INLA
R-INLA 依赖于一个独立的 C 库（GMRF），因此它不在 CRAN 上发布，但你可以在 http://www.r-inla.org 获取所有相关信息。
我将通过两个非常简单的示例来演示如何使用 rjags 进行 MCMC
分析，以及如何使用 R-INLA 进行 INLA 分析。
简单线性回归
\[ y_i=\alpha+\beta x_i+\epsilon_i,\epsilon_i\sim N(0,\sigma) \\ E\left(y_i \mid x_i\right)=\alpha+\beta x_i \]
- Simulated data in R
N<-100#500,5000,25000,10000
x<-rnorm(N ,mean =6, sd = 2)
y<-rnorm(N ,mean =x, sd = 1)
data<- list(x=x,y=y,N=N)- JAGS code
library(R2jags)## Warning: 程序包'R2jags'是用R版本4.4.1 来建造的## 载入需要的程序包：rjags## 载入需要的程序包：coda## Linked to JAGS 4.3.1## Loaded modules: basemod,bugs##
## 载入程序包：'R2jags'## The following object is masked from 'package:coda':
##
## traceplotlibrary(rjags)
library(coda)
model<-function(){
for(i in 1:N){
y[i]~dnorm(mu[i],tau)
mu[i]<-alpha+beta*x[i]
}
alpha~dnorm(0,0.001)
beta~dnorm(0,0.001)
tau~dnorm(0.01,0.01)
}
params<-c("alpha","beta","tau","mu")
jagsmodel<-jags(
data = data,
param = params,
n.chains = 3,
n.iter = 50000,
n.burnin = 5000,
model.file = model
)## module glm loaded## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 100
## Unobserved stochastic nodes: 3
## Total graph size: 407
##
## Initializing model# 提取 alpha 样本值
alpha_samples <- jagsmodel$BUGSoutput$sims.list$alpha
# 绘制频率直方图，设置坐标范围
hist(alpha_samples,
main = "Posterior samples of alpha",
xlab = "alpha",
xlim = c(-4, 2), # 设置 x 轴范围
ylim = c(0, 1.4), # 设置 y 轴范围
breaks = 30, # 可以根据需要调整
col = "lightblue",
border = "black",
freq = FALSE) # 显示频率
- INLA code
library(INLA)## Warning: 程序包'INLA'是用R版本4.4.1 来建造的## 载入需要的程序包：Matrix## 载入需要的程序包：sp## Warning: 程序包'sp'是用R版本4.4.1 来建造的## This is INLA_24.06.27 built 2024-06-27 02:36:04 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - List available models/likelihoods/etc with inla.list.models()
## - Use inla.doc(<NAME>) to access documentation##
## 载入程序包：'INLA'## The following object is masked _by_ '.GlobalEnv':
##
## fresult<-inla(y~x,
family = "gaussian",
data =data,
control.predictor = list(link = 1),
control.compute = list(config = TRUE)
)
result$summary.fixed## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -0.1382005 0.39380920 -0.9120186 -0.1382008 0.6356188 -0.1382008
## x 1.0101877 0.06099011 0.8903445 1.0101877 1.1300306 1.0101877
## kld
## (Intercept) 8.307511e-09
## x 8.307507e-09result$summary.linear.predictor## mean sd 0.025quant 0.5quant 0.975quant mode
## Predictor.001 1.3066653 0.3111000 0.6954858 1.3066652 1.9178458 1.3066652
## Predictor.002 10.4149807 0.2851400 9.8548007 10.4149809 10.9751599 10.4149809
## Predictor.003 6.2123963 0.1161748 5.9841618 6.2123963 6.4406308 6.2123963
## Predictor.004 9.4835274 0.2349452 9.0219592 9.4835275 9.9450950 9.4835275
## Predictor.005 3.2011713 0.2094762 2.7896394 3.2011712 3.6127038 3.2011712
## Predictor.006 3.2160762 0.2087284 2.8060135 3.2160761 3.6261395 3.2160761
## Predictor.007 4.4961011 0.1507802 4.1998819 4.4961010 4.7923206 4.4961010
## Predictor.008 5.7102327 0.1182499 5.4779215 5.7102327 5.9425439 5.7102327
## Predictor.009 7.5568643 0.1456705 7.2706831 7.5568643 7.8430452 7.5568643
## Predictor.010 8.4340682 0.1826031 8.0753301 8.4340683 8.7928059 8.4340683
## Predictor.011 7.9884114 0.1627428 7.6686904 7.9884115 8.3081321 7.9884115
## Predictor.012 5.6252640 0.1193596 5.3907728 5.6252640 5.8597553 5.6252640
## Predictor.013 4.5487792 0.1487701 4.2565090 4.5487791 4.8410496 4.5487791
## Predictor.014 5.6518098 0.1189904 5.4180441 5.6518098 5.8855757 5.6518098
## Predictor.015 6.7947250 0.1234032 6.5522898 6.7947250 7.0371600 6.7947250
## Predictor.016 6.7565741 0.1226353 6.5156475 6.7565741 6.9975005 6.7565741
## Predictor.017 5.3693672 0.1239271 5.1259030 5.3693672 5.6128316 5.3693672
## Predictor.018 8.8520090 0.2027021 8.4537849 8.8520091 9.2502327 8.8520091
## Predictor.019 6.6447492 0.1206093 6.4078028 6.6447492 6.8816956 6.6447492
## Predictor.020 3.7253064 0.1839910 3.3638420 3.7253063 4.0867712 3.7253063
## Predictor.021 3.6582985 0.1871452 3.2906377 3.6582985 4.0259599 3.6582985
## Predictor.022 3.4226608 0.1984921 3.0327079 3.4226607 3.8126141 3.4226607
## Predictor.023 5.5820574 0.1200039 5.3463005 5.5820574 5.8178144 5.5820574
## Predictor.024 8.9756257 0.2088589 8.5653059 8.9756258 9.3859449 8.9756258
## Predictor.025 7.7410529 0.1526485 7.4411630 7.7410530 8.0409426 7.7410530
## Predictor.026 5.8523681 0.1168731 5.6227618 5.8523681 6.0819744 5.8523681
## Predictor.027 10.4254064 0.2857144 9.8640979 10.4254065 10.9867140 10.4254065
## Predictor.028 10.0081342 0.2629188 9.4916095 10.0081344 10.5246583 10.0081344
## Predictor.029 6.8109723 0.1237417 6.5678720 6.8109723 7.0540725 6.8109723
## Predictor.030 8.4212983 0.1820089 8.0637274 8.4212983 8.7788687 8.4212983
## Predictor.031 4.3885287 0.1550062 4.0840071 4.3885286 4.6930506 4.3885286
## Predictor.032 8.1850684 0.1712710 7.8485931 8.1850685 8.5215434 8.1850685
## Predictor.033 5.8746935 0.1167126 5.6454027 5.8746935 6.1039845 5.8746935
## Predictor.034 7.0354539 0.1290912 6.7818442 7.0354539 7.2890634 7.0354539
## Predictor.035 5.2501989 0.1266407 5.0014035 5.2501989 5.4989945 5.2501989
## Predictor.036 6.9976067 0.1281047 6.7459350 6.9976067 7.2492782 6.9976067
## Predictor.037 7.2181188 0.1342962 6.9542835 7.2181189 7.4819539 7.2181189
## Predictor.038 6.1339649 0.1159817 5.9061099 6.1339649 6.3618199 6.1339649
## Predictor.039 5.6729955 0.1187103 5.4397799 5.6729955 5.9062112 5.6729955
## Predictor.040 5.6874568 0.1185267 5.4546020 5.6874568 5.9203117 5.6874568
## Predictor.041 6.9036141 0.1258006 6.6564690 6.9036141 7.1507590 6.9036141
## Predictor.042 5.8533703 0.1168656 5.6237788 5.8533703 6.0829619 5.8533703
## Predictor.043 6.7749631 0.1230007 6.5333188 6.7749631 7.0166074 6.7749631
## Predictor.044 4.9631608 0.1345346 4.6988574 4.9631608 5.2274645 4.9631608
## Predictor.045 6.7282750 0.1220906 6.4884185 6.7282750 6.9681313 6.7282750
## Predictor.046 5.6324576 0.1192576 5.3981670 5.6324576 5.8667484 5.6324576
## Predictor.047 6.1791887 0.1160696 5.9511610 6.1791887 6.4072164 6.1791887
## Predictor.048 3.7628592 0.1822386 3.4048375 3.7628591 4.1208813 3.7628591
## Predictor.049 6.4829814 0.1182988 6.2505741 6.4829814 6.7153885 6.4829814
## Predictor.050 7.7490539 0.1529627 7.4485467 7.7490540 8.0495608 7.7490540
## Predictor.051 4.2559953 0.1604205 3.9408369 4.2559952 4.5711540 4.2559952
## Predictor.052 5.3481068 0.1243851 5.1037426 5.3481067 5.5924710 5.3481067
## Predictor.053 7.5682787 0.1460881 7.2812772 7.5682788 7.8552800 7.5682788
## Predictor.054 5.3538187 0.1242609 5.1096986 5.3538187 5.5979390 5.3538187
## Predictor.055 3.8758317 0.1770369 3.5280291 3.8758316 4.2236346 3.8758316
## Predictor.056 2.8362573 0.2281252 2.3880880 2.8362572 3.2844273 2.8362572
## Predictor.057 5.4325951 0.1226339 5.1916714 5.4325951 5.6735190 5.4325951
## Predictor.058 7.5108824 0.1440094 7.2279646 7.5108824 7.7937998 7.5108824
## Predictor.059 6.6095246 0.1200425 6.3736917 6.6095247 6.8453575 6.6095247
## Predictor.060 6.3549196 0.1170152 6.1250342 6.3549196 6.5848050 6.3549196
## Predictor.061 5.6441153 0.1190953 5.4101434 5.6441153 5.8780873 5.6441153
## Predictor.062 5.5799484 0.1200367 5.3441270 5.5799483 5.8157698 5.5799483
## Predictor.063 6.7526294 0.1225581 6.5118545 6.7526295 6.9934042 6.7526295
## Predictor.064 6.6231976 0.1202584 6.3869406 6.6231976 6.8594545 6.6231976
## Predictor.065 3.9822798 0.1722384 3.6439043 3.9822797 4.3206556 3.9822797
## Predictor.066 5.7973593 0.1173337 5.5668482 5.7973593 6.0278704 5.7973593
## Predictor.067 5.9739353 0.1161852 5.7456805 5.9739353 6.2021900 5.9739353
## Predictor.068 7.5979899 0.1471846 7.3088342 7.5979900 7.8871453 7.5979900
## Predictor.069 7.3708822 0.1391699 7.0974721 7.3708822 7.6442921 7.3708822
## Predictor.070 7.7300881 0.1522194 7.4310412 7.7300882 8.0291347 7.7300882
## Predictor.071 5.4716361 0.1218881 5.2321775 5.4716361 5.7110948 5.4716361
## Predictor.072 7.5626588 0.1458823 7.2760617 7.5626589 7.8492557 7.5626589
## Predictor.073 0.7147409 0.3444747 0.0379941 0.7147407 1.3914886 0.7147407
## Predictor.074 4.7506027 0.1414649 4.4726841 4.7506027 5.0285216 4.7506027
## Predictor.075 5.5608214 0.1203399 5.3244044 5.5608214 5.7972386 5.5608214
## Predictor.076 6.2362364 0.1162716 6.0078118 6.2362364 6.4646610 6.2362364
## Predictor.077 6.6492217 0.1206838 6.4121290 6.6492217 6.8863143 6.6492217
## Predictor.078 6.5788149 0.1195770 6.3438966 6.5788149 6.8137331 6.5788149
## Predictor.079 4.7884690 0.1401696 4.5130951 4.7884689 5.0638432 4.7884689
## Predictor.080 4.9414450 0.1352032 4.6758279 4.9414449 5.2070622 4.9414449
## Predictor.081 5.6543721 0.1189558 5.4206742 5.6543721 5.8880700 5.6543721
## Predictor.082 8.1299191 0.1688389 7.7982219 8.1299192 8.4616159 8.1299192
## Predictor.083 8.5960929 0.1902508 8.2223302 8.5960930 8.9698552 8.5960930
## Predictor.084 7.0520047 0.1295328 6.7975274 7.0520047 7.3064818 7.0520047
## Predictor.085 3.2625895 0.2064024 2.8570964 3.2625894 3.6680832 3.2625894
## Predictor.086 7.8310150 0.1562305 7.5240880 7.8310151 8.1379417 7.8310151
## Predictor.087 5.8222816 0.1171136 5.5922029 5.8222816 6.0523603 5.8222816
## Predictor.088 6.3622815 0.1170756 6.1322774 6.3622815 6.5922855 6.3622815
## Predictor.089 5.1039978 0.1304344 4.8477495 5.1039977 5.3602463 5.1039977
## Predictor.090 5.1069731 0.1303523 4.8508860 5.1069731 5.3630604 5.1069731
## Predictor.091 6.6658350 0.1209652 6.4281894 6.6658350 6.9034805 6.6658350
## Predictor.092 9.3579775 0.2283920 8.9092834 9.3579776 9.8066710 9.3579776
## Predictor.093 10.1929281 0.2729629 9.6566710 10.1929283 10.7291845 10.1929283
## Predictor.094 4.1028760 0.1669336 3.7749221 4.1028759 4.4308302 4.1028759
## Predictor.095 4.5060501 0.1503975 4.2105828 4.5060500 4.8015177 4.5060500
## Predictor.096 6.6173629 0.1201657 6.3812882 6.6173630 6.8534376 6.6173630
## Predictor.097 7.4590426 0.1421782 7.1797224 7.4590427 7.7383625 7.4590427
## Predictor.098 4.6964298 0.1433609 4.4147864 4.6964298 4.9780735 4.6964298
## Predictor.099 5.6292564 0.1193028 5.3948769 5.6292564 5.8636360 5.6292564
## Predictor.100 -0.5086399 0.4146768 -1.3233039 -0.5086401 0.3060254 -0.5086401
## kld
## Predictor.001 6.647494e-09
## Predictor.002 6.647612e-09
## Predictor.003 6.647503e-09
## Predictor.004 6.647557e-09
## Predictor.005 6.647507e-09
## Predictor.006 6.647508e-09
## Predictor.007 6.647544e-09
## Predictor.008 6.647596e-09
## Predictor.009 6.647511e-09
## Predictor.010 6.647306e-09
## Predictor.011 6.647354e-09
## Predictor.012 6.647577e-09
## Predictor.013 6.647435e-09
## Predictor.014 6.647558e-09
## Predictor.015 6.647615e-09
## Predictor.016 6.647492e-09
## Predictor.017 6.647689e-09
## Predictor.018 6.647334e-09
## Predictor.019 6.647617e-09
## Predictor.020 6.647524e-09
## Predictor.021 6.647493e-09
## Predictor.022 6.647522e-09
## Predictor.023 6.647448e-09
## Predictor.024 6.647680e-09
## Predictor.025 6.647710e-09
## Predictor.026 6.647903e-09
## Predictor.027 6.647504e-09
## Predictor.028 6.647497e-09
## Predictor.029 6.647487e-09
## Predictor.030 6.647444e-09
## Predictor.031 6.647463e-09
## Predictor.032 6.647250e-09
## Predictor.033 6.647554e-09
## Predictor.034 6.647492e-09
## Predictor.035 6.647588e-09
## Predictor.036 6.647531e-09
## Predictor.037 6.647601e-09
## Predictor.038 6.647476e-09
## Predictor.039 6.647569e-09
## Predictor.040 6.647444e-09
## Predictor.041 6.647547e-09
## Predictor.042 6.647329e-09
## Predictor.043 6.647632e-09
## Predictor.044 6.647569e-09
## Predictor.045 6.647589e-09
## Predictor.046 6.647720e-09
## Predictor.047 6.647829e-09
## Predictor.048 6.647521e-09
## Predictor.049 6.647435e-09
## Predictor.050 6.647612e-09
## Predictor.051 6.647571e-09
## Predictor.052 6.647560e-09
## Predictor.053 6.647520e-09
## Predictor.054 6.647567e-09
## Predictor.055 6.647525e-09
## Predictor.056 6.647520e-09
## Predictor.057 6.647526e-09
## Predictor.058 6.647735e-09
## Predictor.059 6.647667e-09
## Predictor.060 6.647596e-09
## Predictor.061 6.647747e-09
## Predictor.062 6.647452e-09
## Predictor.063 6.647591e-09
## Predictor.064 6.647284e-09
## Predictor.065 6.647524e-09
## Predictor.066 6.647562e-09
## Predictor.067 6.647375e-09
## Predictor.068 6.647496e-09
## Predictor.069 6.647759e-09
## Predictor.070 6.647677e-09
## Predictor.071 6.647571e-09
## Predictor.072 6.647762e-09
## Predictor.073 6.647491e-09
## Predictor.074 6.647571e-09
## Predictor.075 6.647543e-09
## Predictor.076 6.647876e-09
## Predictor.077 6.647590e-09
## Predictor.078 6.647226e-09
## Predictor.079 6.647441e-09
## Predictor.080 6.647437e-09
## Predictor.081 6.647718e-09
## Predictor.082 6.647448e-09
## Predictor.083 6.647516e-09
## Predictor.084 6.647555e-09
## Predictor.085 6.647516e-09
## Predictor.086 6.647554e-09
## Predictor.087 6.647311e-09
## Predictor.088 6.647481e-09
## Predictor.089 6.647710e-09
## Predictor.090 6.647508e-09
## Predictor.091 6.647290e-09
## Predictor.092 6.647546e-09
## Predictor.093 6.647533e-09
## Predictor.094 6.647464e-09
## Predictor.095 6.647451e-09
## Predictor.096 6.647832e-09
## Predictor.097 6.647566e-09
## Predictor.098 6.647508e-09
## Predictor.099 6.647547e-09
## Predictor.100 6.647487e-09plot(result$marginals.fixed$x,
type = "l",xlab= "beta",ylab = "density",xlim=c(0.8,1.6))
plot(result$marginals.fixed[[1]],
type = "l",xlab= "alpha",ylab = "density",xlim=c(-4,2))
# 设置图形区域
plot(result$marginals.fixed[[1]],
type = "l",
xlab = "Alpha",
ylab = "Density",
xlim = c(-4, 2),
ylim = c(0, 1.4),
col = "blue",
lwd = 2,
main = "Alpha") # 添加标题
# 添加 MCMC 的直方图
hist(alpha_samples,
main = "", # 直方图不需要标题
xlab = "Alpha",
xlim = c(-4, 2),
ylim = c(0, 1.4),
breaks = 30,
col = "lightblue",
border = "black",
freq = FALSE,
add = TRUE) # 添加到同一图形上
# 添加图例
legend("topright", legend = c("INLA", "MCMC"),
col = c("blue", "lightblue"), lwd = 2, fill = c(NA, "lightblue"),
bty = "n")
# 设置图形区域
plot(result$marginals.fixed[[2]],
type = "l",
xlab = "beta",
ylab = "Density",
xlim = c(0.8, 1.6),
ylim = c(0, 9),
col = "blue",
lwd = 2,
main = "beta") # 添加标题
# 添加 MCMC 的直方图
beta_samples <- jagsmodel$BUGSoutput$sims.list$beta
hist(beta_samples,
main = "", # 直方图不需要标题
xlab = "beta",
xlim = c(0.8, 1.6),
ylim = c(0, 9),
breaks = 30,
col = "lightblue",
border = "black",
freq = FALSE,
add = TRUE) # 添加到同一图形上
# 添加图例
legend("topright", legend = c("INLA", "MCMC"),
col = c("blue", "lightblue"), lwd = 2, fill = c(NA, "lightblue"),
bty = "n")
# 设置图形区域
# 提取 tau 的后验边际分布
tau_marginal <- result$marginals.hyperpar[[1]] # 根据顺序提取
# 绘制后验密度
plot(tau_marginal,
type = "l",
xlab = "Tau",
ylab = "Density",
main = "Posterior Density of Tau",
xlim = c(0, 1.5), # 设置 x 轴范围
ylim = c(0, 3), # 设置 y 轴范围
col = "blue",
lwd = 2)
# 添加 MCMC 的直方图
# 提取 tau 的样本
tau_samples <- jagsmodel$BUGSoutput$sims.list$tau
# 绘制 tau 的后验直方图
hist(tau_samples,
main = "Posterior Samples of Tau",
xlab = "Tau",
xlim = c(0, 1.5), # 设置 x 轴范围
ylim = c(0, 3), # 设置 y 轴范围
col = "lightblue",
border = "black",
breaks = 30, # 可以根据需要调整
freq = FALSE,
add = TRUE) # 显示密度
# 添加图例
legend("topright", legend = c("INLA", "MCMC"),
col = c("blue", "lightblue"), lwd = 2, fill = c(NA, "lightblue"),
bty = "n")
###带有iid随机效应的poisson广义线性模型(GLM) \[ \begin{gathered} E\left(y_i \mid x_i\right)=\operatorname{poisson}\left(\mu_i\right) \\ \log \left(\mu_i\right)=\alpha+\beta x_i+\nu_i \\ \nu_i \sim N\left(0, \tau_\nu\right) \end{gathered} \]
- Simulated data in R
N<-100#500,5000,25000,10000
x<-rnorm(N ,mean =5, sd = 1)
nu<-rnorm(N,0,0.1)
mu<-exp(1+0.5*x+nu)
y<-rpois(N,mu)
data<-list(x=x,y=y,nu=nu,N=N)- JAGS code
library(R2jags)
library(rjags)
library(coda)
model<-function(){
for(i in 1:N){
y[i]~dpois(mu[i])
log(mu[i])<-alpha+beta*x[i]+nu[i]
nu[i]~dnorm(0,tau.nu)
}
alpha~dnorm(0,0.001)
beta~dnorm(0,0.001)
tau.nu~dnorm(0.01,0.01)
}
params<-c("alpha","beta","tau.nu","mu")
jagsmodel<-jags(
data = data,
param = params,
n.chains = 3,
n.iter = 50000,
n.burnin = 5000,
model.file = model
)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 200
## Unobserved stochastic nodes: 3
## Total graph size: 607
##
## Initializing modeljagsmodel## Inference for Bugs model at "C:/Users/Admin/AppData/Local/Temp/RtmpoTFluy/model531848c73a53", fit using jags,
## 3 chains, each with 50000 iterations (first 5000 discarded), n.thin = 45
## n.sims = 3000 iterations saved. Running time = 72.98 secs
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## alpha 0.985 0.099 0.791 0.915 0.988 1.052 1.183 1.001 3000
## beta 0.497 0.018 0.461 0.485 0.497 0.510 0.533 1.001 3000
## mu[1] 43.192 0.810 41.600 42.656 43.198 43.729 44.774 1.001 2200
## mu[2] 23.512 0.544 22.478 23.131 23.519 23.890 24.554 1.001 3000
## mu[3] 24.690 0.709 23.326 24.193 24.691 25.177 26.083 1.001 3000
## mu[4] 20.168 0.534 19.147 19.790 20.174 20.535 21.212 1.001 3000
## mu[5] 33.734 0.637 32.512 33.295 33.746 34.179 34.960 1.001 3000
## mu[6] 59.843 1.495 56.865 58.844 59.828 60.829 62.688 1.002 2900
## mu[7] 17.025 0.642 15.796 16.577 17.027 17.461 18.319 1.001 3000
## mu[8] 29.359 0.659 28.088 28.897 29.372 29.821 30.616 1.001 3000
## mu[9] 31.718 0.596 30.571 31.305 31.728 32.132 32.866 1.001 3000
## mu[10] 14.934 0.551 13.883 14.548 14.936 15.307 16.042 1.001 3000
## mu[11] 42.655 0.734 41.213 42.160 42.659 43.146 44.098 1.001 2200
## mu[12] 47.846 0.825 46.227 47.291 47.851 48.400 49.473 1.001 2200
## mu[13] 32.655 0.652 31.406 32.202 32.667 33.110 33.902 1.001 3000
## mu[14] 74.809 2.024 70.831 73.467 74.786 76.158 78.670 1.002 3000
## mu[15] 19.764 0.757 18.317 19.234 19.766 20.275 21.288 1.001 3000
## mu[16] 26.035 0.521 25.035 25.672 26.043 26.400 27.031 1.001 3000
## mu[17] 34.138 0.583 32.991 33.742 34.125 34.546 35.271 1.001 2800
## mu[18] 19.652 0.601 18.500 19.231 19.648 20.067 20.837 1.001 3000
## mu[19] 55.115 1.287 52.561 54.267 55.115 55.956 57.569 1.002 2700
## mu[20] 30.650 0.666 29.382 30.182 30.660 31.116 31.910 1.001 3000
## mu[21] 20.197 0.530 19.186 19.821 20.203 20.564 21.235 1.001 3000
## mu[22] 84.493 2.287 80.000 82.978 84.468 86.017 88.855 1.002 3000
## mu[23] 14.517 0.527 13.513 14.149 14.518 14.876 15.569 1.001 3000
## mu[24] 8.882 0.483 7.979 8.538 8.880 9.201 9.867 1.000 3000
## mu[25] 26.522 0.559 25.448 26.130 26.531 26.915 27.582 1.001 3000
## mu[26] 48.013 0.860 46.321 47.452 48.020 48.588 49.704 1.001 2200
## mu[27] 61.266 1.309 58.694 60.398 61.275 62.130 63.805 1.002 2500
## mu[28] 15.773 0.621 14.586 15.336 15.772 16.189 17.025 1.001 3000
## mu[29] 16.072 0.611 14.903 15.644 16.074 16.486 17.305 1.001 3000
## mu[30] 12.414 0.516 11.434 12.054 12.412 12.754 13.462 1.001 3000
## mu[31] 38.222 0.665 36.913 37.766 38.209 38.686 39.501 1.001 3000
## mu[32] 17.187 0.598 16.049 16.769 17.188 17.593 18.370 1.001 3000
## mu[33] 27.540 0.568 26.449 27.144 27.548 27.940 28.627 1.001 3000
## mu[34] 34.824 0.642 33.588 34.384 34.832 35.267 36.062 1.001 3000
## mu[35] 23.881 0.615 22.712 23.445 23.889 24.307 25.081 1.001 3000
## mu[36] 17.686 0.577 16.588 17.282 17.687 18.082 18.828 1.001 3000
## mu[37] 55.514 1.108 53.350 54.783 55.514 56.241 57.724 1.002 2400
## mu[38] 29.631 0.559 28.558 29.246 29.643 30.022 30.707 1.001 3000
## mu[39] 36.397 0.617 35.183 35.978 36.386 36.826 37.592 1.001 2700
## mu[40] 28.160 0.596 27.018 27.742 28.169 28.577 29.288 1.001 3000
## mu[41] 22.601 0.479 21.683 22.264 22.608 22.936 23.508 1.001 3000
## mu[42] 35.667 0.598 34.486 35.266 35.662 36.066 36.830 1.001 2500
## mu[43] 35.333 0.605 34.140 34.922 35.323 35.755 36.506 1.001 2900
## mu[44] 24.764 0.586 23.646 24.357 24.767 25.175 25.888 1.001 3000
## mu[45] 27.196 0.548 26.142 26.815 27.204 27.582 28.245 1.001 3000
## mu[46] 22.212 0.463 21.319 21.888 22.218 22.537 23.097 1.001 3000
## mu[47] 63.005 1.217 60.627 62.206 63.011 63.807 65.393 1.002 2300
## mu[48] 40.468 0.700 39.099 39.999 40.471 40.937 41.835 1.001 2200
## mu[49] 13.659 0.628 12.478 13.218 13.654 14.075 14.942 1.001 3000
## mu[50] 26.553 0.563 25.475 26.158 26.561 26.947 27.618 1.001 3000
## mu[51] 22.361 0.666 21.083 21.893 22.359 22.818 23.674 1.001 3000
## mu[52] 31.691 0.609 30.524 31.274 31.704 32.113 32.863 1.001 3000
## mu[53] 30.412 0.594 29.268 30.005 30.424 30.821 31.561 1.001 3000
## mu[54] 31.764 0.623 30.564 31.339 31.777 32.195 32.969 1.001 3000
## mu[55] 43.756 0.736 42.315 43.259 43.753 44.251 45.210 1.001 2300
## mu[56] 102.680 3.332 96.111 100.485 102.628 104.908 109.105 1.002 3000
## mu[57] 54.801 1.016 52.814 54.134 54.815 55.477 56.807 1.001 2200
## mu[58] 49.163 1.018 47.169 48.490 49.164 49.840 51.143 1.002 2400
## mu[59] 28.959 0.619 27.770 28.525 28.969 29.392 30.127 1.001 3000
## mu[60] 39.354 0.661 38.049 38.909 39.351 39.797 40.636 1.001 2500
## mu[61] 23.076 0.582 21.974 22.665 23.081 23.478 24.214 1.001 3000
## mu[62] 44.428 0.747 42.964 43.925 44.426 44.930 45.901 1.001 2300
## mu[63] 28.202 0.533 27.182 27.835 28.213 28.575 29.228 1.001 3000
## mu[64] 59.832 1.399 57.055 58.911 59.832 60.746 62.498 1.002 2700
## mu[65] 53.069 0.951 51.199 52.448 53.078 53.705 54.939 1.001 2200
## mu[66] 16.483 0.597 15.344 16.065 16.484 16.890 17.676 1.001 3000
## mu[67] 21.872 0.604 20.710 21.446 21.876 22.283 23.059 1.001 3000
## mu[68] 32.071 0.570 30.963 31.680 32.067 32.467 33.164 1.001 3000
## mu[69] 33.330 0.565 32.231 32.952 33.328 33.713 34.448 1.001 2200
## mu[70] 48.787 1.060 46.707 48.086 48.786 49.486 50.820 1.002 2600
## mu[71] 32.125 0.595 30.975 31.718 32.133 32.533 33.272 1.001 3000
## mu[72] 71.255 1.717 67.851 70.117 71.240 72.382 74.499 1.002 2800
## mu[73] 85.600 3.096 79.567 83.563 85.550 87.675 91.575 1.002 3000
## mu[74] 23.513 0.609 22.357 23.082 23.520 23.933 24.697 1.001 3000
## mu[75] 30.654 0.561 29.574 30.272 30.658 31.042 31.731 1.001 3000
## mu[76] 54.885 1.072 52.792 54.173 54.887 55.588 56.979 1.002 2300
## mu[77] 41.554 0.701 40.187 41.083 41.551 42.026 42.942 1.001 2300
## mu[78] 27.197 0.573 26.096 26.794 27.205 27.599 28.284 1.001 3000
## mu[79] 19.446 0.604 18.289 19.023 19.443 19.862 20.640 1.001 3000
## mu[80] 22.863 0.500 21.906 22.511 22.871 23.215 23.810 1.001 3000
## mu[81] 32.270 0.558 31.168 31.889 32.264 32.658 33.347 1.001 3000
## mu[82] 18.947 0.553 17.882 18.559 18.947 19.326 20.032 1.001 3000
## mu[83] 39.034 0.665 37.729 38.581 39.022 39.499 40.324 1.001 2800
## mu[84] 36.383 0.632 35.139 35.951 36.371 36.823 37.601 1.001 3000
## mu[85] 39.140 0.670 37.835 38.689 39.144 39.590 40.450 1.001 2200
## mu[86] 21.197 0.681 19.899 20.720 21.197 21.666 22.551 1.001 3000
## mu[87] 43.299 0.733 41.859 42.804 43.291 43.807 44.716 1.001 2700
## mu[88] 21.854 0.494 20.903 21.507 21.864 22.200 22.795 1.001 3000
## mu[89] 34.862 0.636 33.636 34.427 34.865 35.303 36.083 1.001 3000
## mu[90] 58.468 1.252 56.012 57.636 58.475 59.293 60.899 1.002 2500
## mu[91] 36.985 0.658 35.707 36.535 36.981 37.442 38.246 1.001 3000
## mu[92] 26.173 0.570 25.085 25.772 26.181 26.572 27.252 1.001 3000
## mu[93] 23.001 0.549 21.959 22.620 23.005 23.386 24.061 1.001 3000
## mu[94] 48.344 0.924 46.544 47.735 48.352 48.955 50.152 1.002 2300
## mu[95] 17.965 0.487 17.027 17.622 17.970 18.297 18.919 1.001 3000
## mu[96] 41.331 0.730 39.898 40.851 41.334 41.820 42.760 1.001 2200
## mu[97] 31.265 0.544 30.195 30.892 31.254 31.644 32.311 1.001 3000
## mu[98] 27.013 0.582 25.900 26.606 27.021 27.420 28.113 1.001 3000
## mu[99] 89.064 2.717 83.738 87.259 89.006 90.877 94.281 1.002 3000
## mu[100] 39.683 0.666 38.373 39.234 39.682 40.124 40.997 1.001 2300
## tau.nu 50.473 5.867 39.567 46.396 50.307 54.373 62.287 1.001 3000
## deviance 490.063 6.286 479.573 485.503 489.436 494.060 503.546 1.001 2300
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 19.8 and DIC = 509.8
## DIC is an estimate of expected predictive error (lower deviance is better).- INLA code
library(INLA)
result<-inla(y~x+f(nu,model = "iid"),
family = "poisson",
data =data,
control.predictor = list(link = 1),
)
result$summary.fixed## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) 0.8937826 0.12784427 0.6405243 0.8944573 1.1432129 0.8944562
## x 0.5148324 0.02401246 0.4678586 0.5147395 0.5623329 0.5147397
## kld
## (Intercept) 1.788962e-08
## x 2.049244e-08result$summary.linear.predictor## mean sd 0.025quant 0.5quant 0.975quant mode
## Predictor.001 3.806746 0.09930834 3.598734 3.812229 3.986701 3.828387
## Predictor.002 3.140168 0.10929541 2.908863 3.146801 3.339568 3.164715
## Predictor.003 3.020573 0.10657320 2.807321 3.020727 3.233221 3.020764
## Predictor.004 3.156025 0.10578196 2.955292 3.151934 3.375492 3.135120
## Predictor.005 3.363703 0.10313752 3.147584 3.368959 3.555523 3.385740
## Predictor.006 4.133310 0.08925268 3.949997 4.136503 4.301232 4.150745
## Predictor.007 2.634538 0.12037184 2.379290 2.641721 2.854290 2.661014
## Predictor.008 3.263434 0.10170626 3.061012 3.263100 3.467382 3.263067
## Predictor.009 3.405521 0.09982746 3.200683 3.408259 3.598432 3.408850
## Predictor.010 2.724181 0.11352650 2.491991 2.726654 2.944818 2.726943
## Predictor.011 3.809706 0.09184671 3.629095 3.808688 3.994664 3.808477
## Predictor.012 3.821166 0.09187044 3.642055 3.819530 4.007368 3.819182
## Predictor.013 3.356152 0.10014315 3.152942 3.357680 3.552738 3.357940
## Predictor.014 4.322509 0.08322823 4.161768 4.320976 4.490396 4.320662
## Predictor.015 2.780334 0.11587511 2.563977 2.774936 3.022801 2.756575
## Predictor.016 3.334975 0.10123143 3.126722 3.337970 3.530043 3.338601
## Predictor.017 3.547026 0.09761851 3.346200 3.550091 3.734429 3.565619
## Predictor.018 2.915799 0.10984110 2.688823 2.919371 3.126383 2.919988
## Predictor.019 4.073744 0.09005545 3.888554 4.077021 4.243344 4.091423
## Predictor.020 3.196062 0.10865219 2.965964 3.202795 3.393576 3.220681
## Predictor.021 2.993212 0.11571781 2.746120 3.001329 3.199332 3.020480
## Predictor.022 4.400163 0.08639015 4.239261 4.397392 4.576151 4.397041
## Predictor.023 2.707959 0.11596217 2.465165 2.713275 2.925310 2.731620
## Predictor.024 2.184472 0.12771942 1.922682 2.187601 2.430651 2.187850
## Predictor.025 3.316989 0.10065130 3.115143 3.317353 3.517264 3.317409
## Predictor.026 3.770511 0.09729780 3.567381 3.775458 3.949744 3.791346
## Predictor.027 4.098116 0.08733827 3.929867 4.096093 4.275471 4.095617
## Predictor.028 2.718902 0.11496338 2.498439 2.715439 2.955912 2.714997
## Predictor.029 2.717242 0.11333403 2.490933 2.717265 2.943496 2.717277
## Predictor.030 2.567989 0.11763318 2.325966 2.571240 2.794484 2.571623
## Predictor.031 3.625287 0.09810390 3.442684 3.620646 3.829468 3.604632
## Predictor.032 2.804132 0.11139708 2.578500 2.805545 3.023325 2.805715
## Predictor.033 3.486284 0.10928862 3.291985 3.479713 3.716003 3.462068
## Predictor.034 3.469590 0.09784343 3.275142 3.469196 3.665784 3.469142
## Predictor.035 3.061056 0.10854926 2.833471 3.066391 3.263813 3.083832
## Predictor.036 2.903399 0.10926173 2.687664 2.902340 3.124007 2.902228
## Predictor.037 4.068665 0.08999755 3.900387 4.065125 4.253962 4.050623
## Predictor.038 3.374190 0.10214429 3.161106 3.378870 3.566005 3.395376
## Predictor.039 3.622843 0.09509041 3.434990 3.622025 3.814167 3.621869
## Predictor.040 3.284617 0.10195474 3.075493 3.287277 3.482040 3.287787
## Predictor.041 3.304370 0.10097410 3.100882 3.305216 3.504257 3.305361
## Predictor.042 3.632261 0.09554180 3.437495 3.634372 3.817934 3.634830
## Predictor.043 3.633117 0.09616392 3.449672 3.629811 3.831211 3.614425
## Predictor.044 3.133636 0.10779573 2.906901 3.139393 3.333392 3.156858
## Predictor.045 3.351625 0.10015151 3.149173 3.352785 3.549128 3.352988
## Predictor.046 3.274840 0.10348331 3.059423 3.279223 3.470432 3.295812
## Predictor.047 3.996144 0.08928415 3.824496 3.993837 4.178099 3.993275
## Predictor.048 3.808451 0.09184859 3.626561 3.807993 3.992328 3.807910
## Predictor.049 2.513378 0.11995676 2.282565 2.510179 2.759852 2.509870
## Predictor.050 3.314636 0.10071256 3.113354 3.314694 3.515735 3.314717
## Predictor.051 3.148805 0.12065752 2.936451 3.140968 3.403854 3.121470
## Predictor.052 3.377695 0.10045377 3.171178 3.380640 3.571326 3.381288
## Predictor.053 3.336289 0.10277981 3.121867 3.340972 3.529437 3.357557
## Predictor.054 3.352426 0.10111043 3.144005 3.355669 3.546573 3.371711
## Predictor.055 3.748317 0.09308456 3.566634 3.746678 3.937032 3.746323
## Predictor.056 4.507484 0.07988345 4.353070 4.506254 4.667870 4.506040
## Predictor.057 3.937165 0.08995977 3.762582 3.935361 4.119667 3.934957
## Predictor.058 4.023169 0.08816857 3.847700 4.023312 4.198030 4.023342
## Predictor.059 3.395452 0.10398576 3.204475 3.389759 3.613625 3.372700
## Predictor.060 3.746099 0.09709811 3.567275 3.741254 3.948323 3.725348
## Predictor.061 3.128188 0.10460738 2.915255 3.130026 3.332982 3.130303
## Predictor.062 3.828446 0.09722175 3.651123 3.823500 4.030839 3.807678
## Predictor.063 3.447052 0.09827591 3.252811 3.446166 3.645093 3.446013
## Predictor.064 4.023838 0.09584266 3.824734 4.028454 4.198614 4.043757
## Predictor.065 3.956640 0.09341587 3.784272 3.952355 4.150163 3.937193
## Predictor.066 2.719431 0.11499182 2.479926 2.724046 2.937010 2.742138
## Predictor.067 3.184635 0.11277230 2.982636 3.177486 3.423054 3.158906
## Predictor.068 3.438093 0.10198923 3.224418 3.443341 3.627409 3.459981
## Predictor.069 3.684640 0.09701902 3.482905 3.689014 3.866239 3.704760
## Predictor.070 4.013327 0.09052147 3.827437 4.016419 4.184870 4.030880
## Predictor.071 3.475603 0.09785467 3.283010 3.474390 3.673460 3.474186
## Predictor.072 4.246564 0.08627777 4.083990 4.243720 4.422952 4.229982
## Predictor.073 4.412879 0.09482332 4.220756 4.415385 4.589514 4.415206
## Predictor.074 3.272302 0.11501839 3.068813 3.264950 3.515185 3.246237
## Predictor.075 3.418092 0.10082841 3.208533 3.422301 3.608730 3.438512
## Predictor.076 3.980133 0.08891771 3.805344 3.979276 4.158664 3.979106
## Predictor.077 3.672653 0.09698618 3.471113 3.676927 3.854622 3.692638
## Predictor.078 3.357288 0.10076289 3.162931 3.354456 3.564170 3.353875
## Predictor.079 3.032094 0.11204496 2.826653 3.025789 3.268111 3.007629
## Predictor.080 3.241777 0.10323333 3.028683 3.245170 3.439807 3.261525
## Predictor.081 3.499013 0.09989471 3.290923 3.503523 3.686611 3.519686
## Predictor.082 2.963475 0.10886237 2.738355 2.967101 3.171990 2.984224
## Predictor.083 3.632443 0.09557796 3.447900 3.629920 3.827989 3.629339
## Predictor.084 3.591707 0.09641365 3.406126 3.588932 3.789396 3.588260
## Predictor.085 3.803222 0.09213098 3.623209 3.801724 3.989716 3.801413
## Predictor.086 2.874287 0.11027750 2.648042 2.877050 3.087892 2.877451
## Predictor.087 3.636933 0.09503087 3.450821 3.635444 3.829465 3.635156
## Predictor.088 3.175841 0.10717157 2.950237 3.181688 3.373998 3.199096
## Predictor.089 3.488354 0.09754231 3.295425 3.487559 3.684742 3.487433
## Predictor.090 4.094071 0.08729403 3.925187 4.092293 4.270929 4.091896
## Predictor.091 3.585531 0.09862916 3.401445 3.580938 3.790808 3.564880
## Predictor.092 3.244986 0.10323148 3.031751 3.248458 3.442790 3.264830
## Predictor.093 3.106600 0.11037326 2.872655 3.113468 3.307271 3.131581
## Predictor.094 3.903795 0.09120531 3.717817 3.906075 4.079681 3.906608
## Predictor.095 3.001588 0.11077743 2.768061 3.007683 3.206192 3.025579
## Predictor.096 3.758190 0.09685323 3.556292 3.762927 3.937579 3.778717
## Predictor.097 3.462734 0.10253094 3.247226 3.468431 3.651239 3.485269
## Predictor.098 3.321899 0.10093028 3.124639 3.320043 3.527321 3.319733
## Predictor.099 4.463055 0.08130344 4.307664 4.461224 4.627430 4.460873
## Predictor.100 3.797335 0.09620548 3.620071 3.792603 3.997512 3.776866
## kld
## Predictor.001 5.638319e-07
## Predictor.002 5.149683e-07
## Predictor.003 2.152302e-07
## Predictor.004 2.824068e-07
## Predictor.005 3.761117e-07
## Predictor.006 2.084967e-07
## Predictor.007 4.940204e-07
## Predictor.008 2.148094e-07
## Predictor.009 2.417207e-07
## Predictor.010 2.228468e-07
## Predictor.011 1.642295e-07
## Predictor.012 1.707932e-07
## Predictor.013 2.217008e-07
## Predictor.014 8.961066e-08
## Predictor.015 3.361003e-07
## Predictor.016 2.514140e-07
## Predictor.017 2.431310e-07
## Predictor.018 2.597599e-07
## Predictor.019 2.153577e-07
## Predictor.020 5.475417e-07
## Predictor.021 8.050249e-07
## Predictor.022 2.600913e-07
## Predictor.023 3.231814e-07
## Predictor.024 1.966722e-07
## Predictor.025 2.139490e-07
## Predictor.026 3.938251e-07
## Predictor.027 1.358287e-07
## Predictor.028 2.367303e-07
## Predictor.029 2.011269e-07
## Predictor.030 2.279240e-07
## Predictor.031 3.326478e-07
## Predictor.032 2.139720e-07
## Predictor.033 9.590480e-07
## Predictor.034 2.039340e-07
## Predictor.035 3.572937e-07
## Predictor.036 2.118528e-07
## Predictor.037 2.433265e-07
## Predictor.038 3.289411e-07
## Predictor.039 1.890706e-07
## Predictor.040 2.438769e-07
## Predictor.041 2.174160e-07
## Predictor.042 2.091068e-07
## Predictor.043 2.391120e-07
## Predictor.044 4.016580e-07
## Predictor.045 2.180079e-07
## Predictor.046 3.079978e-07
## Predictor.047 1.607996e-07
## Predictor.048 1.621625e-07
## Predictor.049 2.155101e-07
## Predictor.050 2.134591e-07
## Predictor.051 1.271453e-06
## Predictor.052 2.484744e-07
## Predictor.053 3.278199e-07
## Predictor.054 2.583005e-07
## Predictor.055 1.802973e-07
## Predictor.056 6.000168e-08
## Predictor.057 1.559340e-07
## Predictor.058 1.242576e-07
## Predictor.059 4.425330e-07
## Predictor.060 3.856250e-07
## Predictor.061 2.298004e-07
## Predictor.062 4.747693e-07
## Predictor.063 2.064258e-07
## Predictor.064 4.923199e-07
## Predictor.065 3.489422e-07
## Predictor.066 2.859692e-07
## Predictor.067 7.470390e-07
## Predictor.068 3.820469e-07
## Predictor.069 3.108737e-07
## Predictor.070 2.020055e-07
## Predictor.071 2.069318e-07
## Predictor.072 1.770757e-07
## Predictor.073 3.715476e-07
## Predictor.074 1.067277e-06
## Predictor.075 3.000122e-07
## Predictor.076 1.350124e-07
## Predictor.077 3.028136e-07
## Predictor.078 2.386168e-07
## Predictor.079 4.651779e-07
## Predictor.080 2.644018e-07
## Predictor.081 3.199934e-07
## Predictor.082 2.640646e-07
## Predictor.083 2.128496e-07
## Predictor.084 2.222150e-07
## Predictor.085 1.714307e-07
## Predictor.086 2.371461e-07
## Predictor.087 1.936766e-07
## Predictor.088 4.162807e-07
## Predictor.089 2.034492e-07
## Predictor.090 1.304434e-07
## Predictor.091 3.253798e-07
## Predictor.092 2.670577e-07
## Predictor.093 5.464742e-07
## Predictor.094 1.774537e-07
## Predictor.095 4.195147e-07
## Predictor.096 3.574937e-07
## Predictor.097 4.419375e-07
## Predictor.098 2.219390e-07
## Predictor.099 9.150100e-08
## Predictor.100 3.779391e-07plot(result$marginals.fixed$x,
type = "l",xlab= "beta")
我没有展示这些分布，但参数估计是可以比较的。再次强调，节省的时间真的很显著。在我的博士期间，我提出了一种多元（多个结果）模型用于预测，如果没有INLA方法，这个模型在我所工作的公共卫生应用中将不具备实际可行性，因为模型必须每天重新拟合。我会在这里引用我关于该模型的论文。虽然我现在的工作倾向于复杂性较低，但我仍然喜欢在许多贝叶斯工作中使用INLA。
我唯一偏好MCMC代码的地方，除了更直观的理论，就是BUGS/JAGS与INLA的语法。虽然INLA看起来不错，因为它类似于R中的glm函数，但有时先验分布对我来说显得过于隐蔽。我发现BUGS代码很好地突出了模型的层次结构，这对教学和学习很有帮助。