¿Por qué es el intervalo creíble Bayesiano en esta regresión polinomial parcial considerando el intervalo de confianza es correcto?

Question

Considere el gráfico siguiente, en el que yo los datos simulados de la siguiente manera. Nos fijamos en un resultado binario $y_{obs}$ de que la verdadera probabilidad 1 se indica por la línea negra. La relación funcional entre la covariable $x$ $p(y_{obs}=1 | x)$ es 3ª polinomiales de orden logístico de enlace (por lo que no es lineal en un doble sentido).

La línea verde es el GLM ajuste de regresión logística donde la $x$ se presentó como el 3 de la orden de polinomio. El discontinua líneas verdes son el 95% de intervalos de confianza alrededor de la predicción de $p(y_{obs}=1 | x, \hat{\beta})$ donde $\hat{\beta}$ el conjunto de los coeficientes de regresión. Yo solía R glm y predict.glm de este.

Del mismo modo, el pruple línea es la media de la parte posterior con un 95% de intervalo creíble para $p(y_{obs}=1 | x, \beta)$ de un Bayesiano modelo de regresión logística utilizando un uniforme de antes. He utilizado el paquete MCMCpack con la función MCMClogit (ajuste B0=0 le da el uniforme de valor informativo previo).

Los puntos rojos indican las observaciones en el conjunto de datos para que $y_{obs}=1$, los puntos negros son las observaciones con $y_{obs}=0$. Tenga en cuenta que es tan común en la clasificación / discretos análisis de la $y$ pero no $p(y_{obs}=1 | x)$ que se observa.

Varias cosas se pueden ver:

1. He simulado en el propósito de que $x$ es escasa en la mano izquierda. Quiero que la confianza y creíble intervalo de obtener amplia, debido a la falta de información (observaciones).
2. Ambas predicciones son sesgados hacia arriba a la izquierda. Este sesgo es causada por los cuatro puntos rojos denotan $y_{obs}=1$ observaciones, que erróneamente sugiere que la verdadera forma funcional iría hasta aquí. El algoritmo tiene información suficiente para concluir que la verdadera forma funcional es doblado hacia abajo.
3. El intervalo de confianza se hace más ancha, como se esperaba, mientras que el intervalo creíble ¿ no. De hecho, el intervalo de confianza incluye el completo espacio de parámetros, como debería debido a la falta de información.

Parece que el intervalo creíble es malo / muy optimista de aquí para una parte de la $x$. Es realmente un comportamiento indeseable para la credibilidad del intervalo para obtener estrecho cuando la información se presenta escasa o está totalmente ausente. Generalmente esto no es como un intervalo creíble reacciona. Puede alguien explicar:

1. ¿Cuáles son las razones para esto?
2. ¿Qué pasos puedo tomar para llegar a un mayor intervalo creíble? (es decir, la que encierra, al menos, la verdadera forma funcional, o mejor se pone tan amplia como el intervalo de confianza)

Código para obtener los intervalos de predicción en el gráfico se imprimen aquí:

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)

---

Datos de acceso: https://pastebin.com/1H2iXiew gracias @DeltaIV y @AdamO

Answer 1

Para un modelo frecuentista, la varianza de la predicción magnifica en proporción al cuadrado de la distancia desde el centroide de $X$. Su método de cálculo de los intervalos de predicción para un Bayesiano GLM utiliza cuantiles empíricos basados en el amueblada curva de probabilidad, pero no tiene en cuenta para $X$'s de apalancamiento.

Un binomio frecuentista GLM no es diferente de un GLM con identidad enlace, salvo que la varianza es proporcional a la media.

Tenga en cuenta que cualquier polinomio representación de logit probabilidades conduce a las predicciones del riesgo que convergen a 0 $X\rightarrow -\infty$ y 1 $X\rightarrow \infty$ o viceversa, dependiendo del signo de la más alta polinomio de orden plazo.

Para frecuentista de la predicción, el cuadrado de la desviación (apalancamiento) aumento proporcional en la varianza de las predicciones que domina esta tendencia. Esta es la razón por la velocidad de convergencia a intervalos de predicción aproximadamente igual a [0, 1] es más rápido que el de la tercera orden de polinomio logit convergencia de probabilidades de 0 o 1 singularmente.

Esto no es así para Bayesiano posterior de módulos de cuantiles. No hay ningún uso explícito de los cuadrados de la desviación, por lo que nos basamos simplemente en la proporción de dominar 0 o 1 tendencias para construir a largo plazo intervalos de predicción.

Esto se hizo evidente por la extrapolación de muy lejos, hacia los extremos de la $X$.

Utilizando el código que he proporcionado anteriormente obtenemos:

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

Así 97.75% del tiempo, el tercer polinomio plazo fue negativo. Esto es verificable a partir de las muestras de Gibbs:

> prop.table(table(gibbs_samps[, 4]< 0))

 FALSE   TRUE 
0.0225 0.9775 

Por lo tanto la predicción de la probabilidad converge a 0 $X$ va al infinito. Si examinamos la SEs de la modelo Bayesiano, nos encontramos con la estimación de la tercer término polinómico es -185.25 con se 108.81 lo que significa que es de 1.70 SDs desde 0, por lo que el uso de probabilidad normal de las leyes, se debe caer por debajo de 0 95.5% del tiempo (no muy diferente de la predicción basada en 10.000 iteraciones). Otra manera de entender este fenómeno.

Por otro lado, la frecuentista encajar los golpes de hasta el 0,1 como se esperaba:

freq <- predict(fit, newdata = data.frame(x=1000), se.fit=T)
plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)

da:

> plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)
     [,1]
[1,]    0
[2,]    1