使用gratia::derivatives计算sz型因子平滑导数时遇错误求助

问题分析 #

使用gratia::derivatives处理带bs="sz"的约束因子平滑项s(Year,Site)时，出现Error in quos(): ! The LHS of := must be a string, not a character vector报错，大概率是旧版gratia对这类带因子的平滑项解析存在bug，或是输入数据的格式不符合函数要求。

解决方案 #

方案1：更新gratia到最新版本 #

旧版本（<0.8.2）的gratia对bs="sz"类型的平滑项导数计算支持不完善，先执行以下命令更新包：

install.packages("gratia")

更新后，用gratia::data_slice生成规范的输入数据（包含所有Site水平和均匀分布的Year序列），再调用derivatives：

# 生成包含所有站点和均匀年份序列的数据集
new_data <- gratia::data_slice(mod, Year = gratia::evenly(Year, n = 100), Site = unique(Site))

# 计算一阶导数及同时置信区间
deriv_site <- derivatives(
  mod,
  term = "s(Year,Site)",
  data = new_data,
  order = 1,
  type = "central",
  interval = "simultaneous",
  n_sim = 10000,
  level = 0.95,
  unconditional = TRUE
)

# 绘图可视化
plot(deriv_site, colour = "Site") +
  ggplot2::geom_hline(yintercept = 0, linetype = "dashed", colour = "black")

方案2：手动数值微分+模拟置信区间 #

如果更新版本后仍报错，可通过手动数值微分结合参数模拟的方式实现需求：

library(dplyr)
library(ggplot2)
library(MASS)

# 生成输入数据
new_data <- gratia::data_slice(mod, Year = gratia::evenly(Year, n = 100), Site = unique(Site))

# 提取s(Year,Site)对应的系数和设计矩阵列
smooth_term <- "s(Year,Site)"
coef_idx <- grep(smooth_term, names(coef(mod)))
pred_mat <- predict(mod, newdata = new_data, type = "lpmatrix")[, coef_idx]
coef_smooth <- coef(mod)[coef_idx]
vcov_smooth <- vcov(mod)[coef_idx, coef_idx]

# 中心差分计算一阶导数
h <- diff(range(new_data$Year)) / 1000  # 设定微分步长
new_data_plus <- new_data %>% mutate(Year = Year + h)
new_data_minus <- new_data %>% mutate(Year = Year - h)

pred_plus <- predict(mod, newdata = new_data_plus, type = "lpmatrix")[, coef_idx]
pred_minus <- predict(mod, newdata = new_data_minus, type = "lpmatrix")[, coef_idx]

deriv_vals <- as.vector((pred_plus %*% coef_smooth - pred_minus %*% coef_smooth) / (2 * h))

# 模拟同时置信区间
set.seed(123)
sim_coef <- mvrnorm(n = 10000, mu = coef_smooth, Sigma = vcov_smooth)
sim_deriv <- (pred_plus %*% t(sim_coef) - pred_minus %*% t(sim_coef)) / (2 * h)

deriv_lower <- apply(sim_deriv, 1, quantile, probs = 0.025)
deriv_upper <- apply(sim_deriv, 1, quantile, probs = 0.975)

# 整理结果
deriv_result <- new_data %>%
  mutate(
    derivative = deriv_vals,
    lower = deriv_lower,
    upper = deriv_upper
  )

# 可视化
ggplot(deriv_result, aes(x = Year, y = derivative, color = Site)) +
  geom_line(linewidth = 0.8) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = Site), alpha = 0.2, color = NA) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  labs(x = "年份", y = "一阶导数", color = "站点", fill = "站点") +
  theme_minimal()

关键注意事项 #

• 确保输入数据包含所有需要分析的Site水平，否则会遗漏部分站点的导数结果
• 使用unconditional = TRUE时，需确保模型计算了无条件协方差矩阵（gam模型默认已支持）
• 模拟同时置信区间时，设置随机种子保证结果可复现

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