貝葉斯地理統(tǒng)計模型R-INLA-4 貝葉斯時空模型
在前述的內(nèi)容中足绅，我們介紹了准夷，如何處理空間的數(shù)據(jù)球碉，利用海拔高度預測降雨量的例子。但是該例子僅僅涉及到的是涉及到回歸方程中，考慮影響因素及空間效應。
那么如果我們的數(shù)據(jù)有時間信息庶诡，如何加入到貝葉斯時空分析呢。譬如每年對某一個地區(qū)進行疾病的發(fā)病率調(diào)查咆课，10年數(shù)據(jù)整合在一起末誓，就可以從時間上或空間上看疾病的變化規(guī)律扯俱，也就會用到貝葉斯時空模型。
下面我們將介紹貝葉斯時空模型基显。該文章中蘸吓，會簡化數(shù)學計算的過程，主要是針對撩幽，在有數(shù)據(jù)的基礎上库继，如何應用貝葉斯時空模型，找出影響因素窜醉，繪制時間變化的空間分布預測圖宪萄。
1.數(shù)據(jù)集
我們模擬一個房屋與面積的地理位置數(shù)據(jù)，從2010年到2015年的房地產(chǎn)數(shù)據(jù)榨惰。因為是模擬數(shù)據(jù)拜英，不具有任何實際意義，僅僅作為展示琅催。
通過簡單的回歸方程居凶，發(fā)現(xiàn)，房屋價格與面積及年份成正相關(guān)藤抡，具有統(tǒng)計學意義侠碧。說明隨時間推遲，房子越值錢缠黍，且面積越大價格也越高弄兜。符合實際
library(inlabru)
library(INLA)
library(AmesHousing)
rm(list = ls())
library(INLA)
data("PRborder")
data(PRprec)
coords <- as.matrix(PRprec[sample(1:nrow(PRprec)), 1:2])
set.seed(1)
df=as.data.frame(coords)
df1= df %>% mutate(area=rnorm(nrow(coords),120,25),
price=area*5+rnorm(nrow(coords),10,5),
year=2010)%>% as.tbl()
df=as.data.frame(coords)
df2= df %>% mutate(area=rnorm(nrow(coords),120,25),
price=area*4+rnorm(nrow(coords),10,5),
year=2011) %>% as.tbl()
df=as.data.frame(coords)
df3= df %>% mutate(area=rnorm(nrow(coords),120,25),
price=area*3+rnorm(nrow(coords),16,5),
year=2012) %>% as.tbl()
df=as.data.frame(coords)
df4= df %>% mutate(area=rnorm(nrow(coords),120,25),
price=area*8+rnorm(nrow(coords),13,5),
year=2013) %>% as.tbl()
df=as.data.frame(coords)
df5= df %>% mutate(area=rnorm(nrow(coords),120,25),
price=area*7+rnorm(nrow(coords),11,5),
year=2014) %>% as.tbl()
df6= df %>% mutate(area=rnorm(nrow(coords),121,15),
price=area*5+rnorm(nrow(coords),12,5),
year=2015)%>% as.tbl()
df=rbind(df1,df2,df3,df4,df5,df6)
# plot
ggplot(df, aes(x=Latitude,y=Longitude)) +
geom_point(aes(color = price)) +
facet_wrap(~year)
# glm
fit_glm=glm(price~area+year,data = df)
summary(fit_glm)
glm(formula = price ~ area + year, data = df)
Deviance Residuals:
Min 1Q Median 3Q Max
-398.46 -158.55 -12.16 122.23 493.73
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -96809.399 3697.070 -26.18 <2e-16 ***
area 5.019 0.131 38.32 <2e-16 ***
year 48.128 1.837 26.20 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 36378.97)
Null deviance: 212678269 on 3695 degrees of freedom
Residual deviance: 134347545 on 3693 degrees of freedom
AIC: 49308
Number of Fisher Scoring iterations: 2
2.INLA
INLA 的步驟包括產(chǎn)生Mesh，建立A project瓷式，將空間位置與三角形網(wǎng)格點對應替饿，計算spde的空間效應。然后整合data到stack里面贸典。這是建立INLA的關(guān)鍵视卢，最后，寫INLA的公式瓤漏，帶入INLA模型腾夯。
2.1 Mesh
下面我們利用時空模型來分析，看看房屋價格隨時間變化蔬充，在空間的分布規(guī)律。
首先班利，根據(jù)空間位置饥漫，計算Matern空間位置效應，產(chǎn)生一個mesh罗标，為降低運行時間
這里的Mesh設定值為10庸队，20积蜻。得到1020個三角形單元。
在繪制Mesh時候彻消，我們用到library(inlabru)包竿拆，可以將mesh結(jié)合到ggplot里面，inlabru也包含了模型結(jié)果的輸出宾尚，可自行探索丙笋。
# Location matrix
coords = df %>% dplyr:::select(Latitude,Longitude) %>%
as.matrix
mesh = inla.mesh.2d(coords, max.edge = c(10, 20))
summary(mesh)
# plot mesh
ggplot()+
gg(mesh)+
geom_point(data=df,aes(x=Latitude,y=Longitude))
2.2 spde
上述中關(guān)鍵一步與前幾期介紹的不一樣的地方是，我們加入了n.group參數(shù)煌贴，用來指定年份御板，先看一下Year有多少年的，然后封裝在s.index里面牛郑。
# spde
Year=length(table(df$year))
spde = inla.spde2.matern(mesh, alpha = 2)
s.index = inla.spde.make.index("w", n.spde = spde$n.spde, n.group = Year)
table(s.index$w.group)
A.price = inla.spde.make.A(mesh = mesh, loc =coords,
group = as.numeric(as.factor(df$year)),
n.group = Year)
2.3 stack
首先將我們的數(shù)據(jù)變成matrix怠肋，提取出Covariates，將year轉(zhuǎn)換成factor因子淹朋。在matrix后的變量笙各，會出現(xiàn)從2010-2014的變量，我們以2010為參照础芍，所以X=data.frame(Xm[,-2])杈抢，來去除2010年這一列。
# matrix
df=df %>% mutate(year=as.factor(year))
Xm <- model.matrix(~ -1 + area+year, data = df)
X=data.frame(Xm[,-2])
head(X)
N=nrow(df)
# satck
stack = inla.stack(tag="est",data=list(y=df$price),
A = list(1,1,A.price),
effects= list(
Intercept = rep(1, N),
X=X,
w = s.index))
dim(inla.stack.A(stack))
2.4 Model fit
建立好數(shù)據(jù)及mesh以后者甲，我們現(xiàn)在寫INLA公式春感，如果不考慮殘差項，那就是固定效應虏缸，現(xiàn)在我們增加殘差項鲫懒，及空間與時間效應。直接寫在放f（）里面刽辙。運行時間有點長窥岩，耐心等待。關(guān)于更多Spatial latent effects,可以參見inla.list.models宰缤。
這里我們使用AR1颂翼，時間自相關(guān)函數(shù)。
# inla
formula = as.formula(paste0("y ~ -1 + Intercept + ", paste0(colnames(X), collapse = " + "),
"+ f(w, model = spde,group = w.group,
control.group = list(model = 'ar1'))"))
formula
fit_inla = inla(formula,
data = inla.stack.data(stack), # Don't forget to change the stack!
control.compute = list(dic = TRUE),
control.predictor = list(A = inla.stack.A(stack)))
summary(fit_inla)
Call:
c("inla(formula = formula, data = inla.stack.data(stack), control.compute =
list(dic = TRUE), ", " control.predictor = list(A = inla.stack.A(stack)))")
Time used:
Pre = 3.05, Running = 162, Post = 2.21, Total = 167
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept -22.645 4.180 -30.852 -22.645 -14.444 -22.645 0
area 5.283 0.029 5.225 5.283 5.340 5.283 0
year2011 -122.201 3.617 -129.308 -122.200 -115.109 -122.196 0
year2012 -235.934 2.893 -241.615 -235.935 -230.257 -235.935 0
year2013 356.431 3.386 349.779 356.433 363.069 356.437 0
year2014 238.186 3.072 232.152 238.187 244.213 238.188 0
year2015 -0.376 3.285 -6.828 -0.375 6.066 -0.373 0
Random effects:
Name Model
w SPDE2 model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 0.001 0.000 0.001 0.001 0.001 0.001
Theta1 for w -3.997 0.049 -4.093 -3.998 -3.899 -4.000
Theta2 for w 0.973 0.067 0.857 0.967 1.117 0.947
GroupRho for w -0.665 0.021 -0.705 -0.665 -0.623 -0.665
Expected number of effective parameters(stdev): 56.81(4.64)
Number of equivalent replicates : 65.06
Deviance Information Criterion (DIC) ...............: 38409.53
Deviance Information Criterion (DIC, saturated) ....: 3937.44
Effective number of parameters .....................: 55.00
Marginal log-Likelihood: -19381.51
Posterior marginals for the linear predictor and
the fitted values are computed
上述結(jié)果慨灭，我們可以看到時間效應的先驗分布及area的估計值朦乏。
2.5 參數(shù)估計
從這個圖，可以看到在我們的INLA模型中氧骤，各個參數(shù)的先驗分布呻疹。主要是Range參數(shù)，可以提供空間相關(guān)性的距離筹陵。R-INLA-參數(shù)介紹
rf <- inla.spde2.result(fit_inla, "w", spde, do.transf = TRUE)
# str(rf)
par(mfrow = c(2, 2), mar = c(3, 3, 1, 0.1), mgp = 2:0)
plot(fit_inla$marginals.hyper[[2]], type = "l", xlab = "beta", ylab = "Density")
plot(rf$marginals.variance.nominal[[1]], type = "l", xlab = "sigma[x]", ylab = "Density")
plot(rf$marginals.kappa[[1]], type = "l", xlab = "kappa", ylab = "Density")
plot(rf$marginals.range.nominal[[1]], type = "l", xlab = "range nominal", ylab = "Density")
par(mfrow=c(1,1))
2.6 空間效應估計
因為隨時間變化刽锤，每一年的空間效應也不一樣镊尺，也就是INLA回歸方程中的殘差在空間上分布不均。我們可以將這種空間的效應畫出來并思。
同前述預測一樣庐氮，先用projector建立空的Grid，然后expand grid到研究區(qū)域宋彼。
然后提取模型中空間效應的的mean與sd弄砍。
# predict
Projection <- inla.mesh.projector(mesh,dims = c(300, 300))
Full.Projection <- expand.grid(x = Projection$x, y = Projection$y)
# get mean for each year
Model=fit_inla
Full.Projection[,paste0("Year",1:Year)] =
apply(matrix(Model$summary.random$w$mean, ncol = Year), 2,
function(x) c(inla.mesh.project(Projection, x)))
head(Full.Projection)
x y Year1 Year2 Year3 Year4 Year5 Year6
1 -28.29342 -55.87092 NA NA NA NA NA NA
2 -28.26937 -55.87092 NA NA NA NA NA NA
3 -28.24533 -55.87092 NA NA NA NA NA NA
4 -28.22129 -55.87092 NA NA NA NA NA NA
5 -28.19725 -55.87092 NA NA NA NA NA NA
6 -28.17321 -55.87092 NA NA NA NA NA NA
# to raster
library(raster)
dfra=c()
for(i in 1:Year){
raster_df=Full.Projection %>% dplyr::select(x,y,paste0("Year",i))
dfr = rasterFromXYZ(raster_df) #Convert first two columns as lon-lat and third as value
#year1=mask(dfr,studyarea)
dfra=raster::stack(dfr,dfra)
}
plot(dfra)
上述的空間分布預測只顯示了空間效應，如何添加Covariate及year宙暇，參見INLA prediction貝葉斯地理統(tǒng)計模型R-INLA-3输枯。
參考
1.Geostatistical data
2.Spatial analysis of geotagged data
3.Spatial and spatio-temporal models with INLA
4.Spatio-temporal modeling of areal data. Lung cancer in Ohio
5. Spatio-temporal modeling with INLA
