关于算法: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.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
## Precision for ID 3712.34 11263.70       3.52     6.94   39903.61 5.18
## 
## Expected number of effective parameters(stdev): 54.95(30.20)
## Number of equivalent replicates : 5.11 
## 
## Deviance Information Criterion (DIC) ...............: 926.93
## Deviance Information Criterion (DIC, saturated) ....: 397.56
## Effective number of parameters .....................: 61.52
## 
## Watanabe-Akaike information criterion (WAIC) ...: 932.63
## Effective number of parameters .................: 57.92
## 
## Marginal log-Likelihood:  -478.93 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

这两个模型预计能够被增加到 SpatialPolygonsDataFrame NY8

NY8$FIXED.EFF <- m1$summary.fitted\[, "mean"\]
NY8$IID.EFF <- m2$summary.fitted\[, "mean"\]
spplot(NY8\[syracuse, \], c("SMR", "FIXED.EFF", "IID.EFF"),
  col.regions = rev(magma(16)))

格子数据波及在不同区域（例如，邻里，城市，省，州等）测量的数据。呈现空间依赖性是因为相邻区域将显示类似的指标变量值。

能够应用poly2nbpackage 中的函数来计算邻接矩阵 spdep。如果其边界至多在某一点上接触，则此性能会将两个区域视为街坊：

这将返回一个 nb 具备邻域构造定义的对象：

NY8.nb

## Neighbour list object:
## Number of regions: 281 
## Number of nonzero links: 1624 
## Percentage nonzero weights: 2.056712 
## Average number of links: 5.779359

另外，当多边形的重心已知时，能够绘制对象：

plot(NY8) 
plot(NY8.nb, coordinates(NY8), add = TRUE, pch = ".", col = "gray")

通常状况是，除了 \（y\_i \）之外，咱们还有许多协变量 \（X\_i \）。因而，咱们可能想对 \（X_i \）_回归_ \（y_i \）。除了协变量，咱们可能还须要思考数据的空间结构。
能够应用不同类型的回归模型来建模晶格数据：

狭义线性模型（具备空间随机效应）。
空间计量经济学模型。

一种常见的办法（对于高斯数据）是应用
具备随机效应的线性回归：

\ [
Y = X \ beta + Zu + \ varepsilon
\]

随机效应的向量 \（u \）被建模为多元正态分布：

\ [
u \ sim N（0，\ sigma ^ 2_u \ Sigma）
\]

\（\ Sigma \）的定义是，它会引起与相邻区域的更高相关性，\（Z \）是随机成果的设计矩阵，而
\（\ varepsilon_i \ sim N（0，\ sigma ^ 2），i = 1，\ ldots，n \）是一个误差项。

在 \（\ Sigma \）中包含空间依赖的办法有很多：

同步自回归（SAR）：

\ [
\ Sigma ^ {-1} = [（I- \ rho W）’（I- \ rho W）]
\]

条件自回归（CAR）：

\ [
\ Sigma ^ {-1} =（I- \ rho W）
\]

（ICAR）：

\ [
\ Sigma ^ {-1} = diag（n_i）– W
\]

\（n_i \）是区域 \（i \）的街坊数。
\（\ Sigma_ {i，j} \）取决于 \（d（i，j）\）的函数。例如：

\ [
\ Sigma_ {i，j} = \ exp \ {-d（i，j）/ \ phi \}
\]

矩阵的“混合”（Leroux 等人的模型）：

\ [
\ Sigma = [（1 – \ lambda）I_n + \ lambda M] ^ {-1}; \ \ lambda \ in（0,1）
\]

第一个示例将基于 ICAR 标准。请留神，应用f- 函数定义空间潜在成果。这将须要一个索引来辨认每个区域中的随机效应，模型的类型和邻接矩阵。为此，将应用稠密矩阵。

## 
## Call:

## Time used:
##     Pre = 0.298, Running = 0.305, Post = 0.0406, Total = 0.644 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.122 0.052     -0.226   -0.122     -0.022 -0.120   0
## AVGIDIST     0.318 0.121      0.075    0.320      0.551  0.324   0
## 
## Random effects:
##   Name     Model
##     ID Besags ICAR model
## 
## Model hyperparameters:
##                  mean   sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID 3.20 1.67       1.41     2.79       7.56 2.27
## 
## Expected number of effective parameters(stdev): 46.68(12.67)
## Number of equivalent replicates : 6.02 
## 
## Deviance Information Criterion (DIC) ...............: 904.12
## Deviance Information Criterion (DIC, saturated) ....: 374.75
## Effective number of parameters .....................: 48.31
## 
## Watanabe-Akaike information criterion (WAIC) ...: 906.77
## Effective number of parameters .................: 44.27
## 
## Marginal log-Likelihood:  -685.70 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Besag，York 和 Mollié模型包含两个潜在的随机效应：ICAR 潜在效应和高斯 iid 潜在效应。线性预测变量 \（\ eta_i \）
为：

\ [
\ eta\_i = \ alpha + \ beta AVGIDIST\_i + u\_i + v\_i
\]

\（u_i \）是 iid 高斯随机效应
\（v_i \）是外在的 CAR 随机效应

## 
## Call:

## Time used:
##     Pre = 0.294, Running = 1, Post = 0.0652, Total = 1.36 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.123 0.052     -0.227   -0.122     -0.023 -0.121   0
## AVGIDIST     0.318 0.121      0.075    0.320      0.551  0.324   0
## 
## Random effects:
##   Name     Model
##     ID BYM model
## 
## Model hyperparameters:
##                                         mean      sd 0.025quant 0.5quant
## Precision for ID (iid component)     1790.06 1769.02     115.76  1265.24
## Precision for ID (spatial component)    3.12    1.36       1.37     2.82
##                                      0.975quant   mode
## Precision for ID (iid component)        6522.28 310.73
## Precision for ID (spatial component)       6.58   2.33
## 
## Expected number of effective parameters(stdev): 47.66(12.79)
## Number of equivalent replicates : 5.90 
## 
## Deviance Information Criterion (DIC) ...............: 903.41
## Deviance Information Criterion (DIC, saturated) ....: 374.04
## Effective number of parameters .....................: 48.75
## 
## Watanabe-Akaike information criterion (WAIC) ...: 906.61
## Effective number of parameters .................: 45.04
## 
## Marginal log-Likelihood:  -425.65 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

该模型是应用矩阵的“混合”（Leroux 等人的模型）
定义的，以定义潜在效应的精度矩阵：

\ [
\ Sigma ^ {-1} = [（1-\ lambda）I_n + \ lambda M]; \ \ lambda \ in（0,1）
\]

为了定义正确的模型，咱们应采纳矩阵 \（C \）如下：

\ [
C = I\_n – M; \ M = diag（n\_i）– W
\]

而后，\（\ lambda_ {max} = 1 \）和

\ [
\ Sigma ^ {-1} =
\ frac {1} {\ tau}（I\_n- \ frac {\ rho} {\ lambda\_ {max}} C）=
\ frac {1} {\ tau}（I\_n- \ rho（I\_n – M））= \ frac {1} {\ tau}（（1- \ rho）I_n + \ rho M）
\]

为了拟合模型，第一步是创立矩阵 \（M \）：

咱们能够查看最大特征值 \（\ lambda_ {max} \）是一个：

max(eigen(Cmatrix)$values)
## \[1\] 1

## \[1\] 1

该模型与平常一样具备性能 inla。留神，\（C \）矩阵应用参数
传递给 f 函数Cmatrix：

## 
## Call:
## Time used:
##     Pre = 0.236, Running = 0.695, Post = 0.0493, Total = 0.98 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode   kld
## (Intercept) -0.128 0.448      -0.91   -0.128      0.656 -0.126 0.075
## AVGIDIST     0.325 0.122       0.08    0.327      0.561  0.330 0.000
## 
## Random effects:
##   Name     Model
##     ID Generic1 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ID 2.720 1.098       1.27    2.489      5.480 2.106
## Beta for ID      0.865 0.142       0.47    0.915      0.997 0.996
## 
## Expected number of effective parameters(stdev): 52.25(13.87)
## Number of equivalent replicates : 5.38 
## 
## Deviance Information Criterion (DIC) ...............: 903.14
## Deviance Information Criterion (DIC, saturated) ....: 373.77
## Effective number of parameters .....................: 53.12
## 
## Watanabe-Akaike information criterion (WAIC) ...: 906.20
## Effective number of parameters .................: 48.19
## 
## Marginal log-Likelihood:  -474.94 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

空间计量经济学是通过对空间建模略有不同的办法开发的。除了应用潜在效应，还能够对空间依赖性进行显式建模。

该模型包含协变量和 误差项 的自回归：

\ [
y = X \ beta + u; u = \ rho Wu + e; e \ sim N（0，\ sigma ^ 2）
\]

\ [
y = X \ beta + \ varepsilon; \ varepsilon \ sim N（0，\ sigma ^ 2（I- \ rho W）^ {-1}（I- \ rho W’）^ {-1}）
\]

该模型包含协变量和响应的自回归：

\ [
y = \ rho W y + X \ beta + e; e \ sim N（0，\ sigma ^ 2）
\]

\ [
y =（I- \ rho W）^ {-1} X \ beta + \ varepsilon; \ \ varepsilon \ sim N（0，\ sigma ^ 2（I- \ rho W）^ {-1}（I- \ rho W’）^ {-1}）
\]

潜在影响`slm`

当初包含一个_试验_所谓的新的潜在影响slm，以合乎以下模型：

\ [
\ mathbf {x} =（I_n- \ rho W）^ {-1}（X \ beta + e）
\]

该模型的元素是：

\（W \）是行标准化的邻接矩阵。
\（\ rho \）是空间自相干参数。
\（X \）是协变量的矩阵，系数为 \（\ beta \）。
\（e \）是具备方差 \（\ sigma ^ 2 \）的高斯 iid 误差。

该 slm 潜成果 的试验，它能够与所述线性预测其余成果组合。

为了定义模型，咱们须要：

X：协变量矩阵
W：行标准化的 邻接矩阵
Q：系数 \（\ beta \）的准确矩阵
范畴 \（\ RHO \），通常由本征值定义

slm潜在作用是通过参数传递 args.sm。在这里，咱们创立了一个具备雷同名称的列表，以将所有必须的值保留在一起：

#Arguments for 'slm'
args.slm = list(
   rho.min = rho.min ,
   rho.max = rho.max,
   W = W,
   X = mmatrix,
   Q.beta = Q.beta
)

此外，还设置了精度参数 \（\ tau \）和空间自相干参数 \（\ rho \）的先验：

#rho 的先验
hyper.slm = list(
   prec = list(prior = "loggamma", param = c(0.01, 0.01)),
      rho = list(initial=0, prior = "logitbeta", param = c(1,1))
)

先前的定义应用具备不同参数的命名列表。参数 prior定义了应用之前 param 及其参数。在此，为精度调配了带有参数 \（0.01 \）和 \（0.01 \）的伽玛先验值，而为空间自相干参数指定了带有参数 \（1 \）和 \（1 \）的 beta 先验值（即 a 区间 \（（（1，1）\））中的平均先验。

## 
## Call:
## Time used:
##     Pre = 0.265, Running = 1.2, Post = 0.058, Total = 1.52 
## Random effects:
##   Name     Model
##     ID SLM model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ID 8.989 4.115      3.709    8.085     19.449 6.641
## Rho for ID       0.822 0.073      0.653    0.832      0.936 0.854
## 
## Expected number of effective parameters(stdev): 62.82(15.46)
## Number of equivalent replicates : 4.47 
## 
## Deviance Information Criterion (DIC) ...............: 902.32
## Deviance Information Criterion (DIC, saturated) ....: 372.95
## Effective number of parameters .....................: 64.13
## 
## Watanabe-Akaike information criterion (WAIC) ...: 905.19
## Effective number of parameters .................: 56.12
## 
## Marginal log-Likelihood:  -477.30 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

系数的预计显示为随机效应的一部分：

round(m.slm$summary.random$ID\[47:48,\], 4)

##    ID   mean     sd 0.025quant 0.5quant 0.975quant   mode kld
## 47 47 0.4634 0.3107    -0.1618   0.4671     1.0648 0.4726   0
## 48 48 0.2606 0.3410    -0.4583   0.2764     0.8894 0.3063   0

空间自相干以外部比例报告（即 0 到 1 之间），并且须要从新缩放：

## Mean            0.644436 
## Stdev           0.145461 
## Quantile  0.025 0.309507 
## Quantile  0.25  0.556851 
## Quantile  0.5   0.663068 
## Quantile  0.75  0.752368 
## Quantile  0.975 0.869702
``````
plot(marg.rho, type = "l", main = "Spatial autocorrelation")

NY8$ICAR <- m.icar$summary.fitted.values\[, "mean"\]
NY8$BYM <- m.bym$summary.fitted.values\[, "mean"\]
NY8$LEROUX <- m.ler$summary.fitted.values\[, "mean"\]
NY8$SLM <- m.slm$summary.fitted.values\[, "mean"\]
spplot(NY8\[syracuse, \], 
  c("FIXED.EFF", "IID.EFF", "ICAR", "BYM", "LEROUX", "SLM"),
  col.regions = rev(magma(16))
)

留神空间模型如何产生绝对危险的更平滑的预计。

为了抉择最佳模型，能够应用下面计算的模型抉择规范：

Bivand, R., E. Pebesma and V. Gómez-Rubio (2013). _Applied spatial data
analysis with R_. Springer-Verlag. New York.

最受欢迎的见解

1.matlab 应用贝叶斯优化的深度学习

2.matlab 贝叶斯隐马尔可夫 hmm 模型实现

3.R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真

4.R 语言中的 block Gibbs 吉布斯采样贝叶斯多元线性回归

5.R 语言中的 Stan 概率编程 MCMC 采样的贝叶斯模型

6.Python 用 PyMC3 实现贝叶斯线性回归模型

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

8.R 语言随机搜寻变量抉择 SSVS 预计贝叶斯向量自回归（BVAR）模型

9.matlab 贝叶斯隐马尔可夫 hmm 模型实现

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

原文链接：http://tecdat.cn/?p=10932

原文出处：拓端数据部落公众号

介绍

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

疾病作图

混合效应模型

泊松回归

具备随机效应的泊松回归

增加点估计以进行映射

晶格数据的空间模型

邻接矩阵

回归模型

线性混合模型

空间随机效应的构造

ICAR 模型

BYM 模型

Leroux 模型

空间计量经济学模型

同步自回归模型（SEM）

空间滞后模型（SLM）

潜在影响`slm`

模型定义

模型拟合

后果汇总

参考文献

Just My Socks（注册教程内含优惠码）

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

原文链接：http://tecdat.cn/?p=10932

原文出处：拓端数据部落公众号

介绍

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

疾病作图

混合效应模型

泊松回归

具备随机效应的泊松回归

增加点估计以进行映射

晶格数据的空间模型

邻接矩阵

回归模型

线性混合模型

空间随机效应的构造

ICAR 模型

BYM 模型

Leroux 模型

空间计量经济学模型

同步自回归模型（SEM）

空间滞后模型（SLM）

潜在影响slm

模型定义

模型拟合

后果汇总

参考文献

Just My Socks（注册教程 内含优惠码）

潜在影响`slm`

Just My Socks（注册教程内含优惠码）