Konfidenzintervalle für Regressionsparameter: Bayesian vs.

15

Bei zwei Arrays x und y, die beide die Länge n haben, passe ich ein Modell y = a + b * x an und möchte ein 95% -Konfidenzintervall für die Steigung berechnen. Dies ist (b - delta, b + delta), wobei b in üblicher Weise und gefunden wird

delta = qt(0.975,df=n-2)*se.slope

und se.slope ist der Standardfehler in der Steigung. Ein Weg, um den Standardfehler der Steigung von R zu erhalten, ist summary(lm(y~x))$coef[2,2].

Nehmen wir nun an, ich schreibe die Wahrscheinlichkeit der Steigung bei x und y, multipliziere diese mit einem "flachen" Vorgänger und benutze eine MCMC-Technik, um eine Probe m aus der posterioren Verteilung zu ziehen. Definieren

lims = quantile(m,c(0.025,0.975))

Meine Frage: ist (lims[[2]]-lims[[1]])/2etwa gleich Delta wie oben definiert?

Anhang Nachfolgend sehen Sie ein einfaches JAGS-Modell, bei dem diese beiden Modelle unterschiedlich zu sein scheinen.

model {
 for (i in 1:N) {
  y[i] ~ dnorm(mu[i], tau)
  mu[i] <- a + b * x[i]
 }
 a ~ dnorm(0, .00001)
 b ~ dnorm(0, .00001)
 tau <- pow(sigma, -2)
 sigma ~ dunif(0, 100)
}

Ich führe folgendes in R aus:

N <- 10
x <- 1:10
y <- c(30.5,40.6,20.5,59.1,52.5,
       96.0,121.4,78.9,112.1,128.4)
lin <- lm(y~x)

#Calculate delta for a 95% confidence interval on the slope
delta.lm <- qt(0.975,df=N-2)*summary(lin)$coef[2,2]

library('rjags')
jags <- jags.model('example.bug', data = list('x' = x,'y' = y,'N' = N),
                   n.chains = 4,n.adapt = 100)
update(jags, 1000)
params <- jags.samples(jags,c('a', 'b', 'sigma'),7500)
lims <- quantile(params$b,c(0.025,0.975))
delta.bayes <- (lims[[2]]-lims[[1]])/2

cat("Classical confidence region: +/-",round(delta.lm, digits=4),"\n")
cat("Bayesian confidence region:  +/-",round(delta.bayes,digits=4),"\n")

Und bekomme:

Klassischer Vertrauensbereich: +/- 4,6939

Bayesianische Vertrauensregion: +/- 5,1605

Die Bayesianische Vertrauensregion wird mehrmals wiederholt und ist durchweg breiter als die klassische. Liegt das an den Prioren, die ich ausgewählt habe?

Ringold
quelle

Antworten:

9

Das "Problem" ist im Prior auf Sigma. Versuchen Sie es mit einer weniger informativen Einstellung

tau ~ dgamma(1.0E-3,1.0E-3)
sigma <- pow(tau, -1/2)

in deiner Zackendatei. Dann aktualisieren Sie eine Reihe

update(10000)

Nehmen Sie die Parameter zur Hand und fassen Sie die Menge Ihres Interesses zusammen. Es sollte einigermaßen gut mit der klassischen Version übereinstimmen.

Klarstellung: Die Aktualisierung dient nur dazu, sicherzustellen, dass Sie vor Ihrer Entscheidung die gewünschten Ziele erreichen, obwohl die Konvergenz von Ketten für Modelle wie dieses mit diffusen Prioritäten und zufälligen Startwerten länger dauert. Bei echten Problemen würden Sie die Konvergenz überprüfen, bevor Sie etwas zusammenfassen, aber Konvergenz ist nicht das Hauptproblem in Ihrem Beispiel, glaube ich nicht.

Conjugateprior
quelle
@Ringold, was hat funktioniert? Der Prior auf Sigma oder das Update? Oder beides? Hast du sie separat getestet?
Neugierig
sollte sein sigma <- pow(tau, -1/2)odersigma <- 1/sqrt(tau)
Curious
@Tomas, ganz richtig. Tippfehler behoben.
Conjugateprior
Obwohl dies offen gesagt die Ursache für den Unterschied sein könnte, da es sich um den Originalcode handelt ...
conjugateprior
6

Wenn Sie aus dem hinteren Teil von b | y und berechne lims (wie du definierst) sollte es dasselbe sein wie (b - delta, b + delta). Insbesondere, wenn Sie die posteriore Verteilung von b | berechnen y unter einem flachen Vorgänger ist es dasselbe wie die klassische Stichprobenverteilung von b.

Weitere Einzelheiten finden Sie unter: Gelman et al. (2003). Bayesianische Datenanalyse. CRC-Presse. Abschnitt 3.6

Bearbeiten:

Ringold, das von Ihnen beobachtete Verhalten stimmt mit der Bayes'schen Idee überein. Das Bayesian Credible Interval (CI) ist im Allgemeinen breiter als die klassischen. Und der Grund ist, wie Sie richtig erraten haben, dass die Hyperprioren die Variabilität aufgrund der unbekannten Parameter berücksichtigt haben.

Für einfache Szenarien wie diese (NICHT ALLGEMEIN):

Bayesianisches CI> Empirisches Bayesianisches CI> Klassisches CI; > == breiter

suncoolsu
quelle
Ich habe mit JAGS Code hinzugefügt, auf den ich anscheinend eine andere Antwort bekomme. Wo ist mein Fehler? Geschieht das wegen der Vorgesetzten?
Ringold
Jetzt bin ich verwirrt. Zuerst haben Sie gesagt, dass die posteriore Verteilung von b | y unter einem flachen Prior der klassischen Stichprobenverteilung von b entspricht. Dann sagten Sie, das Bayesianische CI sei breiter als das klassische. Wie könnte es breiter sein, wenn die Verteilungen gleich sind?
Ringold
Entschuldigung - ich hätte sagen sollen, was @CP in seinen Kommentaren vorgeschlagen hat. Theoretisch ist b | y unter einem flachen Prior und einem klassischen CI dasselbe, aber Sie können dies in JAGS praktisch nur dann tun, wenn Sie einen sehr sehr diffusen Prior wie CP verwenden und viele MCMC-Iterationen verwenden.
Suncoolsu
Ich habe Ihre Konten zusammengeführt, damit Sie Ihre Fragen bearbeiten und Kommentare hinzufügen können. Bitte registrieren Sie Ihr Konto, indem Sie hier klicken: stats.stackexchange.com/users/login ; Mit Ihrer Google Mail OpenID können Sie dies in wenigen Sekunden erledigen und verlieren hier Ihr Konto nicht mehr.
Danke, ich habe mich angemeldet. Und vielen Dank an diejenigen, die diese Frage beantwortet haben. Ich werde das Gamma vor versuchen.
Ringold
5

Für lineare Gauß-Modelle ist es besser, das Bayes-Paket zu verwenden. Es implementiert die semi-konjugierte Familie der Vorfahren, und der Jeffreys-Prior ist ein Grenzfall für diese Familie. Siehe mein Beispiel unten. Dies sind klassische Simulationen, MCMC ist nicht erforderlich.

Ich erinnere mich nicht, ob die Glaubwürdigkeitsintervalle für die Regressionsparameter genau die gleichen sind wie die üblichen Konfidenzintervalle für kleinste Quadrate, aber auf jeden Fall sehr eng.

> # required package
> library(bayesm)
> # data
> age <- c(35,45,55,65,75)
> tension <- c(114,124,143,158,166)
> y <- tension
> # model matrix
> X <- model.matrix(tension~age)
> # prior parameters
> Theta0 <- c(0,0)
> A0 <- 0.0001*diag(2)
> nu0 <- 0
> sigam0sq <- 0
> # number of simulations
> n.sims <- 5000
> # run posterior simulations
> Data <- list(y=y,X=X)
> Prior <- list(betabar=Theta0, A=A0, nu=nu0, ssq=sigam0sq)
> Mcmc <- list(R=n.sims)
> bayesian.reg <- runireg(Data, Prior, Mcmc)
> beta.sims <- t(bayesian.reg$betadraw) # transpose of bayesian.reg$betadraw
> sigmasq.sims <- bayesian.reg$sigmasqdraw
> apply(beta.sims, 1, quantile, probs = c(0.025, 0.975))
[,1] [,2]
2.5% 53.33948 1.170794
97.5% 77.23371 1.585798
> # to be compared with: 
> frequentist.reg <- lm(tension~age)
Stéphane Laurent
quelle
3

Angesichts der Tatsache, dass die einfache lineare Regression analytisch identisch ist zwischen der klassischen und der Bayes'schen Analyse mit der vorherigen von Jeffrey, die beide analytisch sind, erscheint es etwas seltsam, für die Bayes'sche Analyse auf eine numerische Methode wie MCMC zurückzugreifen. MCMC ist nur ein numerisches Integrationswerkzeug, mit dem Bayes'sche Methoden bei komplizierteren Problemen eingesetzt werden können, die analytisch nicht zu lösen sind, genauso wie Newton-Rhapson oder Fisher Scoring numerische Methoden zur Lösung klassischer Probleme, die nicht zu lösen sind.

Die posteriore Verteilung p (b | y) unter Verwendung des vorherigen p (a, b, s) von Jeffrey proportional zu 1 / s (wobei s die Standardabweichung des Fehlers ist) ist eine studentische t-Verteilung mit der Position b_ols, Skala se_b_ols ( ols "für" gewöhnliche kleinste Quadrate "(Schätzung) und n-2 Freiheitsgrade. Die Stichprobenverteilung von b_ols ist aber auch ein Student t mit Position b, Skala se_b_ols und n-2 Freiheitsgraden. Sie sind also identisch, mit der Ausnahme, dass b und b_ols vertauscht wurden. Wenn es also um die Erstellung des Intervalls geht, wird die "est + - Grenze" des Konfidenzintervalls in dem glaubwürdigen Intervall zu einer "est - + Grenze" umgekehrt.

Das Konfidenzintervall und das glaubwürdige Intervall sind also analytisch identisch, und es spielt keine Rolle, welche Methode verwendet wird (vorausgesetzt, es gibt keine zusätzlichen vorherigen Informationen). Nehmen Sie also die Methode, die rechnerisch billiger ist (z. B. die mit weniger Matrixinversionen). Was Ihr Ergebnis mit MCMC zeigt, ist, dass die mit MCMC verwendete Näherung ein glaubwürdiges Intervall ergibt, das im Vergleich zum genauen analytischen glaubwürdigen Intervall zu breit ist. Dies ist wahrscheinlich eine gute Sache (obwohl wir eine bessere Approximation wünschen würden), dass die ungefähre Bayes'sche Lösung konservativer erscheint als die exakte Bayes'sche Lösung.

Wahrscheinlichkeitslogik
quelle
Nicht so seltsam. Ein Grund für die Verwendung einer numerischen Methode, um eine Antwort auf ein Problem zu finden, das analytisch gelöst werden kann, besteht darin, sicherzustellen, dass die Software korrekt verwendet wird.
Ringold
1
Ein weiterer Grund für die Verwendung von Simulationen: Sie erhalten auch hintere Simulationen für jede Funktion der Parameter. Zum Beispiel nutze ich diese Möglichkeit, um ein Glaubwürdigkeitsintervall über die Wahrscheinlichkeit für einen Wert der Kovariate zu erhalten. Ich weiß nicht, wie ich ein häufiges Konfidenzintervall für diese Wahrscheinlichkeit erhalten soll. Pr ( Y > 10 x ) xf(β0,β1,,βp,σ)Pr(Y>10x)x
Stéphane Laurent