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

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

關(guān)鍵詞

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

一、引言

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

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

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

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

4. 繪制空間預(yù)測和高斯隨機場。

二、所需知識和軟件包

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

#?添加dep?=?T表示您還將安裝軟件包依賴項
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ù)集

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

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

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

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

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

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

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

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

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

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

4.2 鄰接矩陣計算

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

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

圖2：鄰接矩陣

4.3 模型公式指定

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

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

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

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

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

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

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

4.5 設(shè)置先驗

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

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

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

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

spplot(obj?=?Fost,?zcol

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

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

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

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

• JanDate（收集樣本的日期）

• Site（從哪個公園收集的）

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

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

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

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

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

嘗試在同一地圖上繪制我們的點和shapefile文件將不起作用，因為它們的坐標(biāo)表示在不同的系統(tǒng)中，無法直接繪制。

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

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

5.1 模型的基本組件

模型的絕對必要組件包括：

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

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

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

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

5.2 投影矩陣（Projector matrix）

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

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

5.3 SPDE

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

5.4 擬合基本空間模型

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

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

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

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

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

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

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

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

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

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

• 網(wǎng)格

• 投影矩陣

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

• 空間索引

• 堆疊（Stack）

• 公式

5.5.1 指定PC先驗

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

5.5.2 空間索引

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

5.5.3 堆疊

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

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

5.5.4 擬合模型

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

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

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

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

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

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

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

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

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

圖12：高斯隨機場的均值和方差

六、繪制空間預(yù)測和高斯隨機場

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

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

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

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

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

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

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

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

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

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

七、結(jié)論

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

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

---

資料獲取

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

點擊文末“閱讀原文”

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

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

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

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

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

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

用SPSS估計HLM多層（層次）線性模型模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

基于R語言的lmer混合線性回歸模型

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

R語言分層線性模型案例

R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗（SAT）建立分層模型

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

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

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

用SPSS估計HLM多層（層次）線性模型模型