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

在本节中,我将重点介绍使用集成嵌套拉普拉斯近似方法的贝叶斯推理。 

由Kaizong Ye,Coin Ge撰写

可以 估计贝叶斯 层次模型的后边缘分布。 鉴于模型类型非常广泛,我们将重点关注用于分析晶格数据的空间模型。

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

为了说明如何与空间模型拟合,将使用纽约白血病数据集。该数据集记录了普查区纽约州北部的许多白血病病例。数据集中的一些变量是:

  • Cases:1978-1982年期间的白血病病例数。
  • POP8:1980年人口。
  • PCTOWNHOME:拥有房屋的人口比例。
  • PCTAGE65P:65岁以上的人口比例。
  • AVGIDIST:到最近的三氯乙烯(TCE)站点的平均反距离。

鉴于有兴趣研究纽约州北部的白血病风险,因此首先要计算预期的病例数。

这是通过计算总死亡率(总病例数除以总人口数)并将其乘以总人口数得出的:


热门课程

R语言数据分析挖掘必知必会

面对扑面而来的数据浪潮,包含Google、Facebook等国际企业,都已采用R语言进行数据分析

探索课程

一旦获得了预期的病例数,就可以使用标准化死亡率(SMR)来获得原始的风险估计,该标准是将观察到的病例数除以预期的病例数得出的:

疾病作图

在流行病学中,重要的是制作地图以显示相对风险的空间分布。在此示例中,我们将重点放在锡拉库扎市以减少生成地图的计算时间。因此,我们用锡拉丘兹市的区域创建索引:

可以使用函数spplot(在包中sp)简单地创建疾病图:

可以轻松创建交互式地图

请注意,先前的地图还包括11个受TCE污染的站点的位置,可以通过缩小看到它。

混合效应模型

泊松回归

我们将考虑的第一个模型是没有潜在随机效应的Poisson模型,因为这将提供与其他模型进行比较的基准。
 

 模型 :

请注意,它的glm功能类似于该功能。在此,参数 E用于预期的案例数。或  设置了其他参数来计算模型参数的边际
(使用control.predictor)并计算一些模型选择标准 (使用control.compute)。

接下来,可以获得模型的摘要:

具有随机效应的泊松回归

可以通过 在线性预测变量中包括iid高斯随机效应,将潜在随机效应添加到模型中,以解决过度分散问题。

现在,该模式的摘要包括有关随机效果的信息:

添加点估计以进行映射

这两个模型估计 可以被添加到 SpatialPolygonsDataFrame NY8 

 ​

晶格数据的空间模型

格子数据涉及在不同区域(例如,邻里,城市,省,州等)测量的数据。出现空间依赖性是因为相邻区域将显示相似的目标变量值。

  

邻接矩阵

可以使用poly2nbpackage中的函数来计算邻接矩阵 spdep。如果其边界 至少在某一点上接触 ,则此功能会将两个区域视为邻居:

这将返回一个nb具有邻域结构定义的对象:

另外, 当多边形的重心 已知时,可以绘制对象:

 ​

回归模型

通常情况是,除了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-函数定义空间潜在效果。这将需要 一个索引来识别每个区域中的随机效应,模型的类型 和邻接矩阵。为此,将使用稀疏矩阵。

BYM模型

Besag,York和Mollié模型包括两个潜在的随机效应:ICAR 潜在效应和高斯iid潜在效应。线性预测变量eta_i
为:

eta_i = \ alpha + \ beta AVGIDIST_i + u_i + v_i

  • u_i 是iid高斯随机效应
  • v_i 是内在的CAR随机效应

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} )是一个:

该模型与往常一样具有功能inla。注意,C 矩阵使用参数
传递给f函数Cmatrix

空间计量经济学模型

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

同步自回归模型(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。在这里,我们创建了一个具有相同名称的列表,以将 所有必需的值保存在一起:

此外,还设置了精度参数tau 和空间 自相关参数rho 的先验:

先前的定义使用具有不同参数的命名列表。参数 prior定义了使用之前param及其参数。在此,为 精度分配了带有参数0.01和0.01的伽玛先验值,而 为空间自相关参数指定了带有参数1 和1 的beta先验值即均匀先验。

模型拟合

系数的估计显示为随机效应的一部分:

空间自相关以内部比例报告(即 0到1 之间),并且需要重新缩放:

 ​

结果汇总

 注意空间模型如何产生相对风险的更平滑的估计。

为了选择最佳模型, 可以使用上面计算的模型选择标准:

 ​

参考文献

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


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498