Selección de modelo bayesiano e intervalo creíble

Tengo un conjunto de datos con tres variables, donde todas las variables son cuantitativas. Vamos a llamarlo , y . Estoy ajustando un modelo de regresión en una perspectiva bayesiana a través de MCMC conx 1 x $y$$x_1$$x_2$rjags

Hice un análisis exploratorio y el diagrama de dispersión de sugiere que se debe usar un término cuadrático. Luego instalé dos modelos$y\times x_2$

(1)$y=\beta_0+\beta_1*x_1+\beta_2*x_2$

(2)$y=\beta_0+\beta_1*x1+\beta_2*x_2+\beta_3*x_1x_2+\beta_4*x_1^2+\beta_5*x_2^2$

En el modelo 1, el tamaño del efecto de cada parámetro no es pequeño y el intervalo creíble del 95% no contiene el valor .$0$

En el modelo 2, el tamaño del efecto de los parámetros y es pequeño y cada uno de los intervalos creíbles para todos los parámetros contiene .β 4$\beta_3$$\beta_4$$0$

¿El hecho de que un intervalo creíble contenga es suficiente para decir que el parámetro no es significativo?$0$

Luego ajusté el siguiente modelo

(3) $y=\beta_0+\beta_1*x_1+\beta_2*x_2+\beta_3*x^2_2$

El tamaño del efecto de cada parámetro no es pequeño, pero con la excepción de $\beta_1$ todos los intervalos creíbles contienen $0$ .

¿Cuál es la forma correcta de hacer una selección de variables en las estadísticas bayesianas?

EDITAR: ¿ Puedo usar Lasso en cualquier modelo de regresión, como el modelo Beta? Estoy usando un modelo con dispersión variable donde donde es un vector. ¿Debo usar Laplace antes en también?

EDIT2: instalé dos modelos, uno con gaussiana para , y otro con Laplace (doble exponencial).δ $\beta_j$$\delta_j$

Las estimaciones para el modelo gaussiano son

            Mean      SD  Naive SE Time-series SE
B[1]     -1.17767 0.07112 0.0007497      0.0007498
B[2]     -0.15624 0.03916 0.0004128      0.0004249
B[3]      0.15600 0.05500 0.0005797      0.0005889
B[4]      0.07682 0.04720 0.0004975      0.0005209
delta[1] -3.42286 0.32934 0.0034715      0.0034712
delta[2]  0.06329 0.27480 0.0028966      0.0028969
delta[3]  1.06856 0.34547 0.0036416      0.0036202
delta[4] -0.32392 0.26944 0.0028401      0.0028138

Los estimativos para el modelo de lazo son

              Mean      SD  Naive SE Time-series SE
B[1]     -1.143644 0.07040 0.0007421      0.0007422
B[2]     -0.160541 0.05341 0.0005630      0.0005631
B[3]      0.137026 0.05642 0.0005947      0.0005897
B[4]      0.046538 0.04770 0.0005028      0.0005134
delta[1] -3.569151 0.27840 0.0029346      0.0029575
delta[2] -0.004544 0.15920 0.0016781      0.0016786
delta[3]  0.411220 0.33422 0.0035230      0.0035629
delta[4] -0.034870 0.16225 0.0017103      0.0017103
lambda    7.269359 5.45714 0.0575233      0.0592808

Las estimaciones para y redujeron mucho en el modelo Lasso, ¿significa que debería eliminar estas variables del modelo?$\delta_2$$\delta_4$

EDITAR3: El modelo con doble exponencial anterior (Lasso) me da valores de desviación, BIC y DIC más grandes que el modelo con anteriores gaussianos e incluso obtengo valores más pequeños después de eliminar el coeficiente de dispersión en el modelo gaussiano.$\delta_2$

Es bien sabido que construir un modelo basado en lo que es significativo (o algún otro criterio como AIC, si un intervalo creíble contiene 0, etc.) es bastante problemático, particularmente si luego hace inferencia como si no hubiera hecho la construcción del modelo. Hacer un análisis bayesiano no cambia eso (ver también https://stats.stackexchange.com/a/201931/86652 ). Es decir, no debe hacer una selección de variables, sino modelar el promedio (o algo que podría obtener algunos coeficientes cero, pero refleja todo el proceso de modelado, como LASSO o red elástica).

La elección del modelo bayesiano se enmarca más típicamente como el promedio del modelo bayesiano. Tienes diferentes modelos, cada uno con una probabilidad previa diferente. Si la probabilidad del modelo posterior para un modelo es lo suficientemente baja, esencialmente está descartando el modelo por completo. Para pesos anteriores iguales para cada modelo y anteriores planos, el promedio del modelo con pesos proporcionales a para cada modelo se aproxima a esto.$\exp(-\text{BIC}/2)$

Alternativamente, puede expresar el promedio del modelo como un previo que es una mezcla entre una masa puntual (el peso de la masa puntual es la probabilidad previa de que el efecto sea exactamente cero = el efecto no está en el modelo) y una distribución continua (p. Ej. puntas y losas anteriores). El muestreo de MCMC puede ser bastante difícil para tal previo.

Carvalho y col. Motive la contracción de la herradura antes sugiriendo que funciona como una aproximación continua a una espiga anterior. También es un caso de integrar el problema en un modelo jerárquico, donde hasta cierto punto el tamaño y la presencia de efectos en algunas variables relajan un poco la evidencia requerida para otras (a través del parámetro de contracción global, esto es un poco como falso descubrimiento control de frecuencia) y, por otro lado, permiten que los efectos individuales se mantengan solos si la evidencia es lo suficientemente clara. Hay una implementación conveniente disponible desde el paquete brms R que se basa en Stan / rstan . Hay una serie de otros antecedentes similares, como la herradura + prior y todo el tema es un área de investigación en curso.

Existen varios métodos formales para la selección de variables bayesianas. Una revisión un poco desactualizada de los métodos de selección de variables bayesianas se presenta en:

Una revisión de los métodos de selección de variables bayesianas: qué, cómo y qué

Una revisión más reciente, que también incluye una comparación de diferentes métodos y el rendimiento de los paquetes R donde se implementan es:

Métodos y herramientas para la selección de variables bayesianas y el promedio de modelos en regresión lineal univariante

Esta referencia es particularmente útil porque apunta a paquetes R específicos donde solo necesita conectar la respuesta y los valores covariables (y en algunos casos los valores de hiperparámetro) para ejecutar la selección de variables.

Otra forma, rápida, sucia y no recomendada, de realizar la selección de variables "bayesianas" es usar la selección por pasos (hacia adelante, hacia atrás, ambas) usando BIC y el comando R stepAIC (), que se puede ajustar para realizar la selección en términos de BIC

https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/stepAIC.html

Otra forma rápida y sucia de probar es usando la relación de densidad Savage-Dickey y la simulación posterior que ya tienes:$\beta_4=0$

https://arxiv.org/pdf/0910.1452.pdf

$β$

También puede encontrar una buena conferencia sobre el tema aquí a las 41:55:

https://vimeo.com/14553953