如何在R中提取GAM（张量积平滑）模型特定周的绘图数据？

我构建了一个包含位置与时间张量积平滑的GAM模型，代码如下：

dat <- gratia::bird_move
dat$longitude <- dat$latitude

mod <- mgcv::gam(count ~ species + te(longitude, latitude, week, d = c(2,1), k = c(8,8)), data = dat)

我知道可以通过plot(mod)获取绘图数据，但如何提取并绘制例如第48周的数据？

运行plot(mod)后得到的绘图数据结构如下：

plotdata <- plot(mod) 
> str(plotdata)
List of 1
 $ :List of 12
  ..$ scale  : logi FALSE
  ..$ se     : logi FALSE
  ..$ m      : num 40
  ..$ nc     : num 3
  ..$ nr     : num 3
  ..$ lo     : num [1:3] 5 5 0
  ..$ hi     : num [1:3] 60 60 52
  ..$ vname  : chr [1:3] "longitude" "latitude" "week"
  ..$ main   : chr "te(longitude,latitude,week,37.46)"
  ..$ exclude: logi [1:14400] FALSE FALSE FALSE FALSE TRUE TRUE ...
  ..$ fit    : num [1:14400, 1] 4.98 4.9 4.83 4.74 NA ...
  ..$ plot.me: logi TRUE

方法一：从plot(mod)返回的数据中提取 #

plot(mod)返回的列表第一个元素就是张量积项的绘图数据，可按以下步骤筛选第48周的数据：

1. 提取核心绘图数据对象

pd <- plotdata[[1]]

2. 生成各变量的网格序列
   根据pd中的范围和点数，生成经度、纬度、周数对应的网格值：

# 经度和纬度各有40个网格点（pd$m=40）
lon_seq <- seq(pd$lo[1], pd$hi[1], length.out = pd$m)
lat_seq <- seq(pd$lo[2], pd$hi[2], length.out = pd$m)
# 计算周数的网格点数：总拟合值数 / 空间网格数（40*40）
week_count <- length(pd$fit) / (pd$m * pd$m)
week_seq <- seq(pd$lo[3], pd$hi[3], length.out = week_count)

3. 定位第48周对应的网格索引
   找到周数序列中最接近48的位置：

target_idx <- which.min(abs(week_seq - 48))

4. 提取该周的拟合值并转为矩阵
   拟合值按「所有空间网格组合 + 周数顺序」排列，提取对应片段并转为空间矩阵：

fit_matrix <- matrix(
  pd$fit[(target_idx-1)*pd$m*pd$m + 1:(pd$m*pd$m)],
  nrow = pd$m, ncol = pd$m
)

5. 绘制空间热图
   用ggplot2可视化：

library(ggplot2)
# 转换为数据框
plot_df <- expand.grid(longitude = lon_seq, latitude = lat_seq)
plot_df$fit <- as.vector(fit_matrix)

ggplot(plot_df, aes(x = longitude, y = latitude, fill = fit)) +
  geom_raster() +
  scale_fill_viridis_c() +
  labs(title = paste("第", round(week_seq[target_idx]), "周拟合值"),
       x = "经度", y = "纬度", fill = "拟合计数") +
  theme_minimal()

方法二：直接用predict生成精确数据（更推荐） #

如果需要精确的第48周拟合值，建议直接通过模型预测生成数据，避免依赖plot的默认网格：

1. 构建预测网格
   生成第48周的空间网格，同时指定物种（模型包含species协变量）：

pred_grid <- expand.grid(
  longitude = seq(min(dat$longitude), max(dat$longitude), length.out = 100),
  latitude = seq(min(dat$latitude), max(dat$latitude), length.out = 100),
  week = 48,
  species = unique(dat$species)[1] # 可替换为其他物种或循环所有物种
)

2. 预测拟合值

pred_grid$fit <- predict(mod, newdata = pred_grid)

3. 可视化

ggplot(pred_grid, aes(x = longitude, y = latitude, fill = fit)) +
  geom_raster() +
  scale_fill_viridis_c() +
  labs(title = "第48周拟合值", x = "经度", y = "纬度", fill = "拟合计数") +
  theme_minimal()

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