Los intervalos de confianza para los parámetros de regresión: Bayesiano vs clásica

Question

Dadas dos matrices x e y, ambos de longitud n, I se ajustan a un modelo y = a + b*x y se desea calcular un intervalo de confianza 95% para la pendiente. Esto es (b - delta, b + delta) donde b se encuentra en la forma habitual y

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

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

Supongamos ahora que escribo la probabilidad de que la pendiente dado de x y de y, multiplica esto por un "plano" del antes y el uso de un MCMC técnica 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]])/2 aproximadamente igual a delta como se definió anteriormente?

Anexo a Continuación es una simple ENTRECORTADO modelo 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)
}

Yo ejecute el 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 obtenemos:

Clásica de la región de confianza: +/- 4.6939

Bayesiana de la región de confianza: +/- 5.1605

Volver a ejecutar esto varias veces, el Bayesiana de la región de confianza es consistentemente mayor que la clásica. Entonces esto es debido a los priores he elegido?

Answer 1

Si la muestra de la parte posterior de b | y y calcular lims (como se definen) debe ser igual a (b - delta, b + delta). Específicamente, si se calcula la distribución posterior de b | y en un plano previo, es igual a la clásica distribución de muestreo de la b.

Para obtener más detalles, consulte: Gelman et al. (2003). Bayesiano De Análisis De Datos. CRC Press. En la sección 3.6

Editar:

Ringold, el comportamiento observado por usted es consistente con el Bayesiano idea. El Bayesiano Creíble Intervalo (CI) es generalmente más amplio que el de los clásicos. Y la razón es, como usted bien pensado, el hyperpriors tenerse en cuenta la variabilidad debido a los parámetros desconocidos.

Para escenarios simples como estas (NO EN GENERAL):

Baysian CI > Bayesiano Empírico CI > Clásica CI ; > == más amplio

Answer 2

Lineal de Gauss modelos es mejor utilizar el bayesm paquete. Se implementa la semi-conjugado de la familia de los priores, y el Jeffreys anterior es un caso límite de esta familia. Ver mi ejemplo a continuación. Estos son clásicos simulaciones, no hay necesidad de utilizar MCMC.

No recuerdo si los intervalos de confianza acerca de los parámetros de regresión son exactamente el mismo que el habitual de los mínimos cuadrados de los intervalos de confianza, 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)

Answer 3

Dado que la regresión lineal simple es analíticamente idénticos entre lo clásico y el análisis Bayesiano con Jeffrey previo, ambos de los cuales son analíticos, parece un poco extraño que recurrir a un método numérico como MCMC para hacer el análisis Bayesiano. MCMC es sólo una integración numérica de la herramienta, lo que permite Bayesiano métodos a ser usados en problemas más complicados que son intratables, de la misma manera como el de Newton-Rhapson o de Fisher Scoring son métodos numéricos para la resolución de problemas clásicos que son insolubles.

La distribución posterior p(b|y) usando el Jeffrey previa de p(a,b,s) 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 ("oles" de "mínimos cuadrados ordinarios" estimación), y n-2 grados de libertad. Pero la distribución de muestreo de b_ols es también una t de student con la ubicación b, escala de se_b_ols, y n-2 grados de libertad. Por lo tanto son idénticas, excepto que b y b_ols han intercambiado, así que cuando se trata de crear el intervalo, el "est +- dependiente" del intervalo de confianza se invierte a una "est -+ enlazado" en el intervalo creíble.

Por lo que el intervalo de confianza y creíble intervalo son idénticos, y no importa qué método se utiliza (siempre que no exista adicionales de información previa) - así que el método es computacionalmente más barato (por ejemplo, el uno con menos de la matriz de inversiones). Lo que su resultado con MCMC muestra es que la particular aproximación utilizada con MCMC da un intervalo creíble que es demasiado amplia en comparación con la exacta analítica creíble intervalo. Esta es probablemente una buena cosa (aunque nos gustaría que la aproximación a estar mejor) que la aproximación Bayesiana de la solución parece más conservadora que la exacta solución Bayesiana.

Answer 4

El "problema" es en el estado en sigma. Trate de un menor de informativos de la configuración

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

en su entrecortado archivo. Luego de la actualización de un montón

update(10000)

agarra los parámetros, y resumir la cantidad de interés. Debe de alinearse razonablemente bien con la versión clásica.

Aclaración: La actualización es sólo para asegurarse de que lleguen a su destino independientemente de la opción de antes de que usted decida, a pesar de que las cadenas de modelos como este difusa de los priores y al azar a partir de los valores que tome más tiempo para converger. En problemas reales, te gustaría comprobar la convergencia antes de resumir nada, pero la convergencia no es el problema principal en su ejemplo, no creo.