原標題：R語言用貝葉斯層次模型進行空間數(shù)據(jù)分析|附代碼數(shù)據(jù)

閱讀全文：http://tecdat.cn/?p=10932

最近我們被客戶要求撰寫關于貝葉斯層次模型的研究報告，包括一些圖形和統(tǒng)計輸出。

在本文中，我將重點介紹使用集成嵌套 拉普拉斯近似方法的貝葉斯推理?？梢怨烙嬝惾~斯 層次模型的后邊緣分布。鑒于模型類型非常廣泛，我們將重點關注用于分析晶格數(shù)據(jù)的空間模型

數(shù)據(jù)集：紐約州北部的白血病

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

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

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

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

NY8$Expected <- NY8$POP8 * rate

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

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

疾病作圖

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

# Subset Syracuse city

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

可以使用函數(shù)spplot（在包中sp）簡單地創(chuàng)建疾病圖：

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

可以輕松創(chuàng)建交互式地圖

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

點擊標題查閱往期相關內容

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

左右滑動查看更多

01

02

03

04

混合效應模型

泊松回歸

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

模型 ：

請注意，它的glm功能類似于該功能。在此，參數(shù) E用于預期的案例數(shù)?；?設置了其他參數(shù)來計算模型參數(shù)的邊際

（使用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高斯隨機效應，將潛在隨機效應添加到模型中，以解決過度分散問題。

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

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)))

晶格數(shù)據(jù)的空間模型

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

鄰接矩陣

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

這將返回一個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 \）之外，我們還有許多協(xié)變量 \（X_i \）。因此，我們可能想對\（X_i \）回歸 \（y_i \）。除了 協(xié)變量，我們可能還需要考慮數(shù)據(jù)的空間結構。

可以使用不同類型的回歸模型來建模晶格數(shù)據(jù)：

廣義線性模型（具有空間隨機效應）。 空間計量經(jīng)濟學模型。

線性混合模型

一種常見的方法（對于高斯數(shù)據(jù)）是使用

具有隨機效應的線性回歸：

\ [

Y = X \ beta + Zu + \ varepsilon

]

隨機效應的向量\（u \）被建模為多元正態(tài)分布：

\ [

u \ sim N（0，\ sigma ^ 2_u \ Sigma）

]

\（\ Sigma \）的定義是，它會引起與相鄰區(qū)域的更高相關性，\（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 \）是區(qū)域\（i \）的鄰居數(shù)。 \（\ Sigma_ {i，j} \）取決于\（d（i，j）\）的函數(shù)。例如：

\ [

\ Sigma_ {i，j} = \ exp \ {-d（i，j）/ \ phi }

]

矩陣的“混合”（Leroux等人的模型）： \ [ \ Sigma = [（1 – \ lambda）I_n + \ lambda M] ^ {-1}; \ \ lambda \ in（0,1） ]

ICAR模型

第一個示例將基于ICAR規(guī)范。請注意， 使用f-函數(shù)定義空間潛在效果。這將需要 一個索引來識別每個區(qū)域中的隨機效應，模型的類型 和鄰接矩陣。為此，將使用稀疏矩陣。

##

## 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）

]

為了擬合模型，第一步是創(chuàng)建矩陣\（M \）：

我們可以檢查最大特征值\（\ lambda_ {max} \）是一個：

max(eigen(Cmatrix)$values)

## [1] 1

## [1] 1

該模型與往常一樣具有功能inla。注意，\（C \）矩陣使用參數(shù)

傳遞給f函數(shù)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

空間計量經(jīng)濟學模型

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

同步自回歸模型（SEM）

該模型包括協(xié)變量和誤差項的自回歸：

\ [

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）

該模型包括協(xié)變量和響應的自回歸：

\ [

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

現(xiàn)在包括一個_實驗_所謂的新的潛在影響slm，以 符合以下模型：

\ [

\ mathbf {x} =（I_n- \ rho W）^ {-1}（X \ beta + e）

]

該模型的元素是：

\（W \）是行標準化的鄰接矩陣。 \（\ rho \）是空間自相關參數(shù)。 \（X \）是協(xié)變量的矩陣，系數(shù)為\（\ beta \）。 \（e \）是具有方差\（\ sigma ^ 2 \）的高斯iid誤差。

該slm潛效果的實驗，它可以 與所述線性預測其他效果組合。

模型定義

為了定義模型，我們需要：

X：協(xié)變量矩陣 W：行標準化的鄰接矩陣 Q：系數(shù)\（\ beta \）的精確矩陣 范圍\（\ RHO \） ，通常由本征值定義

slm潛在作用是通過參數(shù)傳遞 args.sm。在這里，我們創(chuàng)建了一個具有相同名稱的列表，以將 所有必需的值保存在一起：

#Arguments for slm

args.slm = list(

rho.min = rho.min ,

rho.max = rho.max,

W = W,

X = mmatrix,

Q.beta = Q.beta

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

#rho的先驗

hyper.slm = list(

prec = list(

prior = "loggamma", param = c(0.01, 0.01)),

rho = list(initial=0, prior = "logitbeta", param = c(1,1))

先前的定義使用具有不同參數(shù)的命名列表。參數(shù) prior定義了使用之前param及其參數(shù)。在此，為 精度分配了帶有參數(shù)\（0.01 \）和\（0.01 \）的伽瑪先驗值，而 為空間自相關參數(shù)指定了帶有參數(shù)\（1 \） 和\（1 \）的beta先驗值（即a區(qū)間\（（（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

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

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))

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

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

參考文獻

Bivand, R., E. Pebesma and V. Gómez-Rubio (2013). Applied spatial data

analysis with R. Springer-Verlag. New York.

本文摘選 《 R語言使用貝葉斯層次模型進行空間數(shù)據(jù)分析 》 ，點擊“閱讀原文”獲取全文完整代碼數(shù)據(jù)資料。

點擊標題查閱往期內容

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

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

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

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

R語言LME4混合效應模型研究教師的受歡迎程度R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實例

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

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

R語言用lme4多層次（混合效應）廣義線性模型（GLM），邏輯回歸分析教育留級調查數(shù)據(jù)R語言 線性混合效應模型實戰(zhàn)案例

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

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

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

R語言建立和可視化混合效應模型mixed effect model

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

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

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

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

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

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

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

R語言分層線性模型案例

R語言用WinBUGS 軟件對學術能力測驗（SAT）建立分層模型

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

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

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

用SPSS估計HLM多層（層次）線性模型模型返回搜狐，查看更多

責任編輯：