如何用gratia::fitted_values()获取HGAM中全局平滑项的响应预测

假设你的HGAM模型结构类似（基于mgcv构建）：

library(mgcv)
# 示例模型：全局平滑s(OrdDay) + 湖泊分组平滑s(OrdDay, by = Lake)
model <- gam(D ~ s(OrdDay) + s(OrdDay, by = Lake), data = your_dataset, method = "REML")

方法一：直接提取全局平滑项的估计值 #

用gratia::smooth_estimates()可以精准提取指定平滑项的拟合结果，避开分组项的干扰：

library(gratia)

# 查看模型中所有平滑项名称，确认全局项的标识
smooths(model)

# 提取全局平滑s(OrdDay)的拟合值
global_smooth_fit <- smooth_estimates(model, smooth = "s(OrdDay)")

# 结果中的`est`列即为各OrdDay对应的平均预测值
head(global_smooth_fit)

方法二：通过predict()排除分组项效应 #

构造仅包含OrdDay的新数据，用exclude参数屏蔽分组平滑项，得到全局平均预测：

# 生成覆盖所有OrdDay取值的新数据框
new_data <- data.frame(OrdDay = unique(your_dataset$OrdDay))

# 排除分组平滑项，只计算全局项+截距的响应值
global_fitted <- predict(model, newdata = new_data, 
                         exclude = "s(OrdDay):Lake", # 需与模型中分组项名称一致
                         type = "response")

# 组合成带OrdDay的结果数据框
global_preds <- cbind(new_data, avg_fitted = global_fitted)

为什么之前fitted_values()只返回单个湖泊结果？ #

gratia::fitted_values()默认会保留所有模型项的效应，当你未指定Lake变量时，它会自动选取数据集中第一个出现的湖泊值代入计算，因此得到的是单个湖泊的预测而非全局平均。上述两种方法分别从「提取平滑项本身」和「屏蔽分组效应」两个维度解决了这个问题。

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