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

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

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

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

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

• Cases：1978-1982年期間的白血病病例數(shù)。

• POP8：1980年人口。

• PCTOWNHOME：擁有房屋的人口比例。

• PCTAGE65P：65歲以上的人口比例。

• AVGIDIST：到最近的三氯乙烯（TCE）站點的平均反距離。

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

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

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

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

疾病作圖

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

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

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

library(viridis)##?Loading?required?package:?viridisLitespplot(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污染的站點的位置，可以通過縮小看到它。

點擊標題查閱往期相關(guān)內(nèi)容

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

左右滑動查看更多

01

02

03

04

混合效應(yīng)模型

泊松回歸

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

模型 ：

請注意，它的glm功能類似于該功能。在此，參數(shù)?E用于預(yù)期的案例數(shù)?；? 設(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

具有隨機效應(yīng)的泊松回歸

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

現(xiàn)在，該模式的摘要包括有關(guā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ū)域?qū)@示相似的目標變量值。

鄰接矩陣

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

這將返回一個nb具有鄰域結(jié)構(gòu)定義的對象：

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

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

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

回歸模型

通常情況是，除了\（y_i \）之外，我們還有許多協(xié)變量 \（X_i \）。因此，我們可能想對\（X_i \）回歸?\（y_i \）。除了 協(xié)變量，我們可能還需要考慮數(shù)據(jù)的空間結(jié)構(gòu)。
可以使用不同類型的回歸模型來建模晶格數(shù)據(jù)：

• 廣義線性模型（具有空間隨機效應(yīng)）。

• 空間計量經(jīng)濟學(xué)模型。

線性混合模型

一種常見的方法（對于高斯數(shù)據(jù)）是使用
具有隨機效應(yīng)的線性回歸：

\ [
Y = X \ beta + Zu + \ varepsilon
]

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

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

\（\ Sigma \）的定義是，它會引起與相鄰區(qū)域的更高相關(guān)性，\（Z \）是隨機效果的設(shè)計矩陣，而
\（\ varepsilon_i \ sim N（0，\ sigma ^ 2），i = 1，\ ldots，n \）是一個誤差項。

空間隨機效應(yīng)的結(jié)構(gòu)

在\（\ 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ū)域中的隨機效應(yīng)，模型的類型 和鄰接矩陣。為此，將使用稀疏矩陣。

##?##?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é模型包括兩個潛在的隨機效應(yīng)：ICAR 潛在效應(yīng)和高斯iid潛在效應(yīng)。線性預(yù)測變量\（\ eta_i \）
為：

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

• \（u_i \）是iid高斯隨機效應(yīng)

• \（v_i \）是內(nèi)在的CAR隨機效應(yīng)

##?##?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等人的模型）
定義的，以定義潛在效應(yīng)的精度矩陣：

\ [
\ Sigma ^ {-1} = [（1-\ lambda）I_n + \ lambda M]; \ \ lambda \ in（0,1）
]

為了定義正確的模型，我們應(yīng)采用矩陣\（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)濟學(xué)模型

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

同步自回歸模型（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īng)的自回歸：

\ [
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 \）是空間自相關(guān)參數(shù)。

• \（X \）是協(xié)變量的矩陣，系數(shù)為\（\ beta \）。

• \（e \）是具有方差\（\ sigma ^ 2 \）的高斯iid誤差。

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

模型定義

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

• 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è)置了精度參數(shù)\（\ tau \）和空間 自相關(guān)參數(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 \）的伽瑪先驗值，而 為空間自相關(guān)參數(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ù)的估計顯示為隨機效應(yīng)的一部分：

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

空間自相關(guān)以內(nèi)部比例報告（即 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")

結(jié)果匯總

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)生相對風(fēng)險的更平滑的估計。

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

參考文獻

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ù)資料。

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

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的貝葉斯分層混合模型的診斷準確性研究
R語言建立和可視化混合效應(yīng)模型mixed effect model
R語言LME4混合效應(yīng)模型研究教師的受歡迎程度
R語言 線性混合效應(yīng)模型實戰(zhàn)案例
R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
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多層（層次）線性模型模型