基于GAM的时空疾病风险分析：年度差异与预测可视化问题

问题概述 #

使用广义可加模型（GAM）分析疾病时空风险时，初始模型无法有效捕捉2012-2014年的年度差异；将Year转为因子后模型拟合正常，但需要生成各年度的预测结果空间可视化图。

调整后的模型代码：

mod2 = gam(TB ~ offset(log(Population)) + 
             s(Indigenous, k = 10, bs = "cr") + 
             s(Urbanisation, k = 10, bs = "cr") +
             s(Density, k = 10, bs = "cr") + 
             Poverty + 
             s(Poor_Sanitation, k = 10, bs = "cr") + 
             Unemployment + 
             Timeliness +
             as.factor(Year) + 
             Region + 
             s(lon, lat), 
           data = TBdata, 
           family = nb(link = 'log'))

生成年度预测可视化步骤 #

1. 构建预测数据集 #

需为2012/2013/2014三个目标年份生成包含所有模型协变量的数据集，常见两种构建方式：

• 固定协变量均值：用于展示排除其他协变量干扰后的纯年度空间风险差异
• 保留原始协变量变异：用于展示每个空间单元在不同年份的实际风险变化

以下是固定协变量均值的实现代码：

library(mgcv)
library(dplyr)
library(tidyr)

# 提取原始数据中的空间核心信息（经纬度、区域）
spatial_units <- distinct(TBdata, lon, lat, Region)

# 生成跨年度预测数据集，非空间协变量取全局均值
pred_data <- expand_grid(
  spatial_units,
  Year = c(2012, 2013, 2014),
  Indigenous = mean(TBdata$Indigenous, na.rm = TRUE),
  Urbanisation = mean(TBdata$Urbanisation, na.rm = TRUE),
  Density = mean(TBdata$Density, na.rm = TRUE),
  Poverty = mean(TBdata$Poverty, na.rm = TRUE),
  Poor_Sanitation = mean(TBdata$Poor_Sanitation, na.rm = TRUE),
  Unemployment = mean(TBdata$Unemployment, na.rm = TRUE),
  Timeliness = mean(TBdata$Timeliness, na.rm = TRUE),
  Population = mean(TBdata$Population, na.rm = TRUE)  # 建议替换为各空间单元实际人口
)

# 转换Year为因子，与模型拟合格式一致
pred_data$Year <- as.factor(pred_data$Year)

2. 获取模型预测值 #

使用predict()函数生成响应尺度的预测值（即实际发病风险，而非log尺度）：

# 生成响应尺度的预测值
pred_data$TB_pred <- predict(mod2, newdata = pred_data, type = "response")

3. 绘制年度空间预测图 #

点数据可视化（ggplot2）

若数据为点位观测（如区县中心点）：

library(ggplot2)

ggplot(pred_data, aes(x = lon, y = lat, color = TB_pred)) +
  geom_point(size = 3, alpha = 0.8) +
  facet_wrap(~Year, ncol = 2) +
  scale_color_viridis_c(option = "plasma", name = "预测TB风险") +
  theme_minimal() +
  labs(title = "2012-2014年TB时空风险预测", x = "经度", y = "纬度")

面数据可视化（tmap）

若数据为行政区面数据（如Shapefile），先合并空间对象与预测数据：

library(tmap)
library(sf)

# 假设spatial_sf是你的空间面对象
pred_sf <- left_join(spatial_sf, pred_data, by = c("lon", "lat", "Region"))

tm_shape(pred_sf) +
  tm_fill("TB_pred", palette = "plasma", title = "预测TB风险", alpha = 0.7) +
  tm_borders(col = "white", lwd = 0.5) +
  tm_facets(by = "Year") +
  tm_layout(main.title = "2012-2014年TB时空风险预测", legend.outside = TRUE)

注意事项 #

• Population建议使用各空间单元的实际人口数，而非全局均值，避免offset偏差影响预测结果。
• 确保预测数据的lon/lat范围与原始数据一致，超出范围的空间平滑预测可靠性会下降。
• 若需展示预测置信区间，可在predict()中添加se.fit = TRUE，并在可视化时绘制误差范围。

内容的提问来源于stack exchange，提问作者Joe