R语言Stan贝叶斯空间条件自回归CAR模型分析死亡率多维度数据可视化

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

在空间数据分析领域，准确的模型和有效的工具对于研究人员至关重要。本文为区域数据的贝叶斯模型分析提供了一套完整的工作流程，基于Stan这一先进的贝叶斯建模平台构建，帮助客户为空间分析带来了新的解决方案（点击文末“阅读原文”获取完整代码、数据、文档）。

特性

（一）空间回归与疾病映射

本文处理区域单元（如州、县或普查区域）或网络记录的数据统计模型，其中包括空间计量经济学模型。在疾病映射方面，可用于估计小区域（如县）的疾病风险，分析健康结果与其他区域变量的协变关系。例如在公共卫生研究中，研究人员可以利用该功能探究不同地区疾病的分布情况，分析发病率与环境、人口密度等因素的关联。

（二）空间分析工具

本文提供了可视化和测量空间自相关及地图模式的工具，用于探索性分析和模型诊断。在进行空间数据分析时，这些工具可以帮助研究人员直观地了解数据的空间分布特征，判断数据是否存在自相关现象，为后续的建模和分析提供依据。

（三）观测不确定性处理

本文可以将数据可靠性信息，比如社区调查估计的标准误差，纳入到任何模型中。这一特性使得模型能够更好地处理数据的不确定性，提高模型的准确性和可靠性。在实际应用中，数据往往存在各种误差和不确定性，这一功能能够有效考虑这些因素，使分析结果更符合实际情况。

（四）缺失和审查观测值处理

在生命统计和疾病监测系统中，会对低于阈值数量的病例数进行审查。能够对存在审查观测值或缺失观测值的小区域疾病或死亡率风险进行建模。例如在分析某地区的死亡率数据时，可能存在部分数据因各种原因被审查或缺失，可以通过合适的模型对这些数据进行处理，从而得到更准确的死亡率估计。

（五）与RStan生态系统的交互

本文能与许多高质量的R语言贝叶斯建模包轻松交互，这为研究人员提供了更丰富的分析手段和更灵活的建模选择。研究人员可以结合其他相关包的优势，进行更复杂和深入的数据分析。

（六）自定义空间模型

本文提供了在Stan中构建自定义空间或网络模型的工具，满足研究人员对于特定研究问题的个性化建模需求。例如在研究特定地理区域的交通流量网络或生态系统中的物种分布时，研究人员可以根据实际情况构建自定义模型。
在公共卫生研究中，与surveil这个R语言包形成互补，两者结合可以从时空两个维度全面分析公共卫生数据。

示例

（一）数据加载与初步分析

数据集包含了2014 - 2018年期间，55 - 64岁年龄段按性别划分的县人口和死亡率数据。由于公共访问数据的常见情况，部分观测值因疾病控制与预防中心为保护隐私而被审查，从而缺失。
对空间数据进行可视化总结，包括直方图、莫兰散点图和地图。莫兰散点图用于展示数据值与其相邻值的汇总对比，回归线的斜率可衡量自相关程度。
以下代码用于创建邻接矩阵，计算每10,000人的粗死亡率，并进行快速空间诊断：

# 快速空间诊断
diag(mortaate, georgia, w = adjacenix")

执行上述代码后，得到的可视化结果如下：

点击标题查阅往期内容

R语言和STAN,JAGS：用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

左右滑动查看更多

（二）空间条件自回归（CAR）模型拟合

由于县的死亡率和其他健康统计数据在许多情况下是高度不稳定的估计值，不能直接用于公共建议或推断（因为人口规模较小），因此需要使用模型从小区域数据中进行推断。
以下代码使用空间条件自回归（CAR）模型对女性县死亡率数据进行拟合：

# 对女性死亡率数据拟合模型
stacar(deaths.female ~ offse

将拟合好的模型传递给函数，可返回空间模型的一组诊断信息：

使用print方法可以返回模型参数的概率分布摘要，以及来自Stan的马尔可夫链蒙特卡罗（MCMC）诊断信息（均值的蒙特卡罗标准误差se_mean、有效样本大小n_eff和R-hat统计量Rhat）：

print(modefit)

执行sp_g(mode_i, grgia)后，得到的诊断信息可视化结果如下：

（三）提取死亡率估计值与可视化

通过fitted方法提取县死亡率估计值，乘以10,000得到每10,000人的死亡率：

# 每10,000人的死亡率估计值
moraitystimtes <- ftted(de_fit) * 10e3

将估计值放入分箱中用于地图颜色显示，创建地图展示估计值：

oriial_magin = par(mar = rep(1, 4))
# 获取边界
geometry <- sf:st_eoetry(eorgia)
# 绘制地图
plot(geomty,lwd = 0.2,col = colors)
# 添加图例
legbty = 'n')

执行上述代码后，绘制出的地图如下：

利用可信区间，创建点区间图：

# 按死亡率对县进行排序
data <- mortality_estimates\[index, \]
# 收集估计值和95%可信区间
estimate <- data$mean
lower <- data$\`2.5%\`
upper <- data$\`97.5%\`
y <- seq\_along(cony\_name)
x_limit <- (min(lo), max(upper)) %>% round()
# 设置边距
originalarg = r(ma  c(3, , 0)
# 绘制点
plot(esimte,y,pch = 5,col = gray50',bty  'L,axes = FALSE,xlim = _limit,yl = A,xlab = NA)
# 绘制区间
segmens(x0 = lower, x1 = upper,y0 = y,y1  ,col = olor\[index\])
# 添加x轴mit\[1\], x_limit\[2\], by = 20))

执行上述代码后，得到的点区间图如下：