我有一个实验数据，我将藻类生物量的发展视为营养素浓度的函数。生物量（响应变量）和浓度（解释变量）之间的关系或多或少是单峰的，沿 x 轴有一个明显的“最佳”，即生物量峰值。
我已经为数据拟合了一个广义加法模型（GAM），并且对给出最高生物量（即对应于峰值 y 值）的浓度（x 轴上的值）感兴趣。我从 8 种不同浓度的营养物（总共 24 次观察）中复制了 3 次生物量。除了只知道 GAM 在哪个浓度下达到峰值，我还想获得某种对最优位置的不确定性估计。在这里使用引导程序是一种有效的方法吗？
我尝试了以下程序：
对于八种营养浓度中的每一种，从三个重复中替换样本，以便我得到一个包含 24 个观察值的新数据集，但来自同一重复的数据可能会被多次使用。
将 GAM 拟合到数据并估计函数达到峰值的 x 值。
重复很多次，每次都保存估计值，然后计算标准差（或 2.5 和 97.5 个百分位数）。
这是一种有效的方法还是我在这里遇到麻烦？
可用于使用 Simon Wood 的R 的mgcv软件拟合的 GAM 的另一种方法是从拟合的 GAM 中对感兴趣的特征进行后验推断。从本质上讲，这涉及从拟合模型的参数的后验分布进行模拟，预测在精细网格上的响应值地点，找到其中拟合曲线取其最大值，重复许多模拟模型，并计算最优位置的置信度作为模拟模型中最优分布的分位数。
我在下面介绍的内容摘自 Simon Wood课程笔记的第 4 页(pdf)
为了获得类似于生物量示例的内容，我将使用我的coenocliner包沿单个梯度模拟单个物种的丰度。
library("coenocliner")
A0 <- 9 * 10 # max abundance
m <- 45 # location on gradient of modal abundance
r <- 6 * 10 # species range of occurence on gradient
alpha <- 1.5 # shape parameter
gamma <- 0.5 # shape parameter
locs <- 1:100 # gradient locations
pars <- list(m = m, r = r, alpha = alpha,
gamma = gamma, A0 = A0) # species parameters, in list form
set.seed(1)
mu <- coenocline(locs, responseModel = "beta", params = pars, expectation = FALSE)
适合 GAM
library("mgcv")
m <- gam(mu ~ s(locs), method = "REML", family = "poisson")
...在范围内的精细网格上进行预测( locs)...
p <- data.frame(locs = seq(1, 100, length = 5000))
pp <- predict(m, newdata = p, type = "response")
并可视化拟合函数和数据
plot(mu ~ locs)
lines(pp ~ locs, data = p, col = "red")
这产生
5000 个预测位置在这里可能是多余的，当然对于绘图来说也是如此，但是根据您的用例中的拟合函数，您可能需要一个精细的网格来接近拟合曲线的最大值。
现在我们可以从模型的后面进行模拟。首先我们得到矩阵; 矩阵，一旦乘以模型系数，就可以在新位置从模型中得到预测p
Xp <- predict(m, p, type="lpmatrix") ## map coefs to fitted curves
接下来我们收集拟合的模型系数及其（贝叶斯）协方差矩阵
beta <- coef(m)
Vb <- vcov(m) ## posterior mean and cov of coefs
系数是具有均值向量beta和协方差矩阵的多元法线Vb。因此，我们可以从这个多元正态新系数中模拟与拟合模型一致但探索拟合模型中的不确定性的模型。这里我们生成 10000 ( n)` 个模拟模型
n <- 10000
library("MASS") ## for mvrnorm
set.seed(10)
mrand <- mvrnorm(n, beta, Vb) ## simulate n rep coef vectors from posterior
现在我们可以生成n模拟模型的预测，通过应用链接函数的倒数ilink()（p$locs拟合曲线最大点处的值（的值）
opt <- rep(NA, n)
ilink <- family(m)$linkinv
for (i in seq_len(n)) {
pred <- ilink(Xp %*% mrand[i, ])
opt[i] <- p$locs[which.max(pred)]
}
现在我们使用 10,000 个最优分布的概率分位数计算最优的置信区间，每个模拟模型一个
ci <- quantile(opt, c(.025,.975)) ## get 95% CI
对于这个例子，我们有：
> ci
2.5% 97.5%
39.06321 52.39128
我们可以将此信息添加到前面的图中：
plot(mu ~ locs)
abline(v = p$locs[which.max(pp)], lty = "dashed", col = "grey")
lines(pp ~ locs, data = p, col = "red")
lines(y = rep(0,2), x = ci, col = "blue")
points(y = 0, x = p$locs[which.max(pp)], pch = 16, col = "blue")
产生
正如我们所期望的那样，给定数据/观察值，拟合最优值的间隔是非常不对称的。
Simon 课程笔记的幻灯片 5 说明了为什么这种方法可能比引导法更受欢迎。后验模拟的优点是它很快 - 自举 GAM 很慢。引导的另外两个问题是（取自西蒙的笔记！）
应该注意的是，这里执行的后验模拟取决于为模型/样条选择的平滑度参数。这可以解释，但西蒙的笔记表明，如果你真的不厌其烦地去做，这没什么区别。（所以我没来……）
“重复很多次”。我想你是说你总共有 24 个观察值。手动识别看似异常值的观察结果并手动删除这些观察结果可能更容易。如果不挖掘数据集中可能存在的高杠杆点，可能很难通过如此有限的一组观察来获得最大值的完整特征。我的意思是，高（er）杠杆点是那些可能会扭曲你的结果的点，只有 24 个结果很可能会被选中；这会扭曲最大值的分布。
我不确定您正在监测的浓度是由外部环境因素还是由人类的行为决定的。无论如何，传统的敏感性分析或许会有所帮助。用我的行话来说，敏感性分析意味着将特征的值从最优参数集中的值更改 x%。这将有助于评估，例如，如果我们将某物增加 x%，我们将拥有多少生物量