GAM 中部分图的反向变换轴

Question

我创建了一个包含 3 个预测变量的 GAM。由于我构建模型所使用的一些预测变量和响应数据具有非正态分布，因此我对它们进行了对数变换。此外，因为我的响应变量有很多 0，所以我必须向每个数据点添加一个非常小的数字，以便日志对其进行转换。现在我已经完成了模型，但为了解释起见，我认为如果我可以对两个轴进行反向变换，其他人会更容易理解部分效应图；但是，我似乎无法弄清楚如何在部分效果上生成反向变换轴。

这是默认图之一的示例。

BestModel <- gam(Density ~ s(x) + s(y) + s(z), data = df, sp = c(10000, 0.01,10,10000))
plot(BestModel)

这里的两个轴都是对数转换的。我在原始单位中最接近这条曲线的方法是生成一个新的数据集，其中一个预测变量增加，而其他预测变量保持不变。我根据模型的反向变换预测值与反向变换的预测变量进行了绘制。

df2 <- NewDataset
plot(((exp(df2$y)*900)/10000), exp(predict(BestModel, newdata = df2))-0.006993007, 
xlab = "Back Transformed y", ylab = "Back Transformed Density")

这是我期望的曲线，但现在我没有平滑的线或置信区间，而且它看起来不如默认图（BestModel）生成的部分图那么好。

我尝试使用 ggplot2 添加置信区间和更好的线条。为了在这篇文章中进行比较，我在图中留下了点和穿过它们的线。

ggplot(data = df, aes(x = (exp(y)*900)/10000,
    y = exp(predict(BestModel, newdata = df))-0.006993007)) + 
      geom_point()+
      geom_line()+
      labs(x = "y", y = "Density")+
      stat_smooth(method = "gam", se = TRUE)

不幸的是，这似乎并没有产生正确的曲线，切换到“黄土”方法也没有产生正确的曲线。

有没有更简单的方法来单独反向变换每个轴以匹配其原始单位并获得每个预测变量的部分图？

Answer 1

有两种相对简单的方法可以做到这一点

1. 评估协变量已知值的平滑度，将变换应用于这些评估值，然后根据协变量的原始值绘制它们，或者
2. 根据模型进行预测，排除其他模型项，然后根据协变量的原始值绘制预测值

考虑 Simon Wood 的广义加法模型中的这个例子； R 的介绍

#install.packages("gamair")
#remotes::install_github("gavinsimpson/gratia")
library("mgcv")
library("gratia")
library("ggplot2")
library("dplyr")

data(mack, package = "gamair")
mack <- mack |>
  mutate(log.net.area = log(net.area), sqrt.b.depth = sqrt(b.depth))
ctrl <- gam.control(nthreads = 3)
m <- gam(egg.count ~ s(lon, lat, k = 100) + s(sqrt.b.depth) + s(c.dist) +
  s(temp.20m) + offset(log.net.area),
  data = mack, family = tw(), method = "REML", control = ctrl)

海洋深度的平滑度

b.depth

基于协变量的平方根，因此与您的设置相匹配。另请注意，响应数据假定为以模型为条件的 Tweedie 分布。

首先我们创建要评估平滑度或模型的数据：

ds <- data_slice(mack, b.depth = evenly(b.depth), net.area = 1) |>
  mutate(sqrt.b.depth = sqrt(b.depth), log.net.are = log(net.area))

在这里，我使用 gratia 函数在海洋深度上均匀地生成数据，然后将其转换为用于拟合模型的比例，并提供偏移项的值（用于对不同大小的海洋进行标准化）数据中使用的网络）。

选项1

在此选项中，我们正在评估协变量请求值的平滑度，但由于该模型涉及非恒等链接函数，我们需要使用链接函数的逆函数反向转换为响应尺度。我们只评估所需的平滑度，然后添加置信区间、模型常数项，然后再对所有内容进行反向转换。我可以从

b.depth

中提取

ds

的原始值，但在这种情况下，直接使用

b.depth

反向转换到

mutate()

比例会更简单。然后我们进行绘图。

ilink <- inv_link(m)
sm <- smooth_estimates(m, "s(sqrt.b.depth)", data = ds) |>
  add_confint() |>
  add_constant(constant = coef(m)[1]) |>
  transform_fun(fun = ilink) |>
  mutate(b.depth = sqrt.b.depth^2)

sm |>
  ggplot(aes(x = b.depth, y = .estimate)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci), alpha = 0.2,
    fill = "steelblue") +
  geom_line() +
  labs(y = "Estimated egg count", x = "Ocean depth (m)")

选项2

在此选项中，我们根据模型进行预测，但由于模型是可加的，我们只需选择要包含在该预测中的项。这里我选择截距（模型常数项）和海洋深度的平滑度。

gratia::fitted_values()

中的默认设置是在响应尺度上生成预测 - 它通过首先在链接尺度上进行预测，计算可信区间，然后使用直接提取的链接函数的逆函数将所有内容转换为响应尺度来实现此目的从模型来看。

fv <- fitted_values(m , data = ds,
  terms = c("(Intercept)", "s(sqrt.b.depth)"))

fv |>
  ggplot(aes(x = b.depth, y = .fitted)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci), alpha = 0.2,
    fill = "steelblue") +
  geom_line() +
  labs(y = "Estimated egg count", x = "Ocean depth (m)")

两种选择都会产生

y 轴刻度现在更容易解释为预期计数，x 轴处于自然刻度上。但解释是一样的；根据模型中的项，该图显示了预期计数如何随海洋深度变化。