如何绘制GAM模型预测结果以分析2012-2014年时空趋势？

首先，你的当前模型将Year作为主效应纳入，这意味着它假设空间模式（来自s(lon, lat)）在各年份间是一致的，年份仅带来TB计数的全局偏移。如果想要捕捉随时间变化的空间模式（真正的时空趋势），需要调整模型，让空间平滑项与年份交互：

步骤1：调整模型（可选，但推荐用于时空趋势） #

# 包含年份与空间平滑的交互项
mod_spatiotemp <- 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 +
    Region + 
    s(lon, lat, by = as.factor(Year)),  # 按年份拆分空间平滑
  data = TBdata, 
  family = nb(link = 'log')
)

步骤2：构建预测数据集 #

为了生成各年份的空间预测，需要创建覆盖研究区域经纬度的网格，并为每个网格点匹配所有年份。同时将其他协变量设为均值（连续变量）或参考水平（分类变量），以隔离时空效应：

library(tidyverse)
library(mgcv)

# 获取经纬度范围
lon_range <- range(TBdata$lon, na.rm = TRUE)
lat_range <- range(TBdata$lat, na.rm = TRUE)

# 创建经纬度网格（100x100分辨率，可调整）
lon_grid <- seq(lon_range[1], lon_range[2], length.out = 100)
lat_grid <- seq(lat_range[1], lat_range[2], length.out = 100)

# 组合网格与年份
pred_data <- expand.grid(
  lon = lon_grid,
  lat = lat_grid,
  Year = unique(TBdata$Year)
) %>%
  # 设置其他协变量为均值/参考值
  mutate(
    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),
    Region = first(TBdata$Region),  # 使用数据中第一个区域（或替换为众数）
    Population = mean(TBdata$Population, na.rm = TRUE),  # 固定人口为均值，确保跨年份可比
    Year = as.factor(Year)  # 匹配模型中的因子类型
  )

步骤3：生成模型预测 #

使用调整后的模型（或原模型）计算每个网格点的预测值：

# 用调整后的时空模型预测
pred_data$pred_TB <- predict(mod_spatiotemp, newdata = pred_data, type = "response")

# 如果用原模型，替换为：
# pred_data$pred_TB <- predict(mod, newdata = pred_data, type = "response")

步骤4：用ggplot绘制时空趋势 #

使用分面（facet）按年份展示空间预测结果：

ggplot(pred_data, aes(x = lon, y = lat, fill = pred_TB)) +
  geom_raster(interpolate = TRUE) +  # 插值让图像更平滑
  facet_wrap(~Year, ncol = 3) +  # 按年份分面，3列布局
  scale_fill_viridis_c(option = "plasma", name = "预测TB计数") +  # 配色方案
  coord_fixed() +  # 保持经纬度比例一致
  labs(
    x = "经度",
    y = "纬度",
    title = "2012-2014年TB预测计数的时空趋势"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.position = "bottom"
  )

补充说明 #

• 如果使用原模型（无空间-年份交互），各年份的空间图案会完全一致，仅颜色深浅（代表TB计数）随年份变化，这反映的是全局年份效应而非空间模式的变化。
• 若想展示发病率（而非计数），可以将预测值除以Population，并修改填充标签为“预测TB发病率”。

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