我可以使用自舉來估計 GAM 最大值的不確定性嗎？

我有一個實驗數據，我將藻類生物量的發展視為營養素濃度的函數。生物量（響應變量）和濃度（解釋變量）之間的關係或多或少是單峰的，沿 x 軸有一個明顯的“最佳”，即生物量達到峰值。

我已經為數據擬合了一個廣義加法模型（GAM），並且對給出最高生物量（即對應於峰值 y 值）的濃度（x 軸上的值）感興趣。我從 8 種不同濃度的營養物（總共 24 次觀察）中復制了 3 次生物量。除了只知道 GAM 的峰值濃度之外，我還想獲得某種對最優位置的不確定性估計。在這裡使用引導程序是一種有效的方法嗎？

我嘗試了以下程序：

1. 對於八種營養濃度中的每一種，從三個重複中替換樣本，以便我得到一個包含 24 個觀察值的新數據集，但來自同一重複的數據可能會被多次使用。
2. 將 GAM 擬合到數據並估計函數達到峰值的 x 值。
3. 重複很多次，每次都保存估計值，然後計算標準差（或 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 很慢。引導的另外兩個問題是（取自西蒙的筆記！）

1. 對於參數引導，平滑偏差會導致問題，從中模擬的模型是有偏差的，並且對樣本的擬合將更加有偏差。
2. 對於非參數“案例重採樣”，相同數據的複制副本的存在會導致平滑度不足，尤其是基於 GCV 的平滑度選擇。

應該注意的是，這裡執行的後驗模擬取決於為模型/樣條選擇的平滑度參數。這可以解釋，但西蒙的筆記表明，如果你真的不厭其煩地去做，這沒什麼區別。（所以我沒來……）

引用自：https://stats.stackexchange.com/questions/190348