我将使用用于说明 GAM 的经典 4 项数据集进行说明，但将仅模拟来自强非线性项的数据，因为使用单个协变量很容易可视化该过程。$f(x_2)$

library('mgcv')
set.seed(20)

f2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * x)^3 * (1 - x)^10
ysim <- function(n = 500, scale = 2) {
  x <- runif(n)
  e <- rnorm(n, 0, scale)
  f <- f2(x)
  y <- f + e
  data.frame(y = y, x2 = x, f2 = f)
}

df <- ysim()
head(df)

拟合模型：

m <- gam(y ~ s(x2), data = df)

现在，根据估计的样条系数和平滑度参数从响应的后验分布中模拟 50 个新数据点。

set.seed(10)
nsim <- 50
drawx <- function(x, n) runif(n, min = min(x), max = max(x))
newx <- data.frame(x2 = drawx(df[, 'x2'], n = nsim))

在这里，我在的范围内随机均匀地绘制 50 个值。$x_2$

模拟值由下式给出

其中是$\hat{y}$

而是残差标准误差，它通常是根据高斯模型的残差估计的，而不是直接建模的（尽管后者当然是可能的）。我们通过 得到\。$\hat{\sigma}$$\hat{y}$predict()

首先从模型中获取估计的残差

sig2 <- m$sig2

接下来预测 50 个新位置

set.seed(10)
newx <- transform(newx, newy = predict(m, newx, type = "response"))

然后，使用 in的这些值，模拟 50 个来自高斯的新值，其均值和标准差 (sqrt of )$\hat{y}$newynewysig2

newx <- transform(newx, ysim = rnorm(nsim, mean = newy, sd = sqrt(sig2)))

请注意，每个模拟值都是从高斯分布中随机抽取的，其平均值等于模型估计值和方差。$\hat{\sigma}^2$

我们应该检查我们做了什么。首先，从拟合样条生成更多值，以便我们可以看到估计的模型

pred <- with(df, data.frame(x2 = seq(min(x2), max(x2), length = 500)))
pred <- transform(pred, fitted = predict(m, newdata = pred,
                                         type = "response"))

现在我们可以把这一切放在一起

library('ggplot2')
theme_set(theme_bw())

ggplot(df) +
  geom_point(aes(x = x2, y = y), colour = "grey") +
  geom_line(aes(x = x2, y = f2), colour = "forestgreen", size = 1.3) +
  geom_line(aes(x = x2, y = fitted), data = pred, colour = "blue") +
  geom_point(aes(x = x2, y = newy), data = newx, colour = "blue", size = 2) +
  geom_point(aes(x = x2, y = ysim), data = newx, colour = "red",
             size = 2)

粗绿线是真正的函数。灰点是拟合 GAM 的原始数据。蓝线是拟合平滑，估计值。蓝点是我们从模型中模拟的 50 个位置的拟合值。红点是根据估计模型的后验分布中随机抽取的。$f(x_2)$$\hat{f}(x_2)$$y$

这种方法扩展到了多个协变量——只是很难想象正在发生的事情。这也扩展到非高斯响应，例如泊松或二项分布响应。

关键点是：

1. 使用 , 从模型中预测predict()得到响应的期望

   • （这是高斯的平均值，但是泊松的预期计数）

2. 预测或提取响应条件分布所需的任何其他参数的值

   • （对于高斯，如这里所示，有第二个参数，需要完全指定响应的条件分布。对于泊松，例如，只有“均值”，预期计数.)$\sigma^2$

3. 从响应的条件分布模拟

   • （这里我们可以rnorm()直接使用，因为通过 GLM/GAM 估计的分布参数与预期的形式完全相同rnorm()。泊松分布和二项分布也是如此。对于其他分布，情况可能并非如此。例如，带有 Gamma 系列的 GLM/GAM 的参数在您可以将它们插入函数glm()之前需要进行一些转换。gam()rgamma()