【原创】R语言使用贝叶斯 层次模型进行空间数据分析报告论文(附代码数据)
咨询QQ：3025393450
有问题百度搜索“”就可以了
欢迎登陆官网：https://www.360docs.net/doc/b017561476.html,/datablog
R语言使用贝叶斯层次模型进行空间数据分析报告
介绍
在本节中，我将重点介绍使用集成嵌套拉普拉斯近似方法的贝叶斯推理。
可以估计贝叶斯层次模型的后边缘分布。鉴于模型类型非常广泛，我们将重点关注用于分析晶格数据的空间模型。
数据集：纽约州北部的白血病
为了说明如何与空间模型拟合，将使用纽约白血病数据集。该数据集记录了普查区纽约州北部的许多白血病病例。数据集中的一些变量是：
?Cases：1978-1982年期间的白血病病例数。
?POP8：1980年人口。
?PCTOWNHOME：拥有房屋的人口比例。
?PCTAGE65P：65岁以上的人口比例。
?AVGIDIST：到最近的三氯乙烯（TCE）站点的平均反距离。
可以按以下方式加载数据集：
library(spdep)
library(DClusterm)
咨询QQ：3025393450
有问题百度搜索“”就可以了
欢迎登陆官网：https://www.360docs.net/doc/b017561476.html,/datablog
data(NY8)
鉴于有兴趣研究纽约州北部的白血病风险，因此首先要计算预期的病例数。这是通过计算总死亡率（总病例数除以总人口数）并将其乘以总人口数得出的：
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))
咨询QQ：3025393450
有问题百度搜索“”就可以了
欢迎登陆官网：https://www.360docs.net/doc/b017561476.html,/datablog
## Loading required package: viridisLite
可以使用以下tmap软件包轻松创建交互式地图：
library(tmap)
tmap_mode("view")
SMR_map <- tm_shape(NY8[syracuse, ]) +
tm_fill(col = "SMR", alpha = 0.35) +
tm_borders() +
tm_shape(TCE) + tm_dots(col = "red") # Add TCE sites widgetframe::frameWidget(print(SMR_map))
咨询QQ：3025393450
有问题百度搜索“”就可以了
欢迎登陆官网：https://www.360docs.net/doc/b017561476.html,/datablog
请注意，先前的地图还包括11个受TCE污染的站点的位置，可以通过缩小看到它。
混合效应模型
泊松回归
我们将考虑的第一个模型是没有潜在随机效应的Poisson模型，因为这将提供与其他模型进行比较的基准。
模型：
library(INLA)
m1 <- inla(Cases ~ 1 + AVGIDIST,
data = as.data.frame(NY8),
family = "poisson",
E = NY8$Expected, control.predictor = list(compute = TRUE),
https://www.360docs.net/doc/b017561476.html,pute = list(dic = TRUE, waic = TRUE))
请注意，它的glm功能类似于该功能。在此，参数E用于预期的案例数。或设置了其他参数来计算模型参数的边际
（使用control.predictor）并计算一些模型选择标准（使用https://www.360docs.net/doc/b017561476.html,pute）。
接下来，可以获得模型的摘要：
summary(m1)
##
咨询QQ：3025393450
有问题百度搜索“”就可以了
欢迎登陆官网：https://www.360docs.net/doc/b017561476.html,/datablog
## Call:
## c("inla(formula = Cases ~ 1 + AVGIDIST, family = \"poisson\", data
## = as.data.frame(NY8), ", " E = NY8$Expected, https://www.360docs.net/doc/b017561476.html,pute =
## list(dic = TRUE, waic = TRUE), ", " control.predictor =
## list(compute = TRUE))")
## 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高斯随机效应，将潜在随机效应添加到模型中，以解决过度分散问题，
是通过f-function 指定的INLA：
咨询QQ：3025393450
有问题百度搜索“”就可以了
欢迎登陆官网：https://www.360docs.net/doc/b017561476.html,/datablog
NY8$ID <- 1:nrow(NY8)
m2 <- inla(Cases ~ 1 + AVGIDIST + f(ID, model = "iid"),
data = as.data.frame(NY8), family = "poisson",
E = NY8$Expected,
control.predictor = list(compute = TRUE),
https://www.360docs.net/doc/b017561476.html,pute = list(dic = TRUE, waic = TRUE))
现在，该模式的摘要包括有关随机效果的信息：
summary(m2)
##
## Call:
## c("inla(formula = Cases ~ 1 + AVGIDIST + f(ID, model = \"iid\"),
## family = \"poisson\", ", " data = as.data.frame(NY8), E =
## NY8$Expected, https://www.360docs.net/doc/b017561476.html,pute = list(dic = TRUE, ", " waic =
## TRUE), control.predictor = list(compute = TRUE))" )
## 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.347 0.105 0.139 0.346 0.558 0.344 0
##
## Random effects:
## Name Model
## ID IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode