R语言NIMBLE、Stan和INLA贝叶斯平滑及条件空间模型死亡率数据分析：提升疾病风险估计准确性...

全文链接：https://tecdat.cn/?p=40365

在环境流行病学研究中，理解空间数据的特性以及如何通过合适的模型分析疾病的空间分布是至关重要的。本文主要介绍了不同类型的空间数据、空间格点过程的理论，并引入了疾病映射以及对空间风险进行平滑处理的模型。通过这些内容，读者可以深入了解如何通过借用相邻区域的信息来改进风险估计，以及如何运用经验贝叶斯或全贝叶斯方法进行平滑处理等关键问题（点击文末“阅读原文”获取完整代码、数据、文档）。

不同模型下的疾病风险分析

经验贝叶斯和贝叶斯平滑在COPD死亡率分析中的应用

以2010年慢性阻塞性肺疾病（COPD）住院情况为例，对该疾病的发病风险进行分析。有Nl=324Nl=324个地方行政区，每个行政区都有观察到的病例数YlYl和预期病例数ElEl，l=1,...,324l=1,...,324。预期病例数是通过将整体的年龄 - 性别特定发病率应用到每个地区的年龄 - 性别人口结构上，采用间接标准化方法计算得出的。
为了进行经验贝叶斯平滑。以下是相应的R代码：

# 观察数据
Y\_obs <- obs\_data$Y2010
# 偏移量
E\_offset <- exp\_data$E2010
RR\_vals <- eBayes(Y\_obs, E_offset)
plot(RR\_vals$SMR, RR\_vals$RR, xlim = c(0, 2

从平滑结果可以看出，较低和较高的标准化死亡率比（SMRs）都向总体平均值μμ靠近，在这个例子中μ=exp(−0.0309)=0.9696。由于这些地区相对较大且人口众多，平滑效果相对有限。在这个例子中，αα估计值为14.6。
对于全贝叶斯分析，可以使用NIMBLE。

a ~ dgamma(1, 1)beta0 ~ dnorm(0, 10)
})
# 观察数据
y\_obs <- obs\_data$Y2010
# 偏移量
E\_offset <- exp\_data$E2010
N <- length(y_obs)
# 常数列表
constants_list <- list(N = N,E = E_offset
)

通过检查所有可用样本中的最小值来检验有效样本量，以确保该最小值是可接受的。同时，还提供了一些参数链的轨迹图。运行上述代码后，会得到有效样本量和WAIC的相关信息。接下来，对后验分布的样本进行处理，并绘制一些参数的轨迹图：

# 绘制theta的轨迹图
for (i in 1:3)plot(mvSamples\[, c(paste("theta\[", i, "\]", sep = ""))\], bty = "n")

在检查链的收敛性后，可以绘制每个参数的后验均值和95%置信区间：

for (i in 1:N)segments(y\_obs\[i\], post\_fitted$\`95%CI\_low\`\[i\], y\_obs\[i\], post\_fitted$\`95%CI\_upp\`\[i\])
abline(a = 0, b = 1)

Stan模型分析

在Stan中实现相同的模型，首先需要加载必要的包并读取数据：

接下来，定义数据相关的对象并在Stan中运行模型：

# 观察数据
y\_obs <- obs\_data$Y2010
# 偏移量
E\_offset <- exp\_data$E2010
N <- length(y_obs)
# 数据列表
data_list <- list(N = length(y_obs),Y = y_obs,E = E_offset
)

检查一些参数的轨迹图、有效样本量，并获取一些参数的后验摘要：

---

点击标题查阅往期内容

R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化

左右滑动查看更多

01

02

03

04

最后，绘制θiθi和拟合值的后验摘要及其95%后验可信区间：

for (i in 1:N)segments(y\_obs\[i\], summary\_fit$\`5%\`\[i\], y\_obs\[i\], summary\_fit$\`95%\`\[i\])
abline(a = 0, b = 1)

条件空间模型分析

NIMBLE中的条件空间模型拟合

使用NIMBLE拟合具有空间随机效应的Poisson对数正态模型，空间随机效应来自ICAR模型。首先加载必要的包和数据：

# 绘制平均剥夺分数与2010年SMR的散点图ggplot(copd\_data\_df) + geom_point(aes(x = Average.Score, y = SMR2010)) +theme_classic() +【【【【【【【【【【【【【】

为了在Nimble中拟合ICAR模型，需要获取邻域矩阵WW。首先定义一个函数来计算每个区域的邻居数量：

# 定义计算邻居列表的函数
adjlist <- function(W, N) {adj <- 0for(i in 1:N) {for(j in 1:N) {if(W\[i, j\] == 1) {adj <- append(adj, j)}}}adj <- adj\[-1\]return(adj)
}

定义Nimble模型代码

接下来，通过跟踪图和部分参数的有效样本量来再次检查链的收敛性：

截距、平均剥夺指数的固定效应以及ICAR先验分布精度的后验摘要：

潜在效应后验均值的地图：

ggplot() +# 选择空间对象和用于绘图的列

Stan中的条件空间模型拟合

由于Stan没有专门用于ICAR先验的函数：

# 定义计算邻居列表的函数
adjlist = function(W, N) {adj = 0for(i in 1:N) {for(j in 1:N) {if(W\[i, j\] == 1) {adj = append(adj, j)}}}adj = adj\[-1\]return(adj)
}

该模型存储在名为stan的单独文件中

定义邻接矩阵和相邻区域的索引集：

# 创建邻域
W\_nb <- poly2nb(england\_bound, row.names = rownames(england_bound))
# 创建邻接矩阵
W\_mat <- nb2mat(W\_nb, style = "B")
# 定义用于Stan的空间结构
N <- length(unique(england_bound$ID))
neigh <- adjlist(W_mat, N)
numneigh <- apply(W_mat, 2, sum)

部分参数后验样本的跟踪图：

接下来获取截距、与剥夺分数相关的固定效应以及随机效应标准差的后验摘要。riskd表示剥夺程度增加一个单位时的相对风险：

结果表明，剥夺程度增加一个单位，相对风险增加 2.2%。需要注意的是，ICAR先验的标准差指的是条件分布的标准差。
现在绘制潜在效应的后验均值：

ggplot() +# 选择空间对象和用于绘图的列geom\_sf(data = CARan\_merge, aes(fill = summary.mean)) + 
# 更改图例标签

使用CARBayes拟合条件模型

这里考虑使用R包CARBayes对呼吸道入院数据拟合具有空间效应的Poisson对数正态模型。

使用INLA拟合条件模型

现在使用R - INLA对呼吸道入院数据拟合具有空间效应的Poisson对数正态模型。R - INLA使用集成嵌套拉普拉斯近似来近似得到的后验分布。使用R - INLA拟合模型的调用与拟合glm非常相似，潜在效应通过函数f(.)引入。

结论

通过对不同模型（经验贝叶斯、NIMBLE、Stan、CARBayes、INLA）在疾病风险分析和条件空间模型拟合中的应用研究，我们可以看到每种模型都有其独特的优势和适用场景。经验贝叶斯方法能够对风险估计进行平滑处理，减少基于小样本数据估计的不稳定性；NIMBLE和Stan在处理复杂的贝叶斯模型时表现出色，能够对模型参数进行准确估计；CARBayes和INLA则在处理空间效应方面具有优势，能够有效地捕捉疾病发病率的空间相关性。在实际研究中，研究者可以根据数据特点和研究目的选择合适的模型，以更准确地分析环境暴露与疾病发病之间的关联，为公共卫生决策提供有力的支持。

本文中分析的完整数据、代码、文档分享到会员群，扫描下面二维码即可加群！ 

---

资料获取

在公众号后台回复“领资料”，可免费获取数据分析、机器学习、深度学习等学习资料。

点击文末“阅读原文”

获取完整代码、数据、文档。

本文选自《R语言NIMBLE、Stan和INLA贝叶斯平滑及条件空间模型死亡率数据分析：提升疾病风险估计准确性》。

点击标题查阅往期内容

R语言贝叶斯广义线性混合（多层次/水平/嵌套）模型GLMM、逻辑回归分析教育留级影响因素数据

R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例

用SPSS估计HLM多层（层次）线性模型模型

R语言用线性混合效应（多水平/层次/嵌套）模型分析声调高低与礼貌态度的关系

R语言LME4混合效应模型研究教师的受欢迎程度

R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model分析藻类数据实例

R语言混合线性模型、多层次模型、回归模型分析学生平均成绩GPA和可视化

R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例

R语言用lme4多层次（混合效应）广义线性模型（GLM），逻辑回归分析教育留级调查数据

R语言 线性混合效应模型实战案例

R语言混合效应逻辑回归（mixed effects logistic）模型分析肺癌数据

R语言如何用潜类别混合效应模型（LCMM）分析抑郁症状

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言建立和可视化混合效应模型mixed effect model

R语言LME4混合效应模型研究教师的受欢迎程度

R语言 线性混合效应模型实战案例

R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言如何解决线性混合模型中畸形拟合(Singular fit)的问题

基于R语言的lmer混合线性回归模型

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

R语言分层线性模型案例

R语言用WinBUGS 软件对学术能力测验（SAT）建立分层模型

使用SAS，Stata，HLM，R，SPSS和Mplus的分层线性模型HLM

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

SPSS中的多层（等级）线性模型Multilevel linear models研究整容手术数据

用SPSS估计HLM多层（层次）线性模型模型