在ggplot中绘制我GAM中的平滑函数

Question

在ggplot中绘制我GAM中的平滑函数

4

我已经创建了一个GAM模型，并设置了预测，但是在如何绘制模型的任何平滑函数方面遇到了困难。一直在尝试在ggplot中绘制这些函数，但现在添加了月份后，使用参数/美学方面有些困难。看到有些人建议也使用geom_smooth()，但我不确定。如果有人能就此向我提供建议，那将是很棒的。下面是我的数据、模型和预测：

模型

mod = gam(co2 ~ s(timeStep, k = 200, bs = "cs") + s(month, k = 12, bs = "cc"), 
                data = carbonD,
                family = gaussian(link = "identity"))

预测

#create predictions
preds = predict(mod, type = 'terms', se.fit = TRUE) 
#combine our predictions with coefficients
fit = preds$fit + coef(mod)[1]

数据片段

carbonD
       co2 month year timeStep
1   315.42     1 1959        1
2   316.31     2 1959        2
3   316.50     3 1959        3
4   317.56     4 1959        4
5   318.13     5 1959        5
6   318.00     6 1959        6
7   316.39     7 1959        7
8   314.65     8 1959        8
9   313.68     9 1959        9
10  313.18    10 1959       10
11  314.66    11 1959       11
12  315.43    12 1959       12
13  316.27     1 1960       13
14  316.81     2 1960       14
15  317.42     3 1960       15

- Joe

2个回答

1

最简单的方法是使用LOESS和geom_smooth： geom_smooth(method="loess", span=0.5)，可以通过调整span参数来获得更加平滑或波浪形状。请注意保留HTML标签。

- asafpr

这是正确的，但没有绘制OP想要的实际gam模型。 - Allan Cameron

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Allan Cameron · Accepted Answer

有两种方法可以在ggplot中绘制您的精确模型。一种是使用geom_smooth，但是右侧有两个变量时无法使用此方法。实际上，在您的情况下，这是可能的，因为月份可以从时间步骤计算出来，但是现在让我们忽略它，直接使用带和线绘制您的模型预测。

首先，加载所需的软件包并创建模型（请注意，由于我们只有数据的片段，因此我不得不减少节点数）。

library(mgcv)
library(ggplot2)

mod = gam(co2 ~ s(timeStep, k = 4, bs = "cs") + s(month, k = 12, bs = "cc"), 
                data = carbonD,
                family = gaussian(link = "identity"))

现在我们创建一个小数据框，其中包含我们想要进行预测的值，在我们数据范围内有1000个点：

newdata <- data.frame(timeStep = seq(1, 15, length.out = 1000),
                      month = (seq(1, 15, length.out = 1000) - 1) %% 12 + 1)

现在我们进行预测，并使用标准误差拟合来创建上限和下限置信带。

pred <- predict(mod, newdata, type = 'response', se.fit = TRUE) 

newdata$co2   <- pred$fit
newdata$lower <- pred$fit - 1.96 * pred$se.fit
newdata$upper <- pred$fit + 1.96 * pred$se.fit

现在我们可以绘制我们的结果：

ggplot(carbonD, aes(timeStep, co2)) +
  geom_point() +
  geom_ribbon(data = newdata, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_minimal(base_size = 16)

您也可以直接在geom_smooth中使用您的gam，但是您需要能够用y和x来表达模型，其中x是时间步长。您可以通过从时间步长中减去1来获取月份，取得该数字模12，并再次加上1，这样可以避免明确创建预测数据框，代价是使绘图代码更加复杂：

ggplot(carbonD, aes(timeStep, co2)) +
  geom_point() +
  geom_smooth(formula = y ~ s(x, k = 4, bs = "cs") + 
                            s((x - 1) %% 12 + 1, k = 12, bs = "cc"),
              method = "gam", size = 0.7,
              method.args = list(family = gaussian(link = "identity")),
              aes(color = "gam", fill = "confidence interval")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_minimal(base_size = 16)

需要说明的是，我不确定你是否应该同时使用月份和时间步长，因为一个只是另一个的模数。如果您想要分离长期和季节性影响，只使用时间步长或使用年和月份可能会更好。