Prediction
上一期我們介紹了如何來評估INLA模型吼驶，因?yàn)榭臻g位置點(diǎn)的預(yù)測需要驗(yàn)證一铅。那整個流程走完以后疾瓮，最后一步就是對其他地區(qū)進(jìn)行預(yù)測。
我們可以通過計(jì)算新位置的投影矩陣，然后將投影矩陣乘以空間場值，將這些值投影到不同的位置。例如，我們可以如下計(jì)算矩陣newloc中位置處的空間場的后均值：
A test
# projector
newloc <- cbind(c(-90, -78, 18), c(20, 20, 10))
Aproj <- inla.spde.make.A(Mesh, loc = newloc)
a=Aproj %*% fit$summary.random$w$mean
a
[,1] [,2]
[1,] -90 20
[2,] -78 20
[3,] 18 10
研究區(qū)域的空間grid
使用inla.mesh.projector（）和inla.mesh.project（）函數(shù)在不同位置投影空間場值纷捞。首先，我們需要使用inla.mesh.projector（）函數(shù)為新位置計(jì)算投影矩陣键俱。我們可以在參數(shù)loc中指定位置兰绣，也可以通過指定參數(shù)xlim，ylim和dims在網(wǎng)格上計(jì)算位置编振。例如缀辩，我們使用inla.mesh.projector（）計(jì)算覆蓋網(wǎng)格區(qū)域的網(wǎng)格上500 x 500個位置的投影矩陣。
然后踪央，我們使用inla.mesh.project（）函數(shù)將在網(wǎng)格節(jié)點(diǎn)處計(jì)算上述空間網(wǎng)格點(diǎn)上對應(yīng)的均值與標(biāo)準(zhǔn)差臀玄。然后使用expand.grid（）函數(shù)組合proj $ x和proj $ y來獲得網(wǎng)格位置的坐標(biāo)。
rang <- apply(Mesh$loc[, c(1, 2)], 2, range)
proj <- inla.mesh.projector(Mesh,
xlim = rang[, 1], ylim = rang[, 2],
dims = c(500, 500))
mean_s <- inla.mesh.project(proj, fit$summary.random$w$mean)
sd_s <- inla.mesh.project(proj, fit$summary.random$w$mean)
df <- expand.grid(x = proj$x, y = proj$y)
df$mean_s <- as.vector(mean_s)
繪圖
這里只展示如何繪制mean畅蹂，并且將grid數(shù)據(jù)框轉(zhuǎn)換成raster格式健无，就可以利用圖層進(jìn)行crop與mask
ggplot(df, aes(x = x, y = y, fill = mean_s)) +
geom_raster() +
coord_fixed(ratio = 1) + theme_bw()
# to raster
library(raster)
dfr <- rasterFromXYZ(df) #Convert first two columns as lon-lat and third as value
plot(dfr)
這里展示的僅僅是降雨量數(shù)據(jù)預(yù)測，后面還需要涉及到回歸參數(shù)的納入液斜。
當(dāng)然預(yù)測出了系數(shù)累贤，需要考慮影響因素效應(yīng)叠穆，未完待續(xù)。
參考
1.Geostatistical data
2.Spatial analysis of geotagged data
3.Spatial and spatio-temporal models with INLA
