拓端tecdat|R語言使用貝葉斯層次模型進行空間數據分析

本文轉載自查看原文 2020-02-17 17:59 788

原文鏈接：http://tecdat.cn/?p=10932

介紹

在本節中，我將重點介紹使用集成嵌套拉普拉斯近似方法的貝葉斯推理。
可以估計貝葉斯層次模型的后邊緣分布。鑒於模型類型非常廣泛，我們將重點關注用於分析晶格數據的空間模型。

數據集：紐約州北部的白血病

為了說明如何與空間模型擬合，將使用紐約白血病數據集。該數據集記錄了普查區紐約州北部的許多白血病病例。數據集中的一些變量是：

Cases：1978-1982年期間的白血病病例數。
POP8：1980年人口。
PCTOWNHOME：擁有房屋的人口比例。
PCTAGE65P：65歲以上的人口比例。
AVGIDIST：到最近的三氯乙烯（TCE）站點的平均反距離。

鑒於有興趣研究紐約州北部的白血病風險，因此首先要計算預期的病例數。這是通過計算總死亡率（總病例數除以總人口數）並將其乘以總人口數得出的：

rate <- sum(NY8$Cases) / sum(NY8$POP8)
NY8$Expected <- NY8$POP8 * rate

一旦獲得了預期的病例數，就可以使用標准化死亡率（SMR）來獲得原始的風險估計，該標准是將觀察到的病例數除以預期的病例數得出的：

NY8$SMR <- NY8$Cases / NY8$Expected

疾病作圖

在流行病學中，重要的是制作地圖以顯示相對風險的空間分布。在此示例中，我們將重點放在錫拉庫扎市以減少生成地圖的計算時間。因此，我們用錫拉丘茲市的區域創建索引：

# Subset Syracuse city
syracuse <- which(NY8$AREANAME == "Syracuse city")

可以使用函數spplot（在包中sp）簡單地創建疾病圖：

library(viridis)
## Loading required package: viridisLite
spplot(NY8[syracuse, ], "SMR", #at = c(0.6, 0.9801, 1.055, 1.087, 1.125, 13),
   col.regions = rev(magma(16))) #gray.colors(16, 0.9, 0.4))

## Loading required package: viridisLite

可以輕松創建交互式地圖

請注意，先前的地圖還包括11個受TCE污染的站點的位置，可以通過縮小看到它。

混合效應模型

泊松回歸

我們將考慮的第一個模型是沒有潛在隨機效應的Poisson模型，因為這將提供與其他模型進行比較的基准。

模型：

請注意，它的glm功能類似於該功能。在此，參數 E用於預期的案例數。或設置了其他參數來計算模型參數的邊際
（使用control.predictor）並計算一些模型選擇標准（使用control.compute）。

接下來，可以獲得模型的摘要：

summary(m1)

## 
## Call:

## Time used:
##     Pre = 0.368, Running = 0.0968, Post = 0.0587, Total = 0.524 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.065 0.045     -0.155   -0.065      0.023 -0.064   0
## AVGIDIST     0.320 0.078      0.160    0.322      0.465  0.327   0
## 
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 140.25 
## 
## Deviance Information Criterion (DIC) ...............: 948.12
## Deviance Information Criterion (DIC, saturated) ....: 418.75
## Effective number of parameters .....................: 2.00
## 
## Watanabe-Akaike information criterion (WAIC) ...: 949.03
## Effective number of parameters .................: 2.67
## 
## Marginal log-Likelihood:  -480.28 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

具有隨機效應的泊松回歸

可以通過在線性預測變量中包括iid高斯隨機效應，將潛在隨機效應添加到模型中，以解決過度分散問題。

現在，該模式的摘要包括有關隨機效果的信息：

summary(m2)

## 
## Call:

## Time used:
##     Pre = 0.236, Running = 0.315, Post = 0.0744, Total = 0.625 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.126 0.064     -0.256   -0.125     -0.006 -0.122   0
## AVGIDIST     0.347 0.105      0.139    0.346      0.558  0.344   0
## 
## Random effects:
##   Name     Model
##     ID IID model
## 
## Model hyperparameters:
##                     mean       sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID 3712.34 11263.70       3.52     6.94   39903.61 5.18
## 
## Expected number of effective parameters(stdev): 54.95(30.20)
## Number of equivalent replicates : 5.11 
## 
## Deviance Information Criterion (DIC) ...............: 926.93
## Deviance Information Criterion (DIC, saturated) ....: 397.56
## Effective number of parameters .....................: 61.52
## 
## Watanabe-Akaike information criterion (WAIC) ...: 932.63
## Effective number of parameters .................: 57.92
## 
## Marginal log-Likelihood:  -478.93 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

添加點估計以進行映射

這兩個模型估計可以被添加到 SpatialPolygonsDataFrame NY8

NY8$FIXED.EFF <- m1$summary.fitted[, "mean"]
NY8$IID.EFF <- m2$summary.fitted[, "mean"]
spplot(NY8[syracuse, ], c("SMR", "FIXED.EFF", "IID.EFF"),
  col.regions = rev(magma(16)))

晶格數據的空間模型

格子數據涉及在不同區域（例如，鄰里，城市，省，州等）測量的數據。出現空間依賴性是因為相鄰區域將顯示相似的目標變量值。

鄰接矩陣

可以使用poly2nbpackage中的函數來計算鄰接矩陣 spdep。如果其邊界至少在某一點上接觸，則此功能會將兩個區域視為鄰居：

這將返回一個nb具有鄰域結構定義的對象：

NY8.nb

## Neighbour list object:
## Number of regions: 281 
## Number of nonzero links: 1624 
## Percentage nonzero weights: 2.056712 
## Average number of links: 5.779359

另外，當多邊形的重心已知時，可以繪制對象：

plot(NY8) 
plot(NY8.nb, coordinates(NY8), add = TRUE, pch = ".", col = "gray")

回歸模型

通常情況是，除了\（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-函數定義空間潛在效果。這將需要一個索引來識別每個區域中的隨機效應，模型的類型和鄰接矩陣。為此，將使用稀疏矩陣。

## 
## Call:

## Time used:
##     Pre = 0.298, Running = 0.305, Post = 0.0406, Total = 0.644 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.122 0.052     -0.226   -0.122     -0.022 -0.120   0
## AVGIDIST     0.318 0.121      0.075    0.320      0.551  0.324   0
## 
## Random effects:
##   Name     Model
##     ID Besags ICAR model
## 
## Model hyperparameters:
##                  mean   sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID 3.20 1.67       1.41     2.79       7.56 2.27
## 
## Expected number of effective parameters(stdev): 46.68(12.67)
## Number of equivalent replicates : 6.02 
## 
## Deviance Information Criterion (DIC) ...............: 904.12
## Deviance Information Criterion (DIC, saturated) ....: 374.75
## Effective number of parameters .....................: 48.31
## 
## Watanabe-Akaike information criterion (WAIC) ...: 906.77
## Effective number of parameters .................: 44.27
## 
## Marginal log-Likelihood:  -685.70 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

BYM模型

Besag，York和Mollié模型包括兩個潛在的隨機效應：ICAR 潛在效應和高斯iid潛在效應。線性預測變量\（\ eta_i \）
為：

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

\（u_i \）是iid高斯隨機效應
\（v_i \）是內在的CAR隨機效應

## 
## Call:

## Time used:
##     Pre = 0.294, Running = 1, Post = 0.0652, Total = 1.36 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.123 0.052     -0.227   -0.122     -0.023 -0.121   0
## AVGIDIST     0.318 0.121      0.075    0.320      0.551  0.324   0
## 
## Random effects:
##   Name     Model
##     ID BYM model
## 
## Model hyperparameters:
##                                         mean      sd 0.025quant 0.5quant
## Precision for ID (iid component)     1790.06 1769.02     115.76  1265.24
## Precision for ID (spatial component)    3.12    1.36       1.37     2.82
##                                      0.975quant   mode
## Precision for ID (iid component)        6522.28 310.73
## Precision for ID (spatial component)       6.58   2.33
## 
## Expected number of effective parameters(stdev): 47.66(12.79)
## Number of equivalent replicates : 5.90 
## 
## Deviance Information Criterion (DIC) ...............: 903.41
## Deviance Information Criterion (DIC, saturated) ....: 374.04
## Effective number of parameters .....................: 48.75
## 
## Watanabe-Akaike information criterion (WAIC) ...: 906.61
## Effective number of parameters .................: 45.04
## 
## Marginal log-Likelihood:  -425.65 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

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} \）是一個：

max(eigen(Cmatrix)$values)
## [1] 1

## [1] 1

該模型與往常一樣具有功能inla。注意，\（C \）矩陣使用參數
傳遞給f函數Cmatrix：

## 
## Call:
## Time used:
##     Pre = 0.236, Running = 0.695, Post = 0.0493, Total = 0.98 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode   kld
## (Intercept) -0.128 0.448      -0.91   -0.128      0.656 -0.126 0.075
## AVGIDIST     0.325 0.122       0.08    0.327      0.561  0.330 0.000
## 
## Random effects:
##   Name     Model
##     ID Generic1 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ID 2.720 1.098       1.27    2.489      5.480 2.106
## Beta for ID      0.865 0.142       0.47    0.915      0.997 0.996
## 
## Expected number of effective parameters(stdev): 52.25(13.87)
## Number of equivalent replicates : 5.38 
## 
## Deviance Information Criterion (DIC) ...............: 903.14
## Deviance Information Criterion (DIC, saturated) ....: 373.77
## Effective number of parameters .....................: 53.12
## 
## Watanabe-Akaike information criterion (WAIC) ...: 906.20
## Effective number of parameters .................: 48.19
## 
## Marginal log-Likelihood:  -474.94 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

空間計量經濟學模型

空間計量經濟學是通過對空間建模略有不同的方法開發的。除了使用潛在效應，還可以對空間依賴性進行顯式建模。

同步自回歸模型（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。在這里，我們創建了一個具有相同名稱的列表，以將所有必需的值保存在一起：

#Arguments for 'slm'
args.slm = list(
   rho.min = rho.min ,
   rho.max = rho.max,
   W = W,
   X = mmatrix,
   Q.beta = Q.beta
)

此外，還設置了精度參數\（\ tau \）和空間自相關參數\（\ rho \）的先驗：

#Prior on rho
hyper.slm = list(
   prec = list(
      prior = "loggamma", param = c(0.01, 0.01)),
      rho = list(initial=0, prior = "logitbeta", param = c(1,1))
)

先前的定義使用具有不同參數的命名列表。參數 prior定義了使用之前param及其參數。在此，為精度分配了帶有參數\（0.01 \）和\（0.01 \）的伽瑪先驗值，而為空間自相關參數指定了帶有參數\（1 \）和\（1 \）的beta先驗值（即a間隔\（（（1，1）\））中的均勻先驗。

模型擬合

## 
## Call:
## Time used:
##     Pre = 0.265, Running = 1.2, Post = 0.058, Total = 1.52 
## Random effects:
##   Name     Model
##     ID SLM model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ID 8.989 4.115      3.709    8.085     19.449 6.641
## Rho for ID       0.822 0.073      0.653    0.832      0.936 0.854
## 
## Expected number of effective parameters(stdev): 62.82(15.46)
## Number of equivalent replicates : 4.47 
## 
## Deviance Information Criterion (DIC) ...............: 902.32
## Deviance Information Criterion (DIC, saturated) ....: 372.95
## Effective number of parameters .....................: 64.13
## 
## Watanabe-Akaike information criterion (WAIC) ...: 905.19
## Effective number of parameters .................: 56.12
## 
## Marginal log-Likelihood:  -477.30 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

系數的估計顯示為隨機效應的一部分：

round(m.slm$summary.random$ID[47:48,], 4)

##    ID   mean     sd 0.025quant 0.5quant 0.975quant   mode kld
## 47 47 0.4634 0.3107    -0.1618   0.4671     1.0648 0.4726   0
## 48 48 0.2606 0.3410    -0.4583   0.2764     0.8894 0.3063   0

空間自相關以內部比例報告（即 0到1 之間），並且需要重新縮放：

## Mean            0.644436 
## Stdev           0.145461 
## Quantile  0.025 0.309507 
## Quantile  0.25  0.556851 
## Quantile  0.5   0.663068 
## Quantile  0.75  0.752368 
## Quantile  0.975 0.869702

plot(marg.rho, type = "l", main = "Spatial autocorrelation")

結果匯總

NY8$ICAR <- m.icar$summary.fitted.values[, "mean"]
NY8$BYM <- m.bym$summary.fitted.values[, "mean"]
NY8$LEROUX <- m.ler$summary.fitted.values[, "mean"]
NY8$SLM <- m.slm$summary.fitted.values[, "mean"]
spplot(NY8[syracuse, ], 
  c("FIXED.EFF", "IID.EFF", "ICAR", "BYM", "LEROUX", "SLM"),
  col.regions = rev(magma(16))
)

注意空間模型如何產生相對風險的更平滑的估計。

為了選擇最佳模型，可以使用上面計算的模型選擇標准：

參考文獻

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

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 拓端tecdat：R語言用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預測工人工資關於《空間數據分析與R語言實踐》數據分析與挖掘 - R語言：貝葉斯分類算法（案例二）數據分析與挖掘 - R語言：貝葉斯分類算法（案例三）拓端數據tecdat|R語言使用多元AR-GARCH模型衡量市場風險 GIS Experience （十）：Python地理空間數據分析拓端數據tecdat|R語言中的BP神經網絡模型分析學生成績拓端tecdat|R語言POT超閾值模型和極值理論分析 R語言地理空間分析（一）讀寫空間數據文件文件拓端數據tecdat|R語言基於線性回歸的資本資產定價模型（CAPM）