Intervalos de confianza para los parámetros de regresión: Bayesiano vs. clásico

Dadas dos matrices x e y, ambas de longitud n, calzo un modelo y = a + b * x y quiero calcular un intervalo de confianza del 95% para la pendiente. Esto es (b - delta, b + delta) donde b se encuentra de la manera habitual y

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

y se.slope es el error estándar en la pendiente. Una forma de obtener el error estándar de la pendiente de R es summary(lm(y~x))$coef[2,2].

Ahora supongamos que escribo la probabilidad de la pendiente dada x e y, multiplique esto por un "plano" anterior y use una técnica MCMC para extraer una muestra m de la distribución posterior. Definir

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

Mi pregunta: ¿es (lims[[2]]-lims[[1]])/2aproximadamente igual a delta como se definió anteriormente?

Anexo A continuación se muestra un modelo JAGS simple donde estos dos parecen ser diferentes.

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

Ejecuto lo siguiente en 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")

Y obten:

Región de confianza clásica: +/- 4.6939

Región de confianza bayesiana: +/- 5.1605

Al volver a ejecutar esto varias veces, la región de confianza bayesiana es consistentemente más amplia que la clásica. Entonces, ¿esto se debe a los antecedentes que he elegido?

El 'problema' está en el sigma anterior. Pruebe con una configuración menos informativa

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

en tu archivo jags. Luego actualiza un montón

update(10000)

tome los parámetros y resuma su cantidad de interés. Debería alinearse razonablemente bien con la versión clásica.

Aclaración: la actualización es solo para asegurarse de llegar a donde vaya, independientemente de la elección que elija, aunque las cadenas para modelos como este con anteriores difusos y valores iniciales aleatorios tardan más en converger. En problemas reales, verificaría la convergencia antes de resumir algo, pero la convergencia no es el problema principal en su ejemplo, no creo.

Si toma muestras de la parte posterior de b | y y calcular lims (como usted define) debe ser igual a (b - delta, b + delta). Específicamente, si calcula la distribución posterior de b | y bajo un plano anterior, es igual a la distribución de muestreo clásica de b.

Para más detalles, consulte: Gelman et al. (2003) Análisis de datos bayesianos. CRC Press. Sección 3.6

Editar:

Ringold, el comportamiento observado por usted es consistente con la idea bayesiana. El intervalo creíble bayesiano (IC) es generalmente más amplio que los clásicos. Y la razón es, como acertadamente adivinó, los hiperprecios tomaron en cuenta la variabilidad debido a los parámetros desconocidos.

Para escenarios simples como estos (NO EN GENERAL):

CI baysiano> CI empírico bayesiano> CI clásico; > == más ancho

Para los modelos lineales gaussianos, es mejor usar el paquete bayesm. Implementa la familia semi-conjugada de anteriores, y la anterior de Jeffrey es un caso límite de esta familia. Vea mi ejemplo a continuación. Estas son simulaciones clásicas, no hay necesidad de usar MCMC.

No recuerdo si los intervalos de credibilidad sobre los parámetros de regresión son exactamente los mismos que los intervalos de confianza de mínimos cuadrados habituales, pero en cualquier caso están muy cerca.

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

Dado que la regresión lineal simple es analíticamente idéntica entre el análisis clásico y bayesiano con el anterior de Jeffrey, ambos analíticos, parece un poco extraño recurrir a un método numérico como MCMC para hacer el análisis bayesiano. MCMC es solo una herramienta de integración numérica, que permite utilizar métodos bayesianos en problemas más complicados que son analíticamente intratables, al igual que Newton-Rhapson o Fisher Scoring son métodos numéricos para resolver problemas clásicos que son intratables.

La distribución posterior p (b | y) que usa la p (a, b, s) anterior de Jeffrey proporcional a 1 / s (donde s es la desviación estándar del error) es una distribución t de Student con ubicación b_ols, escala se_b_ols (" ols "para la estimación de" mínimos cuadrados ordinarios ") y n-2 grados de libertad. Pero la distribución de muestreo de b_ols también es una t de estudiante con ubicación b, escala se_b_ols y n-2 grados de libertad. Por lo tanto, son idénticos, excepto que b y b_ols se han intercambiado, por lo que cuando se trata de crear el intervalo, el "límite est + -" del intervalo de confianza se invierte a un "límite est - +" en el intervalo creíble.

Por lo tanto, el intervalo de confianza y el intervalo creíble son analíticamente idénticos, y no importa qué método se use (siempre que no haya información previa adicional), así que tome el método que es computacionalmente más barato (por ejemplo, el que tiene menos inversiones de matriz). Lo que muestra su resultado con MCMC es que la aproximación particular utilizada con MCMC proporciona un intervalo creíble que es demasiado amplio en comparación con el intervalo analítico creíble exacto. Probablemente sea algo bueno (aunque nos gustaría que la aproximación sea mejor) que la solución bayesiana aproximada parezca más conservadora que la solución bayesiana exacta.