R语言使用贝叶斯层次模型进行空间数据分析

 介绍

最近我们被客户要求撰写关于贝叶斯层次模型的研究报告，包括一些图形和统计输出。在本节中，我将重点介绍使用集成嵌套 拉普拉斯近似方法的贝叶斯推理。 
可以 估计贝叶斯 层次模型的后边缘分布。 鉴于模型类型非常广泛，我们将重点关注用于分析晶格数据的空间模型。

数据集：纽约州北部的白血病

为了说明如何与空间模型拟合，将使用纽约白血病数据集。该数据集记录了普查区纽约州北部的许多白血病病例。数据集中的一些变量是：

• Cases：1978-1982年期间的白血病病例数。
• POP8：1980年人口。
• PCTOWNHOME：拥有房屋的人口比例。
• PCTAGE65P：65岁以上的人口比例。
• AVGIDIST：到最近的三氯乙烯（TCE）站点的平均反距离。

鉴于有兴趣研究纽约州北部的白血病风险，因此首先要计算预期的病例数。这是通过计算总死亡率（总病例数除以总人口数）并将其乘以总人口数得出的：

rate <- sum(NY8$Cases) / sum(NY8$POP8)
NY8$Expected <- NY8$POP8 * rate

一旦获得了预期的病例数，就可以使用标准化死亡率（SMR）来获得原始的风险估计，该标准是将观察到的病例数除以预期的病例数得出的：

NY8$SMR <- NY8$Cases / NY8$Expected

疾病作图

在流行病学中，重要的是制作地图以显示相对风险的空间分布。在此示例中，我们将重点放在锡拉库扎市以减少生成地图的计算时间。因此，我们用锡拉丘兹市的区域创建索引：

# Subset Syracuse city
syracuse <- which(NY8$AREANAME == "Syracuse city")

可以使用函数spplot（在包中sp）简单地创建疾病图：

library(viridis)
## Loading required package: viridisLite
spplot(NY8[syracuse, ], "SMR", #at = c(0.6, 0.9801, 1.055, 1.087, 1.125, 13),
   col.regions = rev(magma(16))) #gray.colors(16, 0.9, 0.4))

## Loading required package: viridisLite

可以轻松创建交互式地图

请注意，先前的地图还包括11个受TCE污染的站点的位置，可以通过缩小看到它。

混合效应模型

泊松回归

我们将考虑的第一个模型是没有潜在随机效应的Poisson模型，因为这将提供与其他模型进行比较的基准。
 

 模型 ：

请注意，它的glm功能类似于该功能。在此，参数 E用于预期的案例数。或  设置了其他参数来计算模型参数的边际
（使用control.predictor）并计算一些模型选择标准 （使用control.compute）。

接下来，可以获得模型的摘要：

summary(m1)

## 
## Call:

## Time used:
##     Pre = 0.368, Running = 0.0968, Post = 0.0587, Total = 0.524 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.065 0.045     -0.155   -0.065      0.023 -0.064   0
## AVGIDIST     0.320 0.078      0.160    0.322      0.465  0.327   0
## 
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 140.25 
## 
## Deviance Information Criterion (DIC) ...............: 948.12
## Deviance Information Criterion (DIC, saturated) ....: 418.75
## Effective number of parameters .....................: 2.00
## 
## Watanabe-Akaike information criterion (WAIC) ...: 949.03
## Effective number of parameters .................: 2.67
## 
## Marginal log-Likelihood:  -480.28 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

具有随机效应的泊松回归

可以通过 在线性预测变量中包括iid高斯随机效应，将潜在随机效应添加到模型中，以解决过度分散问题。

现在，该模式的摘要包括有关随机效果的信息：

summary(m2)

## 
## Call:

## Time used:
##     Pre = 0.236, Running = 0.315, Post = 0.0744, Total = 0.625 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.126 0.064     -0.256   -0.125     -0.006 -0.122   0
## AVGIDIST     0