Warum ist das Bayes'sche glaubwürdige Intervall in dieser Polynomregression verzerrt, während das Konfidenzintervall korrekt ist?

Betrachten Sie das Diagramm unten, in dem ich Daten wie folgt simuliert habe. Wir betrachten ein binäres Ergebnis für das die wahre Wahrscheinlichkeit, 1 zu sein, durch die schwarze Linie angezeigt wird. Die funktionale Beziehung zwischen einer Kovariate und ist ein Polynom 3. Ordnung mit logistischer Verknüpfung (also in doppelter Hinsicht nicht linear). x p ( y o b s = 1 | x $y_{obs}$$x$$p(y_{obs}=1 | x)$

Die grüne Linie ist die logistische GLM-Regressionsanpassung, bei der als Polynom 3. Ordnung eingeführt wird. Die gestrichelten grünen Linien sind die 95% -Konfidenzintervalle um die Vorhersage , wobei die angepassten Regressionskoeffizienten sind. Ich habe und dafür verwendet.P ( y o b s = 1 | x , β ) $x$$p(y_{obs}=1 | x, \hat{\beta})$$\hat{\beta}$R glmpredict.glm

In ähnlicher Weise ist die Prupellinie der Mittelwert des Seitenzahns mit einem zu 95% glaubwürdigen Intervall für eines Bayes'schen logistischen Regressionsmodells unter Verwendung eines einheitlichen Prior. Ich habe das Paket mit Funktion dafür verwendet (Einstellung gibt den einheitlichen nicht informativen Prior).$p(y_{obs}=1 | x, \beta)$MCMCpackMCMClogitB0=0

Die roten Punkte bezeichnen Beobachtungen im Datensatz, für die , die schwarzen Punkte sind Beobachtungen mit . Beachten Sie, dass wie bei der Klassifizierung / diskreten Analyse üblich aber nicht beobachtet wird.y o b s = 0 y p ( y o b s = 1 | x $y_{obs}=1$$y_{obs}=0$$y$$p(y_{obs}=1 | x)$

Mehrere Dinge können gesehen werden:

1. Ich habe absichtlich simuliert, dass auf der linken Seite spärlich ist. Ich möchte, dass das Vertrauen und das glaubwürdige Intervall hier aufgrund des Mangels an Informationen (Beobachtungen) groß werden.$x$
2. Beide Vorhersagen sind links nach oben voreingenommen. Diese Verzerrung wird durch die vier roten Punkte verursacht, die Beobachtungen bezeichnen, was fälschlicherweise darauf hindeutet, dass die wahre funktionale Form hier hochgehen würde. Der Algorithmus hat nicht genügend Informationen, um zu schließen, dass die wahre Funktionsform nach unten gebogen ist.$y_{obs}=1$
3. Das Konfidenzintervall wird erwartungsgemäß breiter, das glaubwürdige Intervall jedoch nicht . Tatsächlich umfasst das Konfidenzintervall den gesamten Parameterraum, wie es aufgrund fehlender Informationen sein sollte.

Es scheint, dass das glaubwürdige Intervall hier für einen Teil von falsch / zu optimistisch ist . Es ist wirklich unerwünscht, dass das glaubwürdige Intervall eng wird, wenn die Informationen spärlich werden oder vollständig fehlen. Normalerweise reagiert ein glaubwürdiges Intervall nicht so. Kann jemand erklären:$x$

1. Was sind Gründe dafür?
2. Welche Schritte kann ich unternehmen, um ein glaubwürdigeres Intervall zu erreichen? (das heißt, eine, die mindestens die wahre funktionale Form einschließt oder besser so breit wird wie das Konfidenzintervall)

Code zum Abrufen von Vorhersageintervallen in der Grafik wird hier gedruckt:

fit <- glm(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
x_pred <- seq(0, 1, by=0.01)
pred <- predict(fit, newdata = data.frame(x=x_pred), se.fit = T)
plot(plogis(pred$fit), type='l')
matlines(plogis(pred$fit + pred$se.fit %o% c(-1.96,1.96)), type='l', col='black', lty=2)


library(MCMCpack)
mcmcfit <- MCMClogit(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
gibbs_samps <- as.mcmc(mcmcfit)
x_pred_dm <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=x_pred))
gibbs_preds <- apply(gibbs_samps, 1, `%*%`, t(x_pred_dm))
gibbs_pis <- plogis(apply(gibbs_preds, 1, quantile, c(0.025, 0.975)))
matlines(t(gibbs_pis), col='red', lty=2)

Datenzugriff : https://pastebin.com/1H2iXiew Dank @DeltaIV und @AdamO

Für ein frequentistischen Modell, die Varianz der Vorhersage vergrößert im Verhältnis zum Quadrat des Abstandes von dem Schwerpunkt der . Ihre Methode zur Berechnung von Vorhersageintervallen für ein Bayes'sches GLM verwendet empirische Quantile basierend auf der angepassten Wahrscheinlichkeitskurve, berücksichtigt jedoch nicht die Hebelwirkung von$X$$X$

Ein binomialer frequentistischer GLM unterscheidet sich nicht von einem GLM mit Identitätsverknüpfung, außer dass die Varianz proportional zum Mittelwert ist.

Beachten Sie, dass jede Polynomdarstellung von Logit-Wahrscheinlichkeiten zu Risikovorhersagen führt , die je nach Vorzeichen des Ausdrucks der höchsten Polynomordnung gegen 0 als und 1 als oder umgekehrt konvergieren .X → $X\rightarrow -\infty$$X\rightarrow \infty$

Für die häufigere Vorhersage dominiert die proportionale Zunahme der Varianz der Vorhersagen durch quadratische Abweichung (Hebelwirkung) diese Tendenz. Aus diesem Grund ist die Konvergenzrate zu Vorhersageintervallen, die ungefähr [0, 1] entsprechen, schneller als die Polynom-Logit-Konvergenz dritter Ordnung zu Wahrscheinlichkeiten von 0 oder 1.

Dies gilt nicht für Bayes'sche posterior angepasste Quantile. Es gibt keine explizite Verwendung der quadratischen Abweichung, daher verlassen wir uns einfach auf den Anteil dominierender 0- oder 1-Tendenzen, um langfristige Vorhersageintervalle zu konstruieren.

Dies wird durch Extrapolation sehr weit in die Extreme von .$X$

Mit dem oben angegebenen Code erhalten wir:

> x_pred_dom <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=c(1000)))
> gibbs_preds <- plogis(apply(gibbs_samps[1000:10000, ], 1, `%*%`, t(x_pred_dom))) # a bunch of 0/1s basically past machine precision
> prop.table(table(gibbs_preds))
gibbs_preds
         0          1 
0.97733585 0.02266415 
> 

In 97,75% der Fälle war der dritte Polynomterm negativ. Dies ist anhand der Gibbs-Proben überprüfbar:

> prop.table(table(gibbs_samps[, 4]< 0))

 FALSE   TRUE 
0.0225 0.9775 

Daher konvergiert die vorhergesagte Wahrscheinlichkeit gegen 0, wenn gegen unendlich geht. Wenn wir die SEs des Bayes'schen Modells untersuchen, stellen wir fest, dass die Schätzung des dritten Polynomterms -185,25 mit se 108,81 beträgt, was bedeutet, dass es 1,70 SDs von 0 sind. Unter Verwendung normaler Wahrscheinlichkeitsgesetze sollte er in 95,5% der Fälle unter 0 fallen ( keine schrecklich andere Vorhersage basierend auf 10.000 Iterationen). Nur eine andere Art, dieses Phänomen zu verstehen.$X$

Auf der anderen Seite steigt die Frequentist-Passform erwartungsgemäß auf 0,1:

freq <- predict(fit, newdata = data.frame(x=1000), se.fit=T)
plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)

gibt:

> plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)
     [,1]
[1,]    0
[2,]    1