拟合时空耦合模型

library(sf)

library(fdmr)

library(fmesher)

library(INLA)

library(inlabru)

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

mesh <- fm_mesh_2d_inla(sf_data,

  max.edge=c(1, 2)* max_edge,

  offset = c(initial_range, initial_range),

  cutoff=max_edge/7)

plot(mesh)

# 5. 定义 SPDE 模型

prior_range <- initial_range

spde <- INLA::inla.spde2.pcmatern(

  mesh = mesh,

  prior.range = c(prior_range, 0.5),

  prior.sigma = c(1, 0.01)

)

# 6.定义时间索引和离散时间点数量

group_index <- covid19_sf$week

n_groups <- length(unique(covid19_sf$week))

rhoprior <- base::list(theta = list(prior = "pccor1", param = c(0, 0.9)))

# 7.定义模型公式

formula <- cases ~ 0 + Intercept(1) + prevalence +

  f(

    main = st_coordinates,

    model = spde,

    group = group_index,

    ngroup = n_groups,

    control.group = list(

      model = "ar1",

      hyper = rhoprior

    )

  )

# 8.拟合时空耦合模型

inlabru_model <- inlabru::bru(

  formula,

  data = covid19_sf,

  family = "poisson",

  # Population作为泊松模型的期望值

  E = covid19_sf$Population,

  control.family = list(link = "log"),

  # INLA配置选项

  options = list(

    control.inla = list(

      # 使用metis重排序以提高计算效率

      reordering = "metis",

      # 使用经验贝叶斯积分策略

      int.strategy = "eb"

    ),

    # 显示详细输出

    verbose = FALSE,

    # 使用INLA的实验模式

    inla.mode = "experimental"

  )

)

# 模型摘要

summary(inlabru_model)