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

9

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 )yobsxp(yobs=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 , β ) βxp(yobs=1|x,β^)β^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(yobs=1|x,β)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 )yobs=1yobs=0yp(yobs=1|x)

Geben Sie hier die Bildbeschreibung ein

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.yobs=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

Tomka
quelle
Wenn mir jemand erklären könnte, wie man eine Tabelle mit den Daten teilt, kann ich das tun.
Tomka
Sie können dputden Datenrahmen, der die Daten enthält, verwenden und dann die dputAusgabe als Code in Ihren Beitrag aufnehmen.
DeltaIV
1
@ Tomka oh ich verstehe. Ich bin nicht farbenblind, aber es ist sehr schwierig für mich, den grün / blauen Unterschied zu erkennen!
AdamO
1
@ AdamO hoffe, das ist besser
Tomka

Antworten:

6

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 vonX.XX

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 XX

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
AdamO
quelle
Dennoch: Ist das Bayes'sche Modell in Bereichen der Daten , für die es keine Beispiele gesehen hat, nicht zu zuversichtlich ? Ich weiß, dass Bayes'sche Posterioren oder prädiktive Verteilungen oft ein sehr unterschiedliches Verhalten aufweisen (dh eher dem Konf. Intervall). Ich vermute, dass es einige Auswirkungen des Prior gibt. Wenn Sie manipulieren in Sie geben die Genauigkeit eines normalen vor und kann einen großen Einfluss auf die glaubwürdige Intervall beobachten. xB0MCMClogit
Tomka
@tomka Ich weiß nicht, wie ich das genau beantworten soll, da es tangential zu der vorliegenden Frage zu sein scheint. Das Wichtigste ist, darauf hinzuweisen, dass diese Methoden zur Berechnung von PIs nicht wirklich vergleichbar sind, insbesondere da sie sich auf die Extrapolation beziehen. Wenn Sie mit der Bayes'schen Folgerung einen informativen Prior verwenden, gewinnen Sie natürlich an Effizienz, wenn der Prior richtig ist, und verlieren, wenn der Prior falsch ist.
AdamO
Nur um Sie wissen zu lassen, dass ich immer noch über Ihre Antwort nachdenke. Ich finde es immer noch seltsam, dass der hintere Teil nicht durch Verbreiterung auf die Sparsamkeit reagiert. Ich glaube, dass für andere Priors ein besseres Verhalten in der spärlichen Region erreicht werden kann. Ich kann das im Moment nicht genau sagen; Ich werde die Frage vielleicht mit einem Beispiel erweitern, in dem das glaubwürdige Intervall so funktioniert, wie ich es erwarten würde, selbst im Falle einer Extrapolation (ich denke insbesondere an eine normale lineare Bayes'sche Regression). Wenn ich das tue, werde ich es dich wissen lassen.
Tomka