基于R-INLA的时空疾病映射建模：空间图文件构建求助

1. 获取并预处理巴西城市多边形数据 #

先通过geobr获取目标城市的多边形数据，重点保证和你的疾病数据的城市编码（巴西通用IBGE编码）一致：

library(geobr)
library(spdep)
library(INLA)
library(dplyr)

# 获取目标区域的城市多边形（示例为全巴西2020年数据，可替换为特定州编码或年份）
cities_shp <- read_municipality(code_muni = "all", year = 2020)

# 统一编码格式，匹配你的疾病数据字段（假设你的数据用`city_code`存储编码）
cities_shp <- cities_shp %>% 
  mutate(city_code = as.character(code_muni))

2. 生成空间邻接关系（nb对象） #

用poly2nb基于多边形创建邻接列表，可选择"queen"规则（共享顶点也算邻接，适合城市尺度）或"rook"规则（仅共享边界）：

# 创建邻接列表
city_nb <- poly2nb(cities_shp, queen = TRUE)

# 可选：可视化验证邻接关系
plot(cities_shp$geometry)
plot(city_nb, coordinates(cities_shp), add = TRUE, col = "darkred")

3. 转换为INLA兼容的graph文件 #

用nb2INLA将邻接列表转成INLA建模所需的graph格式，文件会自动保存到工作目录：

# 生成INLA可用的graph文件
inla_graph <- nb2INLA("city_spatial_graph", city_nb)

4. 整合时空数据与空间索引 #

你的疾病数据是城市-年-月三元组，需对齐空间邻接图的节点顺序：

# 加载你的疾病数据（示例为CSV格式）
disease_data <- read.csv("your_disease_data.csv")

# 过滤掉不在多边形数据中的城市，确保编码完全匹配
disease_data <- disease_data %>% 
  filter(city_code %in% cities_shp$city_code)

# 为每个城市分配空间节点ID（与graph文件的节点顺序一致）
disease_data <- disease_data %>% 
  mutate(spatial_id = match(city_code, cities_shp$city_code))

# 生成连续时间索引（将年-月转为从起始时间开始的连续月份数）
disease_data <- disease_data %>% 
  mutate(date = as.Date(paste(year, month, "01", sep = "-")),
         time_id = as.integer(date - min(date)) %/% 30 + 1)

5. 构建INLA时空模型 #

将graph参数传入空间随机效应，同时结合时间效应与时空交互效应：

# 泊松分布模型公式（假设病例数为计数数据）
model_formula <- case_count ~ 1 + 
  f(spatial_id, model = "besag", graph = "city_spatial_graph") +  # 空间自相关效应
  f(time_id, model = "rw2") +  # 时间平滑效应
  f(spatial_id:time_id, model = "iid")  # 时空交互效应

# 拟合模型并开启预测计算
inla_model <- inla(model_formula, 
                   data = disease_data,
                   family = "poisson",
                   control.predictor = list(compute = TRUE))

# 查看模型结果与预测值
summary(inla_model)
predicted_cases <- inla_model$summary.fitted.values

关键注意事项 #

• 编码一致性：必须保证geobr返回的IBGE编码和你的疾病数据编码完全匹配，否则会出现邻接关系错位。
• 邻接规则选择：城市尺度建议用queen=TRUE，大区域分析可改用queen=FALSE。
• 时空模型调整：除示例的组合效应，还可尝试model="stcar"等专门的时空自相关模型，根据数据分布调整。

内容的提问来源于stack exchange，提问作者Moisés Augusto