浏览全文：http://tecdat.cn/?p=10932
最近咱们被客户要求撰写对于贝叶斯层次模型的钻研报告，包含一些图形和统计输入。
在本文中，我将重点介绍应用集成嵌套 拉普拉斯近似办法的贝叶斯推理。能够预计贝叶斯 层次模型的后边缘散布。鉴于模型类型十分宽泛，咱们将重点关注用于剖析晶格数据的空间模型
数据集：纽约州北部的白血病
为了阐明如何与空间模型拟合，将应用纽约白血病数据集。该数据集记录了普查区纽约州北部的许多白血病病例。数据集中的一些变量是：
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净化的站点的地位，能够通过放大看到它。
点击题目查阅往期相干内容
[](http://mp.weixin.qq.com/s?__biz=MzA4MDUzOTIxNA==\&mid=2653830410\&idx=1\&sn=8d1f9df1e68e5e6720451be5a67fe779\&chksm=8478262bb30faf3d26e1559c943a5fdfde75b56405fdbff53a2fd56126b68a061652ec48549b\&scene=21#wechat_redirect)R语言用lme4多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据
左右滑动查看更多
01
02
03
04
混合效应模型
泊松回归
咱们将思考的第一个模型是没有潜在随机效应的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模型
第一个示例将基于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
BYM模型
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 模型
该模型是应用矩阵的“混合”（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
空间计量经济学模型
空间计量经济学是通过 对空间建模略有不同的办法开发的。除了应用潜在效应，还能够对空间 依赖性进行显式建模。
同步自回归模型（SEM）
该模型包含协变量和误差项的自回归：
\ [\
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}）\
]
空间滞后模型（SLM）
该模型包含协变量和响应的自回归：
\ [\
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.
本文摘选 《 R语言应用贝叶斯层次模型进行空间数据分析 》 ，点击“浏览原文”获取全文残缺代码数据资料。
点击题目查阅往期内容
R语言贝叶斯狭义线性混合（多层次/程度/嵌套）模型GLMM、逻辑回归剖析教育留级影响因素数据\
R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例\
用SPSS预计HLM多层（档次）线性模型模型\
R语言用线性混合效应（多程度/档次/嵌套）模型剖析腔调高下与礼貌态度的关系\
R语言LME4混合效应模型钻研老师的受欢迎水平R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model剖析藻类数据实例\
R语言混合线性模型、多层次模型、回归模型剖析学生均匀问题GPA和可视化\
R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例\
R语言用lme4多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据R语言 线性混合效应模型实战案例\
R语言混合效应逻辑回归（mixed effects logistic）模型剖析肺癌数据\
R语言如何用潜类别混合效应模型（LCMM）剖析抑郁症状\
R语言基于copula的贝叶斯分层混合模型的诊断准确性钻研\
R语言建设和可视化混合效应模型mixed effect model\
R语言LME4混合效应模型钻研老师的受欢迎水平\
R语言 线性混合效应模型实战案例\
R语言用Rshiny摸索lme4狭义线性混合模型（GLMM）和线性混合模型（LMM）\
R语言基于copula的贝叶斯分层混合模型的诊断准确性钻研\
R语言如何解决线性混合模型中畸形拟合(Singular fit)的问题\
基于R语言的lmer混合线性回归模型\
R语言用WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型\
R语言分层线性模型案例\
R语言用WinBUGS 软件对学术能力测验（SAT）建设分层模型\
应用SAS，Stata，HLM，R，SPSS和Mplus的分层线性模型HLM\
R语言用WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型\
SPSS中的多层（等级）线性模型Multilevel linear models钻研整容手术数据\
用SPSS预计HLM多层（档次）线性模型模型
