乐趣区

关于数据挖掘:R语言用贝叶斯层次模型进行空间数据分析附代码数据

浏览全文: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 多层(档次)线性模型模型

退出移动版