使用INLAspacetime完成时空建模
library(sf)
library(fdmr)
library(fmesher)
library(INLA)
library(inlabru)
library(INLAspacetime)
library(tidyverse)
# 1. 加载数据并转换为 sf 对象
sf_data <- readRDS('data/spatial_dataBris.rds') %>%
st_as_sf(coords = c("LONG", "LAT"), crs = 4326) %>%
mutate(mapp = 0)
covid19_sf <- readRDS('data/covid19_dataBris.rds') %>%
st_as_sf(coords = c("LONG", "LAT"), crs = 4326)
# 2. 计算网格参数,sf_data 的 geometry 类型是 sfc_MULTIPOLYGON
coords <- st_coordinates(sf_data)
initial_range <- diff(range(coords[, "X"])) / 3
max_edge <- initial_range / 2
# 3.mesh
smesh <- fm_mesh_2d_inla(sf_data,
max.edge=c(1, 2)* max_edge,
offset = c(initial_range, initial_range),
cutoff=max_edge/7)
plot(smesh)
nt <- max(covid19_sf$week)
h <- 1
tmesh <- fm_mesh_1d(
loc = seq(1, nt + h/2, h),
degree = 1)
tmesh$n
stmodel <- stModel.define(
smesh = smesh, ## spatial mesh
tmesh = tmesh, ## temporal mesh
model = '202', ## model, see the paper 102,121,202,220
control.priors = list(
prs = c(1, 0.1), ## P(spatial range < 1) = 0.1
prt = c(5, 0), ## temporal range fixed to 5
psigma = c(1, 0.1) ## P(sigma > 1) = 0.1
)
)
components <- cases ~ 1 + prevalence + field(
list(space = st_coordinates(covid19_sf), time = week), model = stmodel)
result <- bru(components,
family = "gaussian",
data = covid19_sf,
options = list(
control.inla = list(
# 使用metis重排序以提高计算效率
reordering = "metis",
# 使用经验贝叶斯积分策略
int.strategy = "eb"
),
verbose = FALSE,
inla.mode = "experimental",
control.compute = list(waic = TRUE,dic = TRUE,cpo = TRUE)
)
)
summary(result)
评论
发表评论