贝叶斯地理统计模型R-INLA-1

贝叶斯地理统计模型INLA

本次博客主要讲述如何使用R-INLA软件进行空间分析，通过随机嵌套偏微分方程方法和集成的嵌套Laplace渐进法可为潜在高斯随机场模型中的边际分布提供准确而有效的估计。近年来已经广泛应用于空间流行病学领域。

由于笔者水平有限，关于理论部分，可前往link，针对数学公式及理论部分，这里不赘述，简化数学公式，强调如何应用，及在R语言里面如何实现。

安装INLA包

INLA官网The R-INLA project
如果在R里面下载速度非常慢，可以去 Index source 下载最新版Windows R-INLA 3.6里面，直接下载安装包

# 稳定版
 install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
 
#测试版
install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/testing"), dep=TRUE)

然后在RStudio里面 Tool->install.packages，选择下载的安装包即可。

简述

空间自相关是地理研究中的涉及到的普遍问题。Tobler的第一地理定律：
“所有事物都与其他事物有关，但附近的事物比远处的事物更相关。”
对于空间和时间上的对象都是如此，通常时间与空间是交互作用的。
我们知道，在流行病中，空间分析主要是对疾病数据进行空间上与时间上描述，找出相关性，绘制疾病风险地图，但是实际上空间分析非常复杂，计算量大且不容易直观体现。再叠加时间元素会让让人望而却步。R-INLA出现给解决此类问题提供了便捷的工具，INLA代表集成嵌套拉普拉斯逼近，我们将进一步了解其含义！

INLA使用确定性贝叶斯方法集成嵌套拉普拉斯近似法。
贝叶斯（Bayesian）=使用贝叶斯定理，与概率论相反。
是基于推断给定确定参数的数据集的概率（涉及设置先验！）。如想了解有关更多详细信息，您可以贝叶斯统计入门教程Bayesian Statistics。

1. 案例数据

我们使用gstat包里面自带的降雨数据，里面包含了467个测量站点信息，每个站点都会监测该点的降雨量，然后包含了该地区的海拔高度的图层，我们根据各个站点提取对应位置的海拔高度，然后将数据分成test与train，test100个点作为模型构建，然后剩余的367个点作为validation验证。

估计该区域的降雨量。

library(INLA)
library(gstat)
library(sf)
library(raster)
data(sic97)

## Data preparation
df_rain = st_as_sf(sic_full)
df_train = st_as_sf(sic_obs)

rast_elev <- raster(demstd)
sic_full@data$altitude=extract(rast_elev,sic_full)

# plot
ggplot()+
  geom_sf(data=df_rain,size=0.2,fill=NA)+
  geom_sf(data=df_train,size=0.8,fill=NA,color="red")

image.png

如何判断空间独立性是在进行空间分析前的首要步骤。可以利用变异函数（variogram ）图来评估残差中的空间（或时间）是否相互性。判断空间独立性有一下两点。

• 1.对于随机数据，几乎没有自动相关性，因此distance非常小，我们可以快速到达顶端。
• 2.对于梯度数据（空间上临近位置数据相互关联），我们的distance很长，而且实际上缓慢到达顶端。

这里主要介绍如何判断空间位置的点是否存相互独立，更多关于Variogram，请参考此处Variograms & Kriging

# set the train and test
df_rain = st_as_sf(sic_full) 
xlat=as.data.frame(st_coordinates(df_rain)/1000)

df_rain=data.frame(ID=df_rain$ID,
              rainfall=df_rain$rainfall,
              altitude=df_rain$altitude,
              lat=xlat$X,lon=xlat$Y)

id=sic_obs@data$ID

test=df_rain %>% dplyr::filter(ID %in% id )
train=df_rain %>% dplyr::filter(!ID %in% id )

# test the spatial dependency
sp_df=df_rain
coordinates(sp_df) = ~lat+lon
plot(variogram(rainfall~1,sp_df))

image.png

由上图可以知道，我们的降雨量信息在该区域是存在空间相关性的，那么下一步我们采用Matern函数来定义点与点之间的相关性（空间位置邻近的点相互影响比较远点点大）。这里关于Matern correlation 不做过多介绍，主要的目的是可以计算空间效应。在一般回归方程中，我们可以加入自变量因变量，现在Matern函数定义好了空间效应，那么回归方程：

下面我们将介绍如何计算100个降雨点之间的空间效应，并纳入Regression model

2. INLA模型

INLA模型中，空间效应的计算是重点，这里利用每个测量点的经纬度信息

2.1 Mesh格点

主要经纬度转换时候，需要变成Matrix。为什么要产生Mesh格点，NLA在计算上很有效，因为它使用SPDE（随机偏微分方程）来估计数据的空间自相关。
这涉及使用离散的采样位置的“Mesh网格”，将其插值以估计空间中的连续过程（请参见非常有用的图）。
主要是利用inla.mesh.2d() 产生Constrained Refined Delaunay Triangulation (CRDT)三角形的网格单元，根据坐标位置产生内部三角形研究区域与外部缓冲区。使得所有的采样点都能落在三角形区域内，然后计算每个三角形是否包含采样点位置信息。（详情请见）
inla.mesh.2d() 常用参数：

• loc 用作初始三角剖分节点的坐标矩阵
• loc.domain 定义空间域的多边形的坐标
• max.edge 内部（和外部）区域的最大允许三角形边缘长度。 值应在与坐标大小相关的比例尺上。 值越低，三角形越多。
• offset 扩大点与内部（和外部）区域的边缘之间的距离的量。 将正值视为绝对距离，将负数视为乘数。
• cutoff 点之间允许的最小距离。 这允许将非常靠近的点放置在同一三角形中。 特别需要注意的是，我们不希望三角形的角度非常锐化，因为三角形在投影时效果会较差。

test_loc = cbind(test$lat, test$lon) #Matrix
train_loc = cbind(train$lat, train$lon)

## 1. MESH
Mesh=inla.mesh.2d(loc.domain = test_loc,
                  max.edge = c(20,100),
                  cutoff = 1)
plot(Mesh)
points(test_loc,col = "red", pch = 1)

summary(Mesh)

100个test点全部包含在488个Mesh网格点中，Vertices：488

image.png

2.2 SPDE model

SPDE模型定义在488（m）个尺寸的网格上，而我们的y（n）有100个点。
我们需要一种将m个网格顶点链接到n个响应的方法。
这是通过投影仪矩阵（A）实现的。
该投影仪矩阵是使用inla.spde.make.A（）函数构建的。链接好好，我们计算Matern相关函数，在网格上构建SPDE模型。 alpha = 2是默认值。

# Making the A matrix
A.test=inla.spde.make.A(Mesh,loc=test_loc)
dim(A.test)

## SPDE
spde=inla.spde2.matern(mesh=Mesh,alpha = 2)
spde$n.spde

## index of spatial field
s.index=inla.spde.make.index(name="w",
                             n.spde = spde$n.spde)

str(s.index)

SPDE模型非常复杂。因此，为了帮助跟踪哪些元素与哪些效果相关，我们可以创建一个索引Index。注意这里的name是w，可以写成spatial feild，意思是每个点对应的空间效应。在这种情况下，我们的空间数据全部在一组中。

2.4 Stack data

在2.1中，我们告知R-INLA我们在网格的哪些顶点具有采样位置，这给了我们投影仪矩阵A.test。 在第2.2节中，我们定义了SPDE模型。 我们需要告知R-INLA，在哪些采样位置我们有y（response）的数据以及在哪里有x（协变量）数据。 由于协变量可能在与响应变量存在于不同位置，因此这一步我们需要整和协变量。使用函数inla.stack

## stack
Xm <- model.matrix(~ -1 + altitude, data = test)

X=data.frame(altitude=Xm[,1])

N=nrow(test)
N
## 5.1 test stack
stack.test=inla.stack(tag="test",
  data=list(y=test$rainfall),
  A = list(1,1,A.test),
  effects= list(
    Intercept = rep(1, N),
    X = X,
    w = s.index))
dim(inla.stack.A(stack.test))

stack.test中包括：

• data 响应变量（y）
  
• A 乘数向量。 通常是一系列1（用于截距，随机效果和固定效果），及指定的空间A矩阵。
  
• effect 效果。 需要分别指定截距，随机效果，模型矩阵和SPDE。 A与effect是相互联系的。 要添加效果，必须在乘数向量上再加上1（在正确的位置）。
  

2.5 Fit model

## 7. fit model
formula = y ~ -1+Intercept+altitude+f(w, model = spde)

fit=inla(formula = formula,
         data = inla.stack.data(stack.test,spde=spde),
         family = "gaussian",
         control.compute = list(dic = TRUE,waic = TRUE),
         control.predictor = list(A = inla.stack.A(stack.test),compute=TRUE)
)

summary(fit)

这里100个test数据显示海拔的系数为0.005（-0.03-0.04）没有统计学意义。
说明海拔对该地区的降雨没有多大影响。

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0.000 31.623 -62.086 -0.001 62.034 0.000 0
altitude 0.005 0.018 -0.030 0.005 0.040 0.005 0

模型参数估计

边际效应；详情见 INLA 参数介绍

模型估计

ggField(fit, Mesh, Groups = 1) +
  scale_fill_brewer(palette = "Blues") 

plot(SpFi.w$marginals.range.nominal[[1]], type = "l",
     xlab = "range nominal", ylab = "Density") 

image.png
image.png

模型验证

未完待续。。。