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

简介: 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模型

第一个示例将基于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 \)的先验:


#Prior on 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.869702plot(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月前
|
数据采集 机器学习/深度学习 数据可视化
R语言从数据到决策:R语言在商业分析中的实践
【9月更文挑战第1天】R语言在商业分析中的应用广泛而深入,从数据收集、预处理、分析到预测模型构建和决策支持,R语言都提供了强大的工具和功能。通过学习和掌握R语言在商业分析中的实践应用,我们可以更好地利用数据驱动企业决策,提升企业的竞争力和盈利能力。未来,随着大数据和人工智能技术的不断发展,R语言在商业分析领域的应用将更加广泛和深入,为企业带来更多的机遇和挑战。
|
1月前
|
数据采集 数据可视化 数据挖掘
R语言在金融数据分析中的深度应用:探索数据背后的市场智慧
【9月更文挑战第1天】R语言在金融数据分析中展现出了强大的功能和广泛的应用前景。通过丰富的数据处理函数、强大的统计分析功能和优秀的可视化效果,R语言能够帮助金融机构深入挖掘数据价值,洞察市场动态。未来,随着金融数据的不断积累和技术的不断进步,R语言在金融数据分析中的应用将更加广泛和深入。
|
2月前
|
机器学习/深度学习 数据采集 数据可视化
R语言在数据科学中的应用实例:探索与预测分析
【8月更文挑战第31天】通过上述实例,我们展示了R语言在数据科学中的强大应用。从数据准备、探索、预处理到建模与预测,R语言提供了完整的解决方案和丰富的工具集。当然,数据科学远不止于此,随着技术的不断发展和业务需求的不断变化,我们需要不断学习和探索新的方法和工具,以更好地应对挑战,挖掘数据的潜在价值。 未来,随着大数据和人工智能技术的普及,R语言在数据科学领域的应用将更加广泛和深入。我们期待看到更多创新的应用实例,为各行各业的发展注入新的动力。
|
2月前
|
数据采集 存储 数据可视化
R语言时间序列分析:处理与建模时间序列数据的深度探索
【8月更文挑战第31天】R语言作为一款功能强大的数据分析工具,为处理时间序列数据提供了丰富的函数和包。从数据读取、预处理、建模到可视化,R语言都提供了灵活且强大的解决方案。然而,时间序列数据的处理和分析是一个复杂的过程,需要结合具体的应用场景和需求来选择合适的方法和模型。希望本文能为读者在R语言中进行时间序列分析提供一些有益的参考和启示。
|
5月前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
5月前
|
机器学习/深度学习 数据可视化
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为
|
2月前
|
资源调度 数据挖掘
R语言回归分析:线性回归模型的构建与评估
【8月更文挑战第31天】线性回归模型是统计分析中一种重要且实用的工具,能够帮助我们理解和预测自变量与因变量之间的线性关系。在R语言中,我们可以轻松地构建和评估线性回归模型,从而对数据背后的关系进行深入的探索和分析。
|
2月前
|
机器学习/深度学习 数据采集
R语言逻辑回归、GAM、LDA、KNN、PCA主成分分类分析预测房价及交叉验证
上述介绍仅为简要概述,每个模型在实施时都需要仔细调整与优化。为了实现高度精确的预测,模型选择与调参是至关重要的步骤,并且交叉验证是提升模型稳健性的有效途径。在真实世界的房价预测问题中,可能还需要结合地域经济、市场趋势等宏观因素进行综合分析。
52 3
|
5月前
|
数据采集 数据可视化
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
|
5月前
利用R语言进行典型相关分析实战
利用R语言进行典型相关分析实战