Przedziały ufności dla parametrów regresji: bayesowska vs. klasyczna

Biorąc pod uwagę dwie tablice x i y o długości n, dopasowuję model y = a + b * x i chcę obliczyć 95% przedział ufności dla nachylenia. Jest to (b - delta, b + delta), gdzie b znajduje się w zwykły sposób i

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

a se.slope to błąd standardowy nachylenia. Jednym ze sposobów uzyskania standardowego błędu nachylenia z R. jest summary(lm(y~x))$coef[2,2].

Przypuśćmy teraz, że piszę prawdopodobieństwo nachylenia podanego xiy, pomnóż to przez „płaski” przedtem i użyj techniki MCMC, aby narysować próbkę m z rozkładu tylnego. Definiować

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

Moje pytanie: czy w (lims[[2]]-lims[[1]])/2przybliżeniu jest równe delcie zdefiniowanej powyżej?

Dodatek Poniżej znajduje się prosty model JAGS, w którym te dwa wydają się różne.

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

Uruchamiam następujące w R:

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

I dostać:

Klasyczny region zaufania: +/- 4,6939

Region zaufania Bayesa: +/- 5.1605

Ponownie uruchamiając to, region zaufania Bayesa jest konsekwentnie szerszy niż region klasyczny. Czy to dlatego, że wybrałem priorytety?

„Problem” dotyczy wcześniej sigmy. Wypróbuj mniej pouczające ustawienie

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

w twoim pliku jags. Następnie zaktualizuj kilka

update(10000)

chwyć parametry i podsumuj swoje zainteresowanie. Powinien dość dobrze pasować do klasycznej wersji.

Wyjaśnienie: Aktualizacja ma na celu upewnienie się, że dotrzesz tam, dokąd zmierzasz, bez względu na to, jaki wybierzesz, chociaż łańcuchy dla modeli takich jak ten z rozproszonymi priorytetami i losowymi wartościami początkowymi zajmują więcej czasu. W prawdziwych problemach sprawdzałbyś zbieżność przed podsumowaniem czegokolwiek, ale nie jest to główny problem w twoim przykładzie.

Jeśli pobierzesz próbkę z tylnej części b | y i oblicz limity (jak zdefiniujesz) powinno być takie same jak (b - delta, b + delta). W szczególności, jeśli obliczysz rozkład tylny b | y pod płaskim przed, jest taki sam jak klasyczny rozkład próbkowania b.

Aby uzyskać więcej informacji, patrz: Gelman i in. (2003). Analiza danych bayesowskich. CRC Press. Sekcja 3.6

Edytować:

Ringold, obserwowane przez ciebie zachowanie jest zgodne z ideą bayesowską. Bayesian Credible Interval (CI) jest ogólnie szerszy niż te klasyczne. A powodem jest, jak słusznie zgadłeś, że hyperpriory wzięły pod uwagę zmienność z powodu nieznanych parametrów.

W przypadku takich prostych scenariuszy (NIE OGÓLNIE):

Baysian CI> Empirical Bayesian CI> Classical CI; > == szerszy

W przypadku liniowych modeli Gaussa lepiej jest użyć pakietu bayesm. Implementuje pół-sprzężoną rodzinę priorów, a przeor Jeffreys jest przypadkiem granicznym tej rodziny. Zobacz mój przykład poniżej. Są to klasyczne symulacje, nie ma potrzeby korzystania z MCMC.

Nie pamiętam, czy przedziały wiarygodności dotyczące parametrów regresji są dokładnie takie same, jak zwykłe przedziały ufności metodą najmniejszych kwadratów, ale w każdym razie są bardzo zbliżone.

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

Biorąc pod uwagę, że prosta regresja liniowa jest analitycznie identyczna między analizą klasyczną i bayesowską z wcześniejszym Jeffreyem, przy czym oba są analityczne, wydaje się nieco dziwne, aby skorzystać z metody numerycznej, takiej jak MCMC, do wykonania analizy bayesowskiej. MCMC to tylko numeryczne narzędzie integracyjne, które pozwala na stosowanie metod bayesowskich w bardziej skomplikowanych problemach, które są trudne do analizy, tak samo jak Newton-Rhapson lub Fisher Scoring to numeryczne metody rozwiązywania klasycznych problemów, które są trudne do rozwiązania.

Rozkład tylny p (b | y) z wykorzystaniem wcześniejszego p (a, b, s) Jeffrey'a proporcjonalnego do 1 / s (gdzie s jest standardowym odchyleniem błędu) jest rozkładem t-studenta z lokalizacją b_ols, skala se_b_ols (" ols „dla„ zwykłych najmniejszych kwadratów ”) i n-2 stopni swobody. Ale rozkład próbkowania b_ols jest również studentem t z lokalizacją b, skalą se_b_ols i n-2 stopniami swobody. Zatem są identyczne, z wyjątkiem tego, że zmieniono b i b_ols, więc jeśli chodzi o tworzenie przedziału, „est + - granica” przedziału ufności zostaje odwrócona do „est - + granica” w wiarygodnym przedziale.

Tak więc przedział ufności i przedział wiarygodności są analitycznie identyczne i nie ma znaczenia, która metoda zostanie użyta (pod warunkiem, że nie ma żadnych dodatkowych wcześniejszych informacji) - więc wybierz metodę, która jest obliczeniowo tańsza (np. Ta z mniejszą inwersją macierzy). Wynik w MCMC pokazuje, że konkretne przybliżenie zastosowane w MCMC daje wiarygodny przedział, który jest zbyt szeroki w porównaniu z dokładnym przedziałem wiarygodności analitycznej. Jest to prawdopodobnie dobra rzecz (chociaż chcielibyśmy, aby przybliżenie było lepsze), że przybliżone rozwiązanie bayesowskie wydaje się bardziej konserwatywne niż dokładne rozwiązanie bayesowskie.