貝葉斯地理統(tǒng)計(jì)模型INLA
本次博客主要講述如何使用R-INLA軟件進(jìn)行空間分析，通過(guò)隨機(jī)嵌套偏微分方程方法和集成的嵌套Laplace漸進(jìn)法可為潛在高斯隨機(jī)場(chǎng)模型中的邊際分布提供準(zhǔn)確而有效的估計(jì)虫几。近年來(lái)已經(jīng)廣泛應(yīng)用于空間流行病學(xué)領(lǐng)域锤灿。
由于筆者水平有限，關(guān)于理論部分辆脸，可前往link但校，針對(duì)數(shù)學(xué)公式及理論部分，這里不贅述啡氢，簡(jiǎn)化數(shù)學(xué)公式状囱，強(qiáng)調(diào)如何應(yīng)用，及在R語(yǔ)言里面如何實(shí)現(xiàn)空执。
安裝INLA包
INLA官網(wǎng)The R-INLA project
如果在R里面下載速度非常慢浪箭，可以去 Index source 下載最新版Windows R-INLA 3.6里面穗椅，直接下載安裝包
# 穩(wěn)定版
install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
#測(cè)試版
install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/testing"), dep=TRUE)
然后在RStudio里面 Tool->install.packages辨绊，選擇下載的安裝包即可。
簡(jiǎn)述
空間自相關(guān)是地理研究中的涉及到的普遍問(wèn)題匹表。Tobler的第一地理定律：
“所有事物都與其他事物有關(guān)门坷，但附近的事物比遠(yuǎn)處的事物更相關(guān)∨鄱疲”
對(duì)于空間和時(shí)間上的對(duì)象都是如此默蚌，通常時(shí)間與空間是交互作用的。
我們知道苇羡，在流行病中绸吸，空間分析主要是對(duì)疾病數(shù)據(jù)進(jìn)行空間上與時(shí)間上描述，找出相關(guān)性设江，繪制疾病風(fēng)險(xiǎn)地圖锦茁，但是實(shí)際上空間分析非常復(fù)雜，計(jì)算量大且不容易直觀體現(xiàn)叉存。再疊加時(shí)間元素會(huì)讓讓人望而卻步码俩。R-INLA出現(xiàn)給解決此類問(wèn)題提供了便捷的工具，INLA代表集成嵌套拉普拉斯逼近歼捏，我們將進(jìn)一步了解其含義稿存！
INLA使用確定性貝葉斯方法集成嵌套拉普拉斯近似法笨篷。
貝葉斯（Bayesian）=使用貝葉斯定理，與概率論相反瓣履。
是基于推斷給定確定參數(shù)的數(shù)據(jù)集的概率（涉及設(shè)置先驗(yàn)Ｂ食帷）。如想了解有關(guān)更多詳細(xì)信息袖迎，您可以貝葉斯統(tǒng)計(jì)入門教程Bayesian Statistics安聘。
1. 案例數(shù)據(jù)
我們使用gstat包里面自帶的降雨數(shù)據(jù)，里面包含了467個(gè)測(cè)量站點(diǎn)信息瓢棒，每個(gè)站點(diǎn)都會(huì)監(jiān)測(cè)該點(diǎn)的降雨量浴韭，然后包含了該地區(qū)的海拔高度的圖層，我們根據(jù)各個(gè)站點(diǎn)提取對(duì)應(yīng)位置的海拔高度脯宿，然后將數(shù)據(jù)分成test與train念颈，test100個(gè)點(diǎn)作為模型構(gòu)建，然后剩余的367個(gè)點(diǎn)作為validation驗(yàn)證连霉。
估計(jì)該區(qū)域的降雨量榴芳。
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")
如何判斷空間獨(dú)立性是在進(jìn)行空間分析前的首要步驟《搴常可以利用變異函數(shù)（variogram ）圖來(lái)評(píng)估殘差中的空間（或時(shí)間）是否相互性窟感。判斷空間獨(dú)立性有一下兩點(diǎn)。
- 1.對(duì)于隨機(jī)數(shù)據(jù)歉井，幾乎沒(méi)有自動(dòng)相關(guān)性柿祈，因此distance非常小，我們可以快速到達(dá)頂端哩至。
- 2.對(duì)于梯度數(shù)據(jù)（空間上臨近位置數(shù)據(jù)相互關(guān)聯(lián)）躏嚎，我們的distance很長(zhǎng)，而且實(shí)際上緩慢到達(dá)頂端菩貌。
這里主要介紹如何判斷空間位置的點(diǎn)是否存相互獨(dú)立卢佣，更多關(guān)于Variogram，請(qǐng)參考此處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))
由上圖可以知道箭阶，我們的降雨量信息在該區(qū)域是存在空間相關(guān)性的虚茶，那么下一步我們采用Matern函數(shù)來(lái)定義點(diǎn)與點(diǎn)之間的相關(guān)性（空間位置鄰近的點(diǎn)相互影響比較遠(yuǎn)點(diǎn)點(diǎn)大）。這里關(guān)于Matern correlation 不做過(guò)多介紹仇参，主要的目的是可以計(jì)算空間效應(yīng)嘹叫。在一般回歸方程中，我們可以加入自變量因變量冈敛，現(xiàn)在Matern函數(shù)定義好了空間效應(yīng)待笑，那么回歸方程：
下面我們將介紹如何計(jì)算100個(gè)降雨點(diǎn)之間的空間效應(yīng)，并納入Regression model
2. INLA模型
INLA模型中抓谴，空間效應(yīng)的計(jì)算是重點(diǎn)暮蹂，這里利用每個(gè)測(cè)量點(diǎn)的經(jīng)緯度信息
2.1 Mesh格點(diǎn)
主要經(jīng)緯度轉(zhuǎn)換時(shí)候寞缝，需要變成Matrix。為什么要產(chǎn)生Mesh格點(diǎn)仰泻，NLA在計(jì)算上很有效荆陆，因?yàn)樗褂肧PDE（隨機(jī)偏微分方程）來(lái)估計(jì)數(shù)據(jù)的空間自相關(guān)。
這涉及使用離散的采樣位置的“Mesh網(wǎng)格”集侯，將其插值以估計(jì)空間中的連續(xù)過(guò)程（請(qǐng)參見(jiàn)非常有用的圖）被啼。
主要是利用inla.mesh.2d() 產(chǎn)生Constrained Refined Delaunay Triangulation (CRDT)三角形的網(wǎng)格單元，根據(jù)坐標(biāo)位置產(chǎn)生內(nèi)部三角形研究區(qū)域與外部緩沖區(qū)棠枉。使得所有的采樣點(diǎn)都能落在三角形區(qū)域內(nèi)浓体，然后計(jì)算每個(gè)三角形是否包含采樣點(diǎn)位置信息。（詳情請(qǐng)見(jiàn)）
inla.mesh.2d() 常用參數(shù)：
- loc 用作初始三角剖分節(jié)點(diǎn)的坐標(biāo)矩陣
- loc.domain 定義空間域的多邊形的坐標(biāo)
- max.edge 內(nèi)部（和外部）區(qū)域的最大允許三角形邊緣長(zhǎng)度辈讶。 值應(yīng)在與坐標(biāo)大小相關(guān)的比例尺上命浴。 值越低，三角形越多贱除。
- offset 擴(kuò)大點(diǎn)與內(nèi)部（和外部）區(qū)域的邊緣之間的距離的量生闲。 將正值視為絕對(duì)距離，將負(fù)數(shù)視為乘數(shù)月幌。
- cutoff 點(diǎn)之間允許的最小距離碍讯。 這允許將非常靠近的點(diǎn)放置在同一三角形中扯躺。 特別需要注意的是捉兴，我們不希望三角形的角度非常銳化，因?yàn)槿切卧谕队皶r(shí)效果會(huì)較差缅帘。
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個(gè)test點(diǎn)全部包含在488個(gè)Mesh網(wǎng)格點(diǎn)中轴术，Vertices：488
2.2 SPDE model
SPDE模型定義在488（m）個(gè)尺寸的網(wǎng)格上难衰，而我們的y（n）有100個(gè)點(diǎn)钦无。
我們需要一種將m個(gè)網(wǎng)格頂點(diǎn)鏈接到n個(gè)響應(yīng)的方法。
這是通過(guò)投影儀矩陣（A）實(shí)現(xiàn)的盖袭。
該投影儀矩陣是使用inla.spde.make.A（）函數(shù)構(gòu)建的失暂。鏈接好好，我們計(jì)算Matern相關(guān)函數(shù)鳄虱，在網(wǎng)格上構(gòu)建SPDE模型弟塞。 alpha = 2是默認(rèn)值。
# 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模型非常復(fù)雜拙已。因此决记，為了幫助跟蹤哪些元素與哪些效果相關(guān)，我們可以創(chuàng)建一個(gè)索引Index倍踪。注意這里的name是w系宫，可以寫成spatial feild索昂，意思是每個(gè)點(diǎn)對(duì)應(yīng)的空間效應(yīng)。在這種情況下扩借，我們的空間數(shù)據(jù)全部在一組中椒惨。
2.4 Stack data
在2.1中，我們告知R-INLA我們?cè)诰W(wǎng)格的哪些頂點(diǎn)具有采樣位置潮罪，這給了我們投影儀矩陣A.test康谆。 在第2.2節(jié)中，我們定義了SPDE模型嫉到。 我們需要告知R-INLA沃暗，在哪些采樣位置我們有y（response）的數(shù)據(jù)以及在哪里有x（協(xié)變量）數(shù)據(jù)。 由于協(xié)變量可能在與響應(yīng)變量存在于不同位置何恶，因此這一步我們需要整和協(xié)變量描睦。使用函數(shù)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īng)變量（y）
- A 乘數(shù)向量。 通常是一系列1（用于截距导而，隨機(jī)效果和固定效果）忱叭，及指定的空間A矩陣。
- effect 效果今艺。 需要分別指定截距韵丑，隨機(jī)效果，模型矩陣和SPDE虚缎。 A與effect是相互聯(lián)系的撵彻。 要添加效果，必須在乘數(shù)向量上再加上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個(gè)test數(shù)據(jù)顯示海拔的系數(shù)為0.005（-0.03-0.04）沒(méi)有統(tǒng)計(jì)學(xué)意義陌僵。
說(shuō)明海拔對(duì)該地區(qū)的降雨沒(méi)有多大影響。
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
模型參數(shù)估計(jì)
邊際效應(yīng)创坞；詳情見(jiàn) INLA 參數(shù)介紹
模型估計(jì)
ggField(fit, Mesh, Groups = 1) +
scale_fill_brewer(palette = "Blues")
plot(SpFi.w$marginals.range.nominal[[1]], type = "l",
xlab = "range nominal", ylab = "Density")
模型驗(yàn)證
未完待續(xù)碗短。。题涨。
