贝叶斯地理统计模型R-INLA-4 贝叶斯时空模型
在前述的内容中，我们介绍了，如何处理空间的数据，利用海拔高度预测降雨量的例子。但是该例子仅仅涉及到的是涉及到回归方程中，考虑影响因素及空间效应。
那么如果我们的数据有时间信息，如何加入到贝叶斯时空分析呢。譬如每年对某一个地区进行疾病的发病率调查，10年数据整合在一起，就可以从时间上或空间上看疾病的变化规律，也就会用到贝叶斯时空模型。
下面我们将介绍贝叶斯时空模型。该文章中，会简化数学计算的过程，主要是针对，在有数据的基础上，如何应用贝叶斯时空模型，找出影响因素，绘制时间变化的空间分布预测图。
1.数据集
我们模拟一个房屋与面积的地理位置数据，从2010年到2015年的房地产数据。因为是模拟数据，不具有任何实际意义，仅仅作为展示。
通过简单的回归方程，发现，房屋价格与面积及年份成正相关，具有统计学意义。说明随时间推迟，房子越值钱，且面积越大价格也越高。符合实际
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 的步骤包括产生Mesh，建立A project，将空间位置与三角形网格点对应，计算spde的空间效应。然后整合data到stack里面。这是建立INLA的关键，最后，写INLA的公式，带入INLA模型。
2.1 Mesh
下面我们利用时空模型来分析，看看房屋价格随时间变化，在空间的分布规律。
首先，根据空间位置，计算Matern空间位置效应，产生一个mesh，为降低运行时间
这里的Mesh设定值为10，20。得到1020个三角形单元。
在绘制Mesh时候，我们用到library(inlabru)包，可以将mesh结合到ggplot里面，inlabru也包含了模型结果的输出，可自行探索。
# 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
上述中关键一步与前几期介绍的不一样的地方是，我们加入了n.group参数，用来指定年份，先看一下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
首先将我们的数据变成matrix，提取出Covariates，将year转换成factor因子。在matrix后的变量，会出现从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
建立好数据及mesh以后，我们现在写INLA公式，如果不考虑残差项，那就是固定效应，现在我们增加残差项，及空间与时间效应。直接写在放f（）里面。运行时间有点长，耐心等待。关于更多Spatial latent effects,可以参见inla.list.models。
这里我们使用AR1，时间自相关函数。
# 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
上述结果，我们可以看到时间效应的先验分布及area的估计值。
2.5 参数估计
从这个图，可以看到在我们的INLA模型中，各个参数的先验分布。主要是Range参数，可以提供空间相关性的距离。R-INLA-参数介绍
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到研究区域。
然后提取模型中空间效应的的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贝叶斯地理统计模型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
