全文鏈接：https://tecdat.cn/?p=40720

本論文旨在為對(duì)空間建模感興趣的研究人員客戶提供使用R-INLA進(jìn)行空間數(shù)據(jù)建模的基礎(chǔ)教程。通過(guò)對(duì)區(qū)域數(shù)據(jù)和地統(tǒng)計(jì)（標(biāo)記點(diǎn)）數(shù)據(jù)的分析，介紹了如何擬合簡(jiǎn)單模型、構(gòu)建和運(yùn)行更復(fù)雜的空間模型，以及繪制空間預(yù)測(cè)和高斯隨機(jī)場(chǎng)等內(nèi)容，幫助讀者掌握空間建模的基本方法和技能，為進(jìn)一步的空間數(shù)據(jù)分析和研究奠定基礎(chǔ)（點(diǎn)擊文末“閱讀原文”獲取完整代碼、數(shù)據(jù)、文檔）。

關(guān)鍵詞

空間建模；R-INLA；區(qū)域數(shù)據(jù)；地統(tǒng)計(jì)數(shù)據(jù)；高斯隨機(jī)場(chǎng)

一、引言

在當(dāng)今的數(shù)據(jù)分析領(lǐng)域，空間數(shù)據(jù)的分析越來(lái)越受到關(guān)注。許多研究領(lǐng)域，如流行病學(xué)、生態(tài)學(xué)和社會(huì)科學(xué)等，都涉及到空間數(shù)據(jù)的處理和分析?？臻g建模的目的是將數(shù)據(jù)的空間結(jié)構(gòu)納入模型中，以更準(zhǔn)確地描述和解釋數(shù)據(jù)。R-INLA（ Integrated Nested Laplace Approximations）是一個(gè)在R語(yǔ)言中用于空間建模的強(qiáng)大工具，它提供了一種高效的方法來(lái)估計(jì)空間自相關(guān)結(jié)構(gòu)和擬合空間模型。
本教程的目標(biāo)是向讀者展示使用R-INLA進(jìn)行空間數(shù)據(jù)建模的基本步驟和方法，旨在成為任何對(duì)空間建模感興趣的人的起點(diǎn)。讀者將能夠掌握以下技能：

1. 學(xué)會(huì)在區(qū)域數(shù)據(jù)上擬合簡(jiǎn)單模型。

2. 了解地統(tǒng)計(jì)（標(biāo)記點(diǎn)）數(shù)據(jù)建模的基礎(chǔ)知識(shí)。

3. 構(gòu)建和運(yùn)行更復(fù)雜的空間模型。

4. 繪制空間預(yù)測(cè)和高斯隨機(jī)場(chǎng)。

二、所需知識(shí)和軟件包

本文假設(shè)讀者具備廣義線性模型（GLMs）和廣義線性混合模型（GLMMs）、貝葉斯統(tǒng)計(jì)以及空間數(shù)據(jù)處理（特別是柵格數(shù)據(jù)處理）的基礎(chǔ)知識(shí)。
在開始本教程之前，需要下載并安裝一些相關(guān)的軟件包，包括RColorBrewer、spdep、sp、rgdal、raster和INLA。其中，INLA的安裝可能需要一些時(shí)間，建議在開始教程之前完成安裝。

#?添加dep?=?T表示您還將安裝軟件包依賴項(xiàng)
install.packages("RColorBrewer",?dep?=?T)
install.packages("spdep",?dep?=?T)
install.packages("sp",?dep?=?T)
install.packages("rgdal",?dep?=?T)
install.packages("raster",?dep?=?T)
#?下載軟件包的最新穩(wěn)定版本
install.packages("INLA",repos?=?c(getOption("repos"),INLA?=?"https://inla.r-inla-download.org/R/stable"),dep?=?T)

三、數(shù)據(jù)集

本文將使用兩個(gè)數(shù)據(jù)集，這些數(shù)據(jù)集來(lái)自于進(jìn)行的實(shí)地調(diào)查。研究的目的是收集市周圍公共綠地中的狐貍糞便（即糞便標(biāo)記），并分析其中的胃腸道寄生蟲。

研究的問(wèn)題是：綠地面積的數(shù)量是否與以下因素顯著相關(guān)：

1. 發(fā)現(xiàn)的狐貍糞便數(shù)量？

2. 每個(gè)糞便中發(fā)現(xiàn)的寄生蟲物種數(shù)量（物種豐富度）？
   使用的數(shù)據(jù)包括發(fā)現(xiàn)的糞便數(shù)量的區(qū)域數(shù)據(jù)（如圖中的六邊形網(wǎng)格）和每個(gè)樣本中發(fā)現(xiàn)的寄生蟲豐富度的點(diǎn)數(shù)據(jù)。

四、區(qū)域數(shù)據(jù)分析

區(qū)域數(shù)據(jù)通常在流行病學(xué)、生態(tài)學(xué)或社會(huì)科學(xué)研究中出現(xiàn)。簡(jiǎn)單來(lái)說(shuō)，區(qū)域數(shù)據(jù)報(bào)告每個(gè)區(qū)域的一個(gè)值（通常是疾病的病例數(shù)），這些區(qū)域可以是行政區(qū)，如郵政編碼區(qū)、議會(huì)區(qū)、地區(qū)等。區(qū)域數(shù)據(jù)的主要特點(diǎn)是每個(gè)區(qū)域都有明確的鄰居，這使得計(jì)算自相關(guān)結(jié)構(gòu)更加容易。區(qū)域數(shù)據(jù)的一個(gè)特殊子集是格網(wǎng)數(shù)據(jù)，它報(bào)告來(lái)自規(guī)則網(wǎng)格單元的區(qū)域數(shù)據(jù)（就像這里的數(shù)據(jù)一樣）。這種類型的區(qū)域數(shù)據(jù)通常更受歡迎，因?yàn)榭臻g被分割成更具可比性的區(qū)域，并且空間離散化更加均勻。然而，這種類型的區(qū)域數(shù)據(jù)很少見(jiàn)，因?yàn)楦窬W(wǎng)數(shù)據(jù)通常是專門從點(diǎn)數(shù)據(jù)構(gòu)建的（在這種情況下，最好直接使用點(diǎn)數(shù)據(jù)），而實(shí)際的區(qū)域數(shù)據(jù)通常來(lái)自于在行政區(qū)級(jí)別進(jìn)行的調(diào)查，這些區(qū)域的形狀自然不規(guī)則。
在INLA中對(duì)區(qū)域數(shù)據(jù)進(jìn)行建模相對(duì)簡(jiǎn)單（至少與點(diǎn)數(shù)據(jù)集相比）。這是因?yàn)閰^(qū)域已經(jīng)有明確的鄰居（您可以通過(guò)查看圖形直接判斷哪些單元相鄰）。這意味著我們需要做的就是將其轉(zhuǎn)換為鄰接矩陣，以INLA能夠理解的方式指定數(shù)據(jù)集的鄰居系統(tǒng)，然后我們就可以直接擬合模型（這與點(diǎn)數(shù)據(jù)集的情況完全不同）。

4.1 數(shù)據(jù)加載和可視化

為了進(jìn)行空間分析，我們將使用一個(gè)經(jīng)過(guò)修改的數(shù)據(jù)集，該數(shù)據(jù)集涉及在市進(jìn)行的為期6個(gè)月的每個(gè)公共綠地調(diào)查中發(fā)現(xiàn)的狐貍糞便數(shù)量。我們構(gòu)建了一個(gè)覆蓋研究區(qū)域的格網(wǎng)，并記錄了每個(gè)區(qū)域中發(fā)現(xiàn)的糞便數(shù)量，以及使用綠地比率。

#?加載格網(wǎng)shapefile文件和狐貍糞便數(shù)據(jù)
require(sp)?#?用于處理空間數(shù)據(jù)的軟件包
require(rgdal)?#?用于處理空間數(shù)據(jù)的軟件包
#?Fox_Lattice是一個(gè)空間對(duì)象，包含根據(jù)數(shù)據(jù)構(gòu)建的多邊形（通常您會(huì)使用行政區(qū)）
Fox_Lattice?<-?readOGR("F.shp")
#Warning?message:
#?忽略此警告消息，這是因?yàn)闆](méi)有為每個(gè)單元分配z值（我們將響應(yīng)值作為數(shù)據(jù)框附加）
require(RColorBrewer)
#?創(chuàng)建一個(gè)顏色調(diào)色板用于圖形
my.palette?<-?brewer.pal(n?=?9,?name?=?"YlOrRd")
#?可視化空間中糞便的數(shù)量
spplot(obj?=?_La,?cuts?=?8)

圖1：空間中狐貍糞便的數(shù)量

4.2 鄰接矩陣計(jì)算

如前所述，INLA需要知道哪些區(qū)域是相鄰的，以便計(jì)算空間自相關(guān)結(jié)構(gòu)。我們通過(guò)計(jì)算鄰接矩陣來(lái)實(shí)現(xiàn)這一點(diǎn)。

這個(gè)矩陣顯示了每個(gè)單元的鄰居關(guān)系。在兩個(gè)軸上都有單元的數(shù)字ID（ZONE_CODE），您可以找到它們與哪些單元相鄰（加上對(duì)角線表示單元與自身相鄰）。例如，您可以用眼睛追蹤單元編號(hào)50，并查看它的鄰居（單元49和51）。每行最多有6個(gè)鄰居（六邊形有6條邊），對(duì)應(yīng)于格網(wǎng)單元的鄰居數(shù)量。請(qǐng)注意，在這種情況下，單元已經(jīng)按字母順序排序，因此它們只與名稱相似的單元相鄰，所以在對(duì)角線周圍有一組相鄰的單元。當(dāng)使用行政區(qū)時(shí)，這個(gè)矩陣可能會(huì)更混亂。

圖2：鄰接矩陣

4.3 模型公式指定

我們還需要指定模型公式。這個(gè)模型將測(cè)試綠地比率（GS_ratio）對(duì)每個(gè)區(qū)域中發(fā)現(xiàn)的狐貍糞便數(shù)量是否有線性影響。我們首先指定模型公式，這實(shí)際上并不會(huì)運(yùn)行我們的模型，我們將在下一步運(yùn)行模型。

formula?<-?Scat\_No?~?1?+?GS\_Ratio?+?#?固定效應(yīng)f(ZONE\_CO,?#?空間效應(yīng)：ZONE\_CODE是格網(wǎng)中每個(gè)區(qū)域的數(shù)字標(biāo)識(shí)符（不適用于因子）graph?=?La)?#?這指定了格網(wǎng)區(qū)域的鄰居關(guān)系

注意：空間效應(yīng)使用BYM（Besag, York和Mollie的模型）進(jìn)行建模，這是通常用于擬合區(qū)域數(shù)據(jù)的模型類型。CAR（條件自回歸）和besag模型是其他選項(xiàng)，但在這里我們將重點(diǎn)關(guān)注BYM，因?yàn)檫@是在處理區(qū)域數(shù)據(jù)時(shí)對(duì)空間效應(yīng)進(jìn)行建模的合適方法。

4.4 模型運(yùn)行和結(jié)果分析

現(xiàn)在我們準(zhǔn)備運(yùn)行我們的模型！

#?最后，我們可以使用inla()函數(shù)運(yùn)行模型
Mod_Lattice?<-?inla(formula,?
family?=?"poisson",?#?因?yàn)槲覀冋谔幚碛?jì)數(shù)數(shù)據(jù)data?=?Lattice_Data,control.compute?=?list(cpo?=?T,?dic?=?T,?waic?=?T))?
#?CPO,?DIC和WAIC度量值都可以通過(guò)在control.compute選項(xiàng)中指定來(lái)計(jì)算
#?如果您想進(jìn)行模型選擇，這些值可以用于此目的
#?查看模型摘要
summary(Mod_Lattice)

我們現(xiàn)在已經(jīng)運(yùn)行了我們的第一個(gè)INLA模型！
在輸出中，您可以找到關(guān)于模型的一些一般信息：運(yùn)行時(shí)間、固定效應(yīng)的摘要、模型選擇標(biāo)準(zhǔn)（如果您在模型中指定了它們），以及任何隨機(jī)效應(yīng)的精度（在這種情況下只是我們的空間組件ZONE_CODE）。重要的是要記住，INLA使用精度（tau = 1/方差），所以精度值越高，方差值越低。
我們可以看到，GS_Ratio對(duì)發(fā)現(xiàn)的糞便數(shù)量有正影響（0.025分位數(shù)和0.075分位數(shù)不跨越零，所以這是一個(gè)“顯著”的正影響），并且ZONE_CODE的iid（隨機(jī)因子效應(yīng)）的精度比空間效應(yīng)低得多，這意味著在這種情況下，使用ZONE_CODE作為標(biāo)準(zhǔn)的因子隨機(jī)效應(yīng)可能就足夠了。

4.5 設(shè)置先驗(yàn)

我們還可以通過(guò)在公式中指定超參數(shù)（先驗(yàn)分布的參數(shù)）的先驗(yàn)來(lái)設(shè)置先驗(yàn)。INLA使用精度（tau = 1/方差），所以默認(rèn)情況下，非常低的精度對(duì)應(yīng)于非常高的方差。請(qǐng)記住，需要為模型的線性預(yù)測(cè)器指定先驗(yàn)（因此需要根據(jù)數(shù)據(jù)分布進(jìn)行轉(zhuǎn)換），在這種情況下，它們遵循對(duì)數(shù)伽馬分布（因?yàn)檫@是一個(gè)泊松模型）。

隨機(jī)（空間）效應(yīng)的后驗(yàn)均值也可以計(jì)算并繪制在格網(wǎng)上。為此，我們需要提取格網(wǎng)中每個(gè)單元的空間效應(yīng)的后驗(yàn)均值（使用emarginal()函數(shù)），然后將其添加到原始shapefile中，以便我們可以映射它。
這表示在考慮了模型中包含的協(xié)變量后，響應(yīng)變量在空間中的分布?？梢詫⑵湟暈楦鶕?jù)模型的響應(yīng)變量在空間中的“真實(shí)分布”（顯然，這僅與我們擁有的模型一樣好，如果估計(jì)不佳、我們有缺失數(shù)據(jù)或我們未能在模型中包含重要的協(xié)變量，它將受到影響）。

#?計(jì)算區(qū)域數(shù)量
Nares?<-?length(Lat_Data\[,1\])
#?選擇每個(gè)區(qū)域的空間隨機(jī)效應(yīng)的后驗(yàn)邊緣分布
#?這些對(duì)應(yīng)于空間隨機(jī)效應(yīng)（zeta）的邊緣分布的前347個(gè)（單元數(shù)量）項(xiàng)
zone.ix?<-?Mod\_Lat$ZONE\_CODE\[1:Nareas\]
#?對(duì)每個(gè)區(qū)域的邊緣進(jìn)行指數(shù)運(yùn)算，以將其返回到原始值（請(qǐng)記住，這是一個(gè)泊松模型，因此模型的所有組件都進(jìn)行了對(duì)數(shù)變換）

現(xiàn)在我們準(zhǔn)備創(chuàng)建一個(gè)顏色調(diào)色板并制作我們的地圖！

spplot(obj?=?Fost,?zcol

圖3：根據(jù)我們的模型繪制的空間后驗(yàn)均值，顯示狐貍糞便的數(shù)量
同樣，我們可以繪制與后驗(yàn)均值相關(guān)的不確定性。與任何建模一樣，不僅要考慮均值，還要考慮我們對(duì)該均值的信心程度。

圖4：根據(jù)我們的模型繪制的空間后驗(yàn)均值的不確定性
請(qǐng)注意，后驗(yàn)均值在不確定性最高的地方最高。我們有一些區(qū)域，響應(yīng)變量達(dá)到非常高的數(shù)值，這是因?yàn)檫@些區(qū)域缺少GS數(shù)據(jù)（GS = 0），所以模型對(duì)此進(jìn)行了補(bǔ)償；然而，這些區(qū)域也是我們不確定性最高的地方，因?yàn)槟Ｐ蜔o(wú)法產(chǎn)生準(zhǔn)確的估計(jì)。

五、地統(tǒng)計(jì)（標(biāo)記點(diǎn)）數(shù)據(jù)建模

對(duì)于本分析，我們將使用地統(tǒng)計(jì)數(shù)據(jù)，也稱為標(biāo)記點(diǎn)數(shù)據(jù)。這是最常見(jiàn)的空間數(shù)據(jù)類型之一。它包括點(diǎn)（帶有相關(guān)坐標(biāo)），每個(gè)點(diǎn)都有一個(gè)附加的值，通常是我們感興趣的響應(yīng)變量的測(cè)量值。其思想是這些點(diǎn)是在空間中無(wú)處不在的平滑空間過(guò)程的實(shí)現(xiàn)，而這些點(diǎn)只是這個(gè)過(guò)程的樣本（我們永遠(yuǎn)無(wú)法對(duì)整個(gè)過(guò)程進(jìn)行采樣，因?yàn)樵谶B續(xù)空間中有無(wú)限個(gè)點(diǎn)）。
一個(gè)經(jīng)典的例子是土壤pH值：這是土壤的一種屬性，它無(wú)處不在，但我們只會(huì)在某些位置測(cè)量它。通過(guò)將我們收集的值與其他測(cè)量值聯(lián)系起來(lái)，我們可以發(fā)現(xiàn)土壤pH值取決于降水水平或植被類型，并且（有足夠的信息）我們可以重建潛在的空間過(guò)程。
我們通常對(duì)理解潛在過(guò)程（哪些變量影響它？它在空間和時(shí)間上如何變化？）以及重建它（通過(guò)生成模型預(yù)測(cè)）感興趣。
在這個(gè)例子中，我們將使用與生成空間數(shù)據(jù)相同的點(diǎn)（狐貍糞便），但我們將關(guān)注每個(gè)糞便中發(fā)現(xiàn)的寄生蟲物種數(shù)量（Spp_Rich）。數(shù)據(jù)集包括每個(gè)點(diǎn)的位置（每個(gè)點(diǎn)都是在調(diào)查中發(fā)現(xiàn)的糞便），但我們?cè)诖烁信d趣的是對(duì)每個(gè)糞便中發(fā)現(xiàn)的寄生蟲物種數(shù)量進(jìn)行建模。這意味著數(shù)據(jù)集中的每個(gè)點(diǎn)都有一個(gè)附加的值（一個(gè)標(biāo)記，因此稱為標(biāo)記點(diǎn)過(guò)程），這就是我們感興趣的建模內(nèi)容。在這種情況下，點(diǎn)沒(méi)有明確的鄰居，因此我們需要構(gòu)建空間的人工離散化，并告訴INLA離散化的鄰居結(jié)構(gòu)。
數(shù)據(jù)集還包含與每個(gè)樣本相關(guān)的許多其他變量：

• JanDate（收集樣本的日期）

• Site（從哪個(gè)公園收集的）

• Greenspace variability（GS_Var），這是一個(gè)分類變量，用于測(cè)量不同綠地類型的數(shù)量（低、中、高）
  在這種情況下，我們將把胃腸道寄生蟲的物種豐富度建模為綠地比率的函數(shù)，同時(shí)考慮空間效應(yīng)和上述其他協(xié)變量。

當(dāng)將點(diǎn)數(shù)據(jù)集轉(zhuǎn)換為空間對(duì)象時(shí)，我們需要指定一個(gè)坐標(biāo)參考系統(tǒng)（CRS）。此數(shù)據(jù)集的坐標(biāo)以東坐標(biāo)/北坐標(biāo)表示，并使用國(guó)家網(wǎng)格（BNG）進(jìn)行投影。如果您使用多個(gè)可能不在同一坐標(biāo)系統(tǒng)中的shapefile文件，這一點(diǎn)很重要，它們將必須相應(yīng)地進(jìn)行投影。
注意：CRS的選擇應(yīng)基于研究區(qū)域的范圍。

• 小區(qū)域：對(duì)于小區(qū)域（如此處的區(qū)域），東坐標(biāo) - 北坐標(biāo)系統(tǒng)是最好的。它們有效地在平面上表示坐標(biāo)（不考慮地球曲率以及因此投影形狀的修改）。

• 中等規(guī)模研究：對(duì)于中等規(guī)模的研究（國(guó)家級(jí)別/多個(gè)國(guó)家級(jí)別），我們應(yīng)該使用緯度 - 經(jīng)度，因?yàn)檫@將考慮更現(xiàn)實(shí)的地圖形狀。

• 大陸和全球規(guī)模研究：最后，對(duì)于在大陸和全球規(guī)模進(jìn)行的研究，我們應(yīng)該使用弧度，并在考慮地球曲率的情況下擬合網(wǎng)格。

嘗試在同一地圖上繪制我們的點(diǎn)和shapefile文件將不起作用，因?yàn)樗鼈兊淖鴺?biāo)表示在不同的系統(tǒng)中，無(wú)法直接繪制。

圖5：混合不同的坐標(biāo)系統(tǒng)會(huì)導(dǎo)致錯(cuò)誤的圖形！
然而，如果我們使用spTransform()函數(shù)更改Scot_Shape的CRS，我們可以正確地將狐貍糞便點(diǎn)和shapefile文件一起映射。

圖6：
現(xiàn)在數(shù)據(jù)已正確加載，我們可以開始將地統(tǒng)計(jì)INLA模型所需的所有組件組合在一起。我們將首先擬合一個(gè)簡(jiǎn)單的基礎(chǔ)模型，其中僅包含一個(gè)截距和空間效應(yīng)，然后在此基礎(chǔ)上增加模型的復(fù)雜性。

5.1 模型的基本組件

模型的絕對(duì)必要組件包括：

• 網(wǎng)格（Mesh）：與區(qū)域數(shù)據(jù)不同，點(diǎn)數(shù)據(jù)沒(méi)有明確的鄰居，因此我們必須計(jì)算空間中每個(gè)可能點(diǎn)之間的自相關(guān)結(jié)構(gòu)，這顯然是不可能的。因此，第一步是離散化空間以創(chuàng)建一個(gè)網(wǎng)格，該網(wǎng)格將創(chuàng)建人工（但有用）的鄰居集，以便我們可以計(jì)算點(diǎn)之間的自相關(guān)。INLA使用三角形網(wǎng)格，因?yàn)樗`活，可以適應(yīng)不規(guī)則的空間。有幾種選項(xiàng)可用于調(diào)整網(wǎng)格。

理想情況下，我們的目標(biāo)是擁有一個(gè)規(guī)則的網(wǎng)格，內(nèi)部有一層三角形，沒(méi)有聚集，并且在外層有一個(gè)平滑的、較低密度的三角形。

圖7
第三個(gè)網(wǎng)格似乎是最規(guī)則的，并且適合這個(gè)數(shù)據(jù)集。

圖8：這是要使用的最佳網(wǎng)格。
注意：您可以看到網(wǎng)格延伸到海岸線之外進(jìn)入海洋。由于我們?cè)噲D評(píng)估綠地比率對(duì)狐貍寄生蟲物種的影響，將屬于海洋的區(qū)域包含在網(wǎng)格中沒(méi)有意義。有兩種可能的解決方案：第一種是使用這個(gè)網(wǎng)格運(yùn)行模型，然后簡(jiǎn)單地忽略模型為海洋區(qū)域提供的結(jié)果。第二種是修改網(wǎng)格以反映海岸線。

5.2 投影矩陣（Projector matrix）

現(xiàn)在我們已經(jīng)構(gòu)建了我們的網(wǎng)格，我們需要將數(shù)據(jù)點(diǎn)與網(wǎng)格頂點(diǎn)相關(guān)聯(lián)。投影矩陣使用網(wǎng)格頂點(diǎn)作為明確的鄰居，為模型提供數(shù)據(jù)集的鄰域結(jié)構(gòu)。
如前所述，地統(tǒng)計(jì)數(shù)據(jù)沒(méi)有明確的鄰居，因此我們需要使用網(wǎng)格對(duì)空間進(jìn)行人工離散化。投影矩陣將點(diǎn)投影到網(wǎng)格上，其中每個(gè)頂點(diǎn)都有明確指定的鄰居。如果數(shù)據(jù)點(diǎn)落在頂點(diǎn)上（頂點(diǎn)是多邊形的每個(gè)角點(diǎn)，這里是三角形），那么它將直接與相鄰的頂點(diǎn)相關(guān)聯(lián)（就像圖中的藍(lán)色點(diǎn)）。然而，如果數(shù)據(jù)點(diǎn)落在網(wǎng)格三角形內(nèi)（深紅色點(diǎn)），它的權(quán)重將根據(jù)與每個(gè)頂點(diǎn)的接近程度在三個(gè)頂點(diǎn)之間分配（帶有深色邊框的紅色、橙色和黃色點(diǎn)）。然后，原始數(shù)據(jù)點(diǎn)將根據(jù)定義三角形的頂點(diǎn)的鄰居擁有更多的“偽鄰居”，權(quán)重的分配方式與這些頂點(diǎn)類似（但是，每個(gè)數(shù)據(jù)點(diǎn)的總權(quán)重始終為1）。

圖9：投影矩陣如何創(chuàng)建鄰居的圖形表示。
投影矩陣會(huì)自動(dòng)計(jì)算每個(gè)點(diǎn)的鄰域的權(quán)重向量，并通過(guò)將網(wǎng)格和數(shù)據(jù)點(diǎn)的位置提供給函數(shù)來(lái)計(jì)算。

5.3 SPDE

SPDE（隨機(jī)偏微分方程）是Matérn協(xié)方差函數(shù)的數(shù)學(xué)解，它實(shí)際上允許INLA有效地計(jì)算網(wǎng)格頂點(diǎn)處數(shù)據(jù)集的空間自相關(guān)結(jié)構(gòu)。它計(jì)算網(wǎng)格頂點(diǎn)之間的相關(guān)結(jié)構(gòu)（然后將通過(guò)使用投影矩陣計(jì)算的向量進(jìn)行加權(quán)，以計(jì)算適用于實(shí)際數(shù)據(jù)集的相關(guān)矩陣）。

5.4 擬合基本空間模型

我們將首先擬合一個(gè)僅包含截距和空間效應(yīng)的模型，以展示如何編寫此代碼。這個(gè)模型只是測(cè)試空間自相關(guān)對(duì)寄生蟲物種豐富度的影響，而不包括任何其他協(xié)變量。
需要記住的一件事是，INLA語(yǔ)法使用格式f(Covariate Name, model = Effect Type)對(duì)非線性效應(yīng)進(jìn)行編碼。對(duì)于空間效應(yīng)，模型名稱是您為SPDE指定的名稱（在這種情況下為spde1）。

#首先，我們指定公式
formula_p1?<-?y?~?-1?+?Intercept?+f(spatd1,?model?=?se1)?#?這指定了空間隨機(jī)效應(yīng)。名稱（sield1）由您選擇，但需要與您將在模型中包含的名稱相同

我們已經(jīng)有了公式，現(xiàn)在準(zhǔn)備運(yùn)行模型！

我們可以探索固定效應(yīng)和隨機(jī)效應(yīng)的結(jié)果。

#?我們可以通過(guò)使用以下方式訪問(wèn)固定（這里只有截距）和隨機(jī)效應(yīng)的摘要
round(Mod_Point1$summary.fixed,3)
round(Mod_Point1$summary.hyperpar\[1,\],3)

我們還可以使用emarginal()函數(shù)計(jì)算隨機(jī)項(xiàng)的方差（請(qǐng)記住，INLA使用精度，因此我們不能直接提取方差）。
注意：INLA提供了許多函數(shù)來(lái)處理后驗(yàn)邊緣。在本教程中，我們僅使用emarginal()（它計(jì)算函數(shù)的期望，并用于將精度轉(zhuǎn)換為方差等），但值得知道的是，有一整套用于邊緣操作的函數(shù)，例如從邊緣采樣、轉(zhuǎn)換它們或計(jì)算摘要統(tǒng)計(jì)信息。

現(xiàn)在回到提取我們的隨機(jī)項(xiàng)方差。

我們可以在此提取的兩個(gè)最重要的東西是范圍參數(shù)（kappa）、標(biāo)稱方差（sigma）和范圍（r，自相關(guān)降至0.1以下的半徑）。這些是空間自相關(guān)的重要參數(shù)：Kappa越高，空間自相關(guān)結(jié)構(gòu)越平滑（并且范圍越高）。較短的范圍表示緊密相鄰的點(diǎn)之間自相關(guān)的急劇增加和更強(qiáng)的自相關(guān)效應(yīng)。

5.5 構(gòu)建和運(yùn)行更復(fù)雜的空間模型

通常，我們感興趣的是擬合包含協(xié)變量的模型，并且關(guān)注這些協(xié)變量在考慮空間自相關(guān)的情況下如何影響響應(yīng)變量。在這種情況下，我們需要在模型構(gòu)建中添加另一個(gè)步驟。我們將保留之前使用的相同網(wǎng)格（Mesh3）和投影矩陣（A_point），并在此基礎(chǔ)上繼續(xù)進(jìn)行。我將簡(jiǎn)要提及對(duì)模型的各種自定義（例如時(shí)空建模）。
現(xiàn)在，我們將擴(kuò)展我們的模型以包含所有可用組件：

• 網(wǎng)格

• 投影矩陣

• 相關(guān)結(jié)構(gòu)說(shuō)明符（SPDE），包括對(duì)空間結(jié)構(gòu)的PC先驗(yàn)

• 空間索引

• 堆疊（Stack）

• 公式

5.5.1 指定PC先驗(yàn)

我們可以為空間項(xiàng)提供先驗(yàn)。一種特殊類型的先驗(yàn)（懲罰復(fù)雜度或PC先驗(yàn)）可以施加在SPDE上。這些先驗(yàn)被廣泛使用，因?yàn)樗鼈?#xff08;正如其名稱所示）懲罰模型的復(fù)雜性。實(shí)際上，它們使空間模型向基礎(chǔ)模型（沒(méi)有空間項(xiàng)的模型）收縮。為此，我們應(yīng)用弱信息先驗(yàn)，懲罰小范圍和大方差。有關(guān)PC先驗(yàn)如何工作的更詳細(xì)理論解釋，請(qǐng)查看Fulgstad等人（2018）的論文。

5.5.2 空間索引

一個(gè)有用的步驟包括構(gòu)建空間索引。這將為SPDE模型提供所有必需的元素。這并不是嚴(yán)格必要的，除非您想要?jiǎng)?chuàng)建多個(gè)空間場(chǎng)（例如特定年份的空間場(chǎng)）。重復(fù)的數(shù)量將產(chǎn)生iid（獨(dú)立同分布）的重復(fù)（方差將在各個(gè)級(jí)別上均勻分布，這相當(dāng)于GLM中的標(biāo)準(zhǔn)因子效應(yīng)），而組的數(shù)量將產(chǎn)生相關(guān)的重復(fù)（組的每個(gè)級(jí)別將依賴于前一個(gè)/后一個(gè)級(jí)別）。
以下顯示的是索引的默認(rèn)設(shè)置（未指定重復(fù)或組）：

5.5.3 堆疊

堆疊因其特別難以處理而聞名，但簡(jiǎn)而言之，它提供了將在模型中使用的所有元素。它包括數(shù)據(jù)、協(xié)變量（包括線性和非線性協(xié)變量）以及每個(gè)協(xié)變量的索引。需要記住的一件事是，堆疊不會(huì)自動(dòng)包含截距，因此需要明確指定。

注意：在這種情況下，截距被擬合為在空間中是常數(shù)（它與空間效應(yīng)一起擬合，這意味著在網(wǎng)格的每個(gè)n.spde頂點(diǎn)處它始終為1）。不一定是這種情況，如果您想擬合截距在整個(gè)數(shù)據(jù)集中是常數(shù)（并因此受到空間效應(yīng)的影響），您可以將其與其他協(xié)變量的列表一起編碼，但請(qǐng)記住，然后您需要將截距指定為Intercept = rep(1, n.dat)，其中n.dat是數(shù)據(jù)集中的數(shù)據(jù)點(diǎn)數(shù)量（而不是網(wǎng)格頂點(diǎn)的數(shù)量）。

5.5.4 擬合模型

在公式中，我們指定每個(gè)協(xié)變量應(yīng)該具有的效應(yīng)類型。線性變量以標(biāo)準(zhǔn)GLM方式指定，而隨機(jī)效應(yīng)和非線性效應(yīng)需要使用f(Cov Name, model = Effect Type)格式指定，類似于我們到目前為止看到的空間效應(yīng)項(xiàng)。

最后，我們準(zhǔn)備運(yùn)行模型。這包括堆疊（要包含哪些數(shù)據(jù)）、公式（如何對(duì)協(xié)變量進(jìn)行建模）以及模型的詳細(xì)信息（例如計(jì)算模型選擇工具或進(jìn)行預(yù)測(cè)）。此模型測(cè)試GS_ratio（綠地比率）和GS變異性對(duì)寄生蟲物種豐富度的影響，同時(shí)考慮空間自相關(guān)、時(shí)間自相關(guān)以及樣本采集的地點(diǎn)（以考慮重復(fù)采樣）。

現(xiàn)在我們可以制作一些圖來(lái)可視化我們感興趣的一些變量的影響。

綠地的數(shù)量（GS Ratio）顯然與物種豐富度呈正相關(guān)，但效果相當(dāng)線性，所以我們可能希望在下一個(gè)模型中考慮將其擬合為線性效應(yīng)（這樣做我們不會(huì)丟失太多信息）。

plot(Modummary.ranDate\[,1:2\],type?=?"l",lwd?=?2,

圖11：根據(jù)我們的模型結(jié)果可視化效應(yīng)
現(xiàn)在我們可以提取關(guān)于空間場(chǎng)的一些進(jìn)一步信息。

我們可能也有興趣可視化高斯隨機(jī)場(chǎng)（GRF）。如前所述，GRF表示在考慮模型中的所有協(xié)變量后，響應(yīng)變量在空間中的變化。它可以被視為“響應(yīng)變量在空間中的真實(shí)分布”。
然而，這也可能反映模型中缺少重要的協(xié)變量，檢查GRF的空間分布可以揭示缺少哪些協(xié)變量。例如，如果海拔與響應(yīng)變量呈正相關(guān)，但未包含在模型中，我們可能會(huì)在海拔較高的區(qū)域看到較高的后驗(yàn)均值。熟悉地形的研究人員將能夠識(shí)別這一點(diǎn)并相應(yīng)地改進(jìn)模型。

我們需要為GRF的均值和方差創(chuàng)建空間對(duì)象。

xmean_ras和xsd_ras是柵格項(xiàng)，可以使用writeRaster()函數(shù)在R之外（包括在GIS軟件中）導(dǎo)出、存儲(chǔ)和操作。
現(xiàn)在我們可以繪制GRF（我使用了與區(qū)域數(shù)據(jù)相同的配色方案）：

圖12：高斯隨機(jī)場(chǎng)的均值和方差

六、繪制空間預(yù)測(cè)和高斯隨機(jī)場(chǎng)

最后，我將展示如何從INLA模型生成空間預(yù)測(cè)。這將涉及一些柵格和矩陣的操作。本質(zhì)上，這歸結(jié)為創(chuàng)建一個(gè)我們沒(méi)有值但希望使用模型估計(jì)為響應(yīng)變量生成預(yù)測(cè)的空間坐標(biāo)網(wǎng)格（考慮數(shù)據(jù)的空間自相關(guān)結(jié)構(gòu)）。

圖13：綠地
為了使用INLA生成預(yù)測(cè)，我們需要生成一個(gè)數(shù)據(jù)集（在我們希望預(yù)測(cè)的位置附加坐標(biāo)），并為其附加一系列缺失的觀測(cè)值（在R中編碼為NA）。當(dāng)缺失的觀測(cè)值在響應(yīng)變量中時(shí)，INLA會(huì)自動(dòng)計(jì)算相應(yīng)線性預(yù)測(cè)器和擬合值的預(yù)測(cè)分布。
使用INLA語(yǔ)法，可以通過(guò)擬合一個(gè)響應(yīng)變量設(shè)置為NA的堆疊，然后將此堆疊與估計(jì)堆疊（與我們到目前為止使用的類似）連接來(lái)生成模型預(yù)測(cè)。然后，我們可以提取預(yù)測(cè)響應(yīng)變量的值，并使用inla.mjector()函數(shù)將這些值投影到網(wǎng)格頂點(diǎn)上（就像我們之前繪制GRF時(shí)所做的那樣）。
首先，我們將綠地?cái)?shù)量（GS ratio）的柵格值轉(zhuǎn)換為矩陣，然后將坐標(biāo)重新分配到一個(gè)ncol X nrow單元的矩陣中（列數(shù)和行數(shù)）。

接下來(lái)，我們需要?jiǎng)?chuàng)建一個(gè)ncol X nrow單元的網(wǎng)格，其中包含我們希望投影模型預(yù)測(cè)的點(diǎn)的坐標(biāo)。

Seqid?<-?seq(from?=?GS_Prnt@xmin,to?=?GS_Pred

現(xiàn)在我們已經(jīng)有了預(yù)測(cè)的坐標(biāo)網(wǎng)格，接下來(lái)要?jiǎng)?chuàng)建一個(gè)新的堆疊用于預(yù)測(cè)，這個(gè)堆疊將包含預(yù)測(cè)點(diǎn)的信息以及與估計(jì)模型相同的結(jié)構(gòu)（協(xié)變量等）。

現(xiàn)在我們可以使用連接后的堆疊來(lái)運(yùn)行模型并獲取預(yù)測(cè)結(jié)果。

接下來(lái)，我們將預(yù)測(cè)的均值和標(biāo)準(zhǔn)差轉(zhuǎn)換為柵格對(duì)象，以便進(jìn)行可視化。

現(xiàn)在我們可以繪制預(yù)測(cè)的均值和標(biāo)準(zhǔn)差的柵格圖，以直觀地展示寄生蟲物種豐富度的預(yù)測(cè)情況。

#?繪制預(yù)測(cè)均值的柵格圖
par(mfrow?=?c(1,?1),?mar?=?c(2,?2,?1,?1))
plot(prean

通過(guò)以上步驟，我們完成了從構(gòu)建模型到生成預(yù)測(cè)以及可視化預(yù)測(cè)結(jié)果的整個(gè)過(guò)程。這些預(yù)測(cè)結(jié)果可以幫助我們更好地理解在不同空間位置上寄生蟲物種豐富度的可能情況，同時(shí)考慮了空間自相關(guān)和其他協(xié)變量的影響。

七、結(jié)論

在本文中，我們深入探討了使用INLA（集成嵌套拉普拉斯近似）在R中進(jìn)行空間數(shù)據(jù)建模的過(guò)程，特別是針對(duì)地統(tǒng)計(jì)（標(biāo)記點(diǎn)）數(shù)據(jù)。我們從數(shù)據(jù)加載和準(zhǔn)備開始，詳細(xì)介紹了如何處理坐標(biāo)參考系統(tǒng)，確保數(shù)據(jù)的正確投影和可視化。
我們構(gòu)建了空間模型的基本組件，包括網(wǎng)格、投影矩陣和SPDE，通過(guò)逐步演示，展示了如何根據(jù)數(shù)據(jù)特點(diǎn)選擇合適的參數(shù)來(lái)構(gòu)建這些組件。在模型擬合方面，我們先從一個(gè)簡(jiǎn)單的僅包含截距和空間效應(yīng)的基礎(chǔ)模型入手，逐步擴(kuò)展到包含多個(gè)協(xié)變量的復(fù)雜模型，包括如何指定PC先驗(yàn)、空間索引和堆疊等關(guān)鍵步驟。
在模型結(jié)果的分析和可視化方面，我們不僅展示了如何提取固定效應(yīng)和隨機(jī)效應(yīng)的摘要信息，還介紹了如何計(jì)算隨機(jī)項(xiàng)的方差，以及如何可視化非線性效應(yīng)、高斯隨機(jī)場(chǎng)和空間預(yù)測(cè)結(jié)果。這些可視化工具幫助我們直觀地理解模型的輸出，發(fā)現(xiàn)變量之間的關(guān)系以及空間上的變化模式。

本文中分析的完整數(shù)據(jù)、代碼、文檔分享到會(huì)員群，掃描下面二維碼即可加群！?

---

資料獲取

在公眾號(hào)后臺(tái)回復(fù)“領(lǐng)資料”，可免費(fèi)獲取數(shù)據(jù)分析、機(jī)器學(xué)習(xí)、深度學(xué)習(xí)等學(xué)習(xí)資料。

點(diǎn)擊文末“閱讀原文”

獲取完整代碼、數(shù)據(jù)、文檔。

本文選自《R-INLA實(shí)現(xiàn)綠地與狐貍寄生蟲數(shù)據(jù)空間建模：含BYM、SPDE模型及PC先驗(yàn)應(yīng)用可視化》。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語(yǔ)言貝葉斯INLA空間自相關(guān)、混合效應(yīng)、季節(jié)空間模型、SPDE、時(shí)空分析野生動(dòng)物數(shù)據(jù)可視化

R語(yǔ)言貝葉斯廣義線性混合（多層次/水平/嵌套）模型GLMM、邏輯回歸分析教育留級(jí)影響因素?cái)?shù)據(jù)

R語(yǔ)言線性混合效應(yīng)模型（固定效應(yīng)&隨機(jī)效應(yīng)）和交互可視化3案例

用SPSS估計(jì)HLM多層（層次）線性模型模型

R語(yǔ)言用線性混合效應(yīng)（多水平/層次/嵌套）模型分析聲調(diào)高低與禮貌態(tài)度的關(guān)系

R語(yǔ)言LME4混合效應(yīng)模型研究教師的受歡迎程度

R語(yǔ)言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實(shí)例

R語(yǔ)言混合線性模型、多層次模型、回歸模型分析學(xué)生平均成績(jī)GPA和可視化

R語(yǔ)言線性混合效應(yīng)模型（固定效應(yīng)&隨機(jī)效應(yīng)）和交互可視化3案例

R語(yǔ)言用lme4多層次（混合效應(yīng)）廣義線性模型（GLM），邏輯回歸分析教育留級(jí)調(diào)查數(shù)據(jù)

R語(yǔ)言 線性混合效應(yīng)模型實(shí)戰(zhàn)案例

R語(yǔ)言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)

R語(yǔ)言如何用潛類別混合效應(yīng)模型（LCMM）分析抑郁癥狀

R語(yǔ)言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究

R語(yǔ)言建立和可視化混合效應(yīng)模型mixed effect model

R語(yǔ)言LME4混合效應(yīng)模型研究教師的受歡迎程度

R語(yǔ)言 線性混合效應(yīng)模型實(shí)戰(zhàn)案例

R語(yǔ)言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）

R語(yǔ)言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究

R語(yǔ)言如何解決線性混合模型中畸形擬合(Singular fit)的問(wèn)題

基于R語(yǔ)言的lmer混合線性回歸模型

R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型

R語(yǔ)言分層線性模型案例

R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)（SAT）建立分層模型

使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM

R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型

SPSS中的多層（等級(jí)）線性模型Multilevel linear models研究整容手術(shù)數(shù)據(jù)

用SPSS估計(jì)HLM多層（層次）線性模型模型