Kann ich Bootstrapping verwenden, um die Unsicherheit in einem Maximalwert eines GAM abzuschätzen?

Ich habe Daten aus einem Experiment, in dem ich die Entwicklung von Algenbiomasse als Funktion der Konzentration eines Nährstoffs untersuche. Die Beziehung zwischen Biomasse (der Antwortvariablen) und der Konzentration (der erklärenden Variablen) ist mehr oder weniger unimodal, mit einem klaren "Optimum" entlang der x-Achse, wo die Biomasse ihren Höhepunkt erreicht.

Ich habe ein generalisiertes additives Modell (GAM) an die Daten angepasst und interessiere mich für die Konzentration (den Wert auf der x-Achse), die die höchste Biomasse ergibt (dh dem maximalen y-Wert entspricht). Ich habe drei Replikate von Biomasse aus jeweils 8 verschiedenen Konzentrationen des Nährstoffs (insgesamt 24 Beobachtungen). Ich möchte nicht nur wissen, bei welcher Konzentration das GAM seinen Höhepunkt erreicht, sondern auch eine Art Unsicherheitsschätzung der Position des Optimums erhalten. Ist es ein gültiger Ansatz, hier Bootstrapping zu verwenden?

Ich habe das folgende Verfahren versucht:

1. Für jede der acht Nährstoffkonzentrationen Probe mit Ersatz aus den drei Replikaten, damit ich einen neuen Datensatz mit 24 Beobachtungen erhalte, wobei jedoch Daten aus demselben Replikat mehr als einmal verwendet werden können.

2. Passen Sie ein GAM an die Daten an und schätzen Sie den x-Wert, an dem die Funktion ihren Höhepunkt erreicht.

3. Wiederholen Sie dies mehrmals, speichern Sie die Schätzung jedes Mal und berechnen Sie die Standardabweichung (oder die 2,5- und 97,5-Perzentile).

Ist das ein gültiger Ansatz oder bin ich hier in Schwierigkeiten?

Ein alternativer Ansatz, der für GAMs verwendet werden kann, die mit der mgcv- Software von Simon Wood für R angepasst wurden, besteht darin, aus dem angepassten GAM eine posteriore Inferenz für das interessierende Merkmal durchzuführen . Dies beinhaltet im Wesentlichen die Simulation der posterioren Verteilung der Parameter des angepassten Modells, wobei Werte der Antwort über ein feines Gitter von vorhergesagt werden$x$ Standorte, die finden $x$ Wenn die angepasste Kurve ihren Maximalwert annimmt, wiederholen Sie dies für viele simulierte Modelle und berechnen Sie ein Vertrauen für den Ort der Optima als Quantile der Verteilung der Optima aus den simulierten Modellen.

Das Fleisch von dem, was ich unten präsentieren wurde von Seite 4 von Simon Wood cribbed Skriptum (pdf)

Um ein Beispiel für Biomasse zu erhalten, simuliere ich mit meinem Coenocliner- Paket die Häufigkeit einer einzelnen Art entlang eines einzelnen Gradienten .

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)

Passen Sie das GAM an

library("mgcv")
m <- gam(mu ~ s(locs), method = "REML", family = "poisson")

... auf einem feinen Gitter über den Bereich von vorhersagen $x$( locs) ...

p  <- data.frame(locs = seq(1, 100, length = 5000))
pp <- predict(m, newdata = p, type = "response")

und visualisieren Sie die angepasste Funktion und die Daten

plot(mu ~ locs)
lines(pp ~ locs, data = p, col = "red")

Dies erzeugt

Die 5000 Vorhersagepositionen sind hier wahrscheinlich übertrieben und sicherlich für die Darstellung, aber abhängig von der angepassten Funktion in Ihrem Anwendungsfall benötigen Sie möglicherweise ein feines Gitter, um sich dem Maximum der angepassten Kurve zu nähern.

Jetzt können wir vom hinteren Teil des Modells aus simulieren. Zuerst bekommen wir die$Xp$Matrix; Die Matrix, die einmal mit Modellkoeffizienten multipliziert wird, liefert Vorhersagen aus dem Modell an neuen Ortenp

Xp <- predict(m, p, type="lpmatrix") ## map coefs to fitted curves

Als nächstes sammeln wir die angepassten Modellkoeffizienten und ihre (Bayes'sche) Kovarianzmatrix

beta <- coef(m)
Vb   <- vcov(m) ## posterior mean and cov of coefs

Die Koeffizienten sind eine multivariate Normalen mit mittlerem Vektor betaund Kovarianzmatrix Vb. Daher können wir aus diesem multivariaten normalen neuen Koeffizienten für Modelle simulieren, die mit dem angepassten übereinstimmen, aber die Unsicherheit im angepassten Modell untersuchen. Hier erzeugen wir 10000 ( n) `simulierte Modelle

n <- 10000
library("MASS") ## for mvrnorm
set.seed(10)
mrand <- mvrnorm(n, beta, Vb) ## simulate n rep coef vectors from posterior

Jetzt können wir Vorhersagen für die nsimulierten Modelle generieren , von der Skala des linearen Prädiktors zur Antwortskala transformieren, indem wir die Umkehrung der Verknüpfungsfunktion ( ilink()) anwenden und dann die berechnen$x$Wert (Wert von p$locs) am Maximalpunkt der angepassten Kurve

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)]
}

Nun berechnen wir das Konfidenzintervall für die Optima unter Verwendung von Wahrscheinlichkeitsquantilen der Verteilung von 10.000 Optima, eines pro simuliertem Modell

ci <- quantile(opt, c(.025,.975)) ## get 95% CI

Für dieses Beispiel haben wir:

> ci
    2.5%    97.5% 
39.06321 52.39128

Wir können diese Informationen zum früheren Plot hinzufügen:

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")

was produziert

Wie wir es angesichts der Daten / Beobachtungen erwarten würden, ist das Intervall auf den angepassten Optima ziemlich asymmetrisch.

Folie 5 von Simons Kursnotizen zeigt, warum dieser Ansatz dem Bootstrapping vorgezogen werden könnte. Vorteile der posterioren Simulation sind, dass sie schnell ist - das Bootstrapping von GAMs ist langsam. Zwei zusätzliche Probleme beim Bootstrapping sind (aus Simons Notizen!)

1. Für das parametrische Bootstrapping verursacht die Glättungsverzerrung Probleme, das simulierte Modell ist vorgespannt und die Anpassungen an die Stichproben sind noch stärker vorgespannt.
2. Bei nicht parametrischem "Case-Resampling" führt das Vorhandensein von Replikatkopien derselben Daten zu einer Unterglättung, insbesondere bei der GCV-basierten Glättungsauswahl.

Es ist zu beachten, dass die hier durchgeführte posteriore Simulation von den gewählten Glättungsparametern für das Modell / den Spline abhängig ist. Dies kann erklärt werden, aber Simons Notizen legen nahe, dass dies kaum einen Unterschied macht, wenn Sie sich tatsächlich die Mühe machen, dies zu tun. (also habe ich nicht hier ...)

"Wiederholen Sie eine große Anzahl von Malen". Ich denke, Sie sagen, Sie haben insgesamt 24 Beobachtungen. Es ist wahrscheinlich einfacher, die Beobachtungen, die als Ausreißer erscheinen, von Hand zu identifizieren und von Hand zu entfernen. Es kann schwierig sein, mit einer derart begrenzten Anzahl von Beobachtungen die vollständige Charakterisierung des Maximums zu erreichen, ohne die hohen Hebelpunkte auszubaggern, die Sie möglicherweise in Ihrem Datensatz haben. Was ich damit sagen will, ist, dass hohe (er) Hebelpunkte diejenigen sind, die Ihre Ergebnisse verzerren könnten und nur 24 Ergebnisse haben, was eine ziemlich hohe Wahrscheinlichkeit ergibt, dass diese ausgewählt werden; und das würde die Verteilung der max.

Ich bin nicht sicher, ob die von Ihnen überwachten Konzentrationen durch externe Umweltfaktoren oder durch die Handlung eines Menschen bestimmt werden. In jedem Fall würde vielleicht eine traditionelle Sensitivitätsanalyse helfen. Sensitivitätsanalyse bedeutet in meinem Jargon, den Wert eines Features um x% gegenüber dem Wert im optimalen Parametersatz zu ändern. Das würde helfen zu beurteilen, z. B. wie viel mehr Biomasse wir hätten, wenn wir etwas um x% erhöhen würden