如何跨时间处理虚拟变量？基于GAM模拟冷暖冬季预测昆虫种群

Hey Heather, 我来帮你搞定这个虚拟数据集构建的问题！用GAM做昆虫种群预测，要模拟暖冬/寒冬场景，核心就是精准调整冬季时段的温度，同时保持其他季节的温度模式和原有数据一致，这样才能单独看冬季温度对种群的影响。下面给你一步步拆解，还有具体的R代码示例：

首先得从你的现有数据里提取季节/月份的温度基准，这样虚拟数据的非冬季部分才符合真实情况：

library(lubridate)
library(mgcv)

# 假设你的原始数据集叫dat，包含date（日期）和temperature（温度）列
# 给数据加月份和季节标签（冬季定义可以根据你的研究区域调整，比如北半球是12-2月）
dat$month <- month(dat$date)
dat$season <- case_when(
  dat$month %in% c(12, 1, 2) ~ "winter",
  dat$month %in% c(3, 4, 5) ~ "spring",
  dat$month %in% c(6, 7, 8) ~ "summer",
  TRUE ~ "autumn"
)

# 计算每个月份的平均温度（用来填充虚拟数据的非冬季月份）
monthly_temp_avg <- tapply(dat$temperature, dat$month, mean, na.rm = TRUE)

# 统计原有冬季温度的关键值：用来定义寒冬/暖冬
winter_temps <- dat$temperature[dat$season == "winter"]
winter_avg <- mean(winter_temps, na.rm = TRUE)
# 寒冬可以取原有冬季温度的5%分位数（也就是比95%的冬季都冷）
winter_cold_threshold <- quantile(winter_temps, 0.05, na.rm = TRUE)
# 暖冬取95%分位数（比95%的冬季都暖）
winter_warm_threshold <- quantile(winter_temps, 0.95, na.rm = TRUE)

先生成你要预测的未来1年的时间轴，比如按月度数据（如果是按天/周，把by="month"改成by="day"或by="week"就行）：

# 生成未来12个月的日期序列（起始日期可以根据你的需求改）
future_dates <- seq.Date(from = as.Date("2024-01-01"), length.out = 12, by = "month")

# 构建基础虚拟数据框，包含日期、月份、季节
new_dat_base <- data.frame(date = future_dates)
new_dat_base$month <- month(new_dat_base$date)
new_dat_base$season <- case_when(
  new_dat_base$month %in% c(12, 1, 2) ~ "winter",
  new_dat_base$month %in% c(3, 4, 5) ~ "spring",
  new_dat_base$month %in% c(6, 7, 8) ~ "summer",
  TRUE ~ "autumn"
)

核心逻辑：非冬季月份用原有数据的月均温，冬季月份替换成寒冬/暖冬的阈值：

# --- 寒冬场景数据集 ---
new_dat_cold <- new_dat_base
new_dat_cold$temperature <- ifelse(
  new_dat_cold$season == "winter",
  winter_cold_threshold,  # 冬季用寒冬阈值
  monthly_temp_avg[as.character(new_dat_cold$month)]  # 非冬季用对应月份的平均温度
)

# --- 暖冬场景数据集 ---
new_dat_warm <- new_dat_base
new_dat_warm$temperature <- ifelse(
  new_dat_warm$season == "winter",
  winter_warm_threshold,  # 冬季用暖冬阈值
  monthly_temp_avg[as.character(new_dat_warm$month)]
)

假设你的GAM模型已经构建好了（比如包含温度的平滑项、时间/季节的平滑项），直接用predict()函数就行：

# 假设你的模型是gam_model，比如：
# gam_model <- gam(population ~ s(temperature) + s(month, bs = "cc"), data = dat)

# 预测寒冬场景的种群数量
new_dat_cold$predicted_pop <- predict(gam_model, newdata = new_dat_cold, type = "response")

# 预测暖冬场景的种群数量
new_dat_warm$predicted_pop <- predict(gam_model, newdata = new_dat_warm, type = "response")

• 自定义寒冬/暖冬定义：如果分位数不符合你的研究需求，也可以直接给冬季温度减/加固定值（比如winter_avg - 3表示比平均冬季冷3℃），只要符合你对“寒冬”的定义就行。
• 其他协变量的处理：如果你的GAM模型里还有其他协变量（比如降水、植被覆盖），要确保这些变量在寒冬/暖冬场景中保持一致（比如用历史均值填充），这样才能单独隔离温度的影响。
• 外推风险：如果你的寒冬/暖冬温度超出了原始数据的温度范围，GAM的平滑项外推可能会有不确定性，建议检查模型的拟合范围，或者在模型中用xt=list(boundary=c(min_temp, max_temp))限制平滑项的外推范围。

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