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、时空分析野生动物数据可视化

左右滑动查看更多

最后，绘制θ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 <- 0
 for(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 = 0
 for(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则在处理空间效应方面具有优势，能够有效地捕捉疾病发病率的空间相关性。在实际研究中，研究者可以根据数据特点和研究目的选择合适的模型，以更准确地分析环境暴露与疾病发病之间的关联，为公共卫生决策提供有力的支持。