Regresión lineal bajo el enfoque Bayesiano

1. Contenido

Introducción
Regresión lineal clásica
Regresión lineal Bayesiana
Criterios para el análisis Bayesiano
Conclusión
Referencias

2. Objetivos del taller

Revisar la teoría sobre la estadística bayesiana.
Revisar conceptos de regresión lineal bajo el enfoque clásico y bayesiano.
Aplicar los modelos de regresión a la data BostonHousing usando R.

3. Introducción

Para implementar los modelos estadísticos existen dos enfoques o marcos generales:

Enfoque frecuentista: Los datos muestreados de la población se consideran aleatorios y los valores de los parámetros de la población, conocidos como hipótesis nulas, como fijos (pero desconocidos). Para estimar así esta hipótesis nula buscamos los parámetros muestrales que maximicen la probabilidad de los datos.
El enfoque bayesiano: Este enfoque se basa en el teorema de Bayes, por ejemplo, si tenemos un parámetro \(\theta\) de una población y tenemos algunos datos muestreados al azar de esta población \(D\), la probabilidad a posteriori será entonces:

\[\begin{equation} \overbrace{p(\theta/D)}^{Posterior}=\frac{\overbrace{p(D/\theta)}^{Likelihood}.\overbrace{p(\theta)}^{Prior}}{\underbrace{p(D)}_{Evidence}} \end{equation}\]

La \(p(D)\) es la probabilidad de los datos disponibles independientemente del parámetro \(\theta\).

4. Preparación de datos

Para simplificar, usaremos los datos de BostonHousing de la librería mlbench. Para obtener más detalles sobre estos datos, ejecute este comando ?BostonHousing, después de llamar a la librería.

Los datos de BostonHousingconsite en datos de vivienda para 506 secciones censales de Boston del censo de 1970. El marco de BostonHousingdatos contiene los datos originales de Harrison y Rubinfeld (1979). Los datos originales son 506 observaciones sobre 14 variables, siendo medv la variable objetivo:

Las variables mas importantes son:

chas: Variable ficticia de Charles River (= 1 si el tramo limita con el río; 0 en caso contrario).
dis: Distancias ponderadas a cinco centros de empleo de Boston.
medv: Valor medio corregido de viviendas ocupadas por los propietarios.
age: Proporción de unidades ocupadas por sus propietarios construidas antes de 1940.

library(mlbench)
data("BostonHousing")
str(BostonHousing)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : num  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ b      : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Para comprender bien cómo funciona la regresión lineal bayesiana, usaremos solo tres características, dos variables numéricas age , dis y una categórica chas, con la variable objetivo medv que es el valor mediano de las viviendas ocupadas por sus propietarios.

bost <- BostonHousing[,c("medv","age","dis","chas")]
summary(bost)

##       medv            age              dis         chas   
##  Min.   : 5.00   Min.   :  2.90   Min.   : 1.130   0:471  
##  1st Qu.:17.02   1st Qu.: 45.02   1st Qu.: 2.100   1: 35  
##  Median :21.20   Median : 77.50   Median : 3.207          
##  Mean   :22.53   Mean   : 68.57   Mean   : 3.795          
##  3rd Qu.:25.00   3rd Qu.: 94.08   3rd Qu.: 5.188          
##  Max.   :50.00   Max.   :100.00   Max.   :12.127

5. Regresión lineal clásica

La regresión lineal es el modelo básico para el modelamiento estadístico que consiste en determinar relaciones de dependencia de tipo lineal entre una variable dependiente o respuesta, respecto de una o varias variables explicativas o independiente.

\[\begin{equation} y=\beta_1 + \beta_2 X_{2i} +\beta_3 X_{3i} +...+\beta_k X_{ki} + \epsilon_i \end{equation}\]

Donde, asumimos que el error aleatorio \(\epsilon_i\) es IID (Independiente e Idénticamente Distribuida) con media cero y varianza contante \(\sigma^2\), es decir,

\[\begin{equation} \epsilon_i \mathrel{\mathop{\sim}\limits^{\rm iid}}\textsf{Normal}(0, \sigma^2). \end{equation}\]

Para la estimación de los parámetro se hace uso del método de Mínimos Cuadrados Ordinarios (MCO), el método de Máxima Verosimilitud, entre otros. Para el ajuste del modelo se puede hacer uso de la función lmdel R.

Para resaltar la diferencia entre la regresión bayesiana y la regresión lineal clásica (enfoque frecuentista), primero ajustemos esta última a nuestros datos.

model_freq <- lm(medv~., data=bost)
summary(model_freq)

## 
## Call:
## lm(formula = medv ~ ., data = bost)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.961  -4.985  -1.776   2.140  31.868 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 32.73979    2.24705  14.570  < 2e-16 ***
## age         -0.14280    0.01981  -7.208 2.09e-12 ***
## dis         -0.24616    0.26514  -0.928    0.354    
## chas1        7.51297    1.46466   5.130 4.16e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.317 on 502 degrees of freedom
## Multiple R-squared:  0.187,  Adjusted R-squared:  0.1822 
## F-statistic:  38.5 on 3 and 502 DF,  p-value: < 2.2e-16

Utilizando el p-valor de cada regresor, todos los regresores son significativos excepto la variable dis . Dado que la variable chas es categórica con dos niveles, el coeficiente de chas1 es la diferencia entre el precio madiano de las casas en los límites del río charles y el de las demás, por lo que el precio medio de las primeras es mayor alrededor de 7.513.

6. Regresión lineal bayesiana

Considerando \(\boldsymbol{\theta}=\left(\theta_{0}, \ldots, \theta_{k}\right)^{\top}\) como el vector de parámetros a estimar, \(\mathbf{y}=\left(y_{1}, \ldots, y_{n}\right)^{\top}\) un vector \(n \times 1\) donde es una realización de la variable aleatoria cuya distribución es \(p\left(y_{i} \mid \boldsymbol{\theta}\right)\). La función de verosimilitud es dado por \(y_{i}\):

\[\begin{equation} \mathcal{L}(\boldsymbol{\theta} \mid \mathbf{y})=\prod_{i=1}^{n} p\left(y_{i} \mid \boldsymbol{\theta}\right). \end{equation}\]

Toda la información de las observaciones \(y_{i}\) sobre \(\boldsymbol{\theta}\) se incluye en la ecuación anterior. La dificultad para estimar \(\boldsymbol{\theta}\) se vuelve un problema de optimización de maximizar la función de probabilidad (o su logaritmo). En contraste, bajo la metodología bayesiana la estimación de \(\boldsymbol{\theta}\) viene dada por la distribución posterior conjunta, que está definida por la Teorema de Bayes:

\[\begin{equation} p(\boldsymbol{\theta} \mid \mathbf{y})=\frac{\mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta})}{\int_{\Theta} \mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta}) d \boldsymbol{\theta}} \end{equation}\]

donde \(\Theta\) representa el espacio paramétrico de \(\boldsymbol{\theta}\) y \(p(\boldsymbol{\theta})\) es la distribución a priori. La ecuación anterior se puede escribir como:

\[\begin{equation} p(\boldsymbol{\theta} \mid \mathbf{y}) \propto \mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta}) \end{equation}\]

donde \(\int_{\Theta} \mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta}) d \boldsymbol{\theta}\) es la distribución marginal de \(\mathbf{y}\) que no depende de \(\boldsymbol{\theta}.\) La distribución a posteriori \(p(\boldsymbol{\theta} \mid \mathbf{y})\) proporciona toda la información que uno puede tener sobre \(\boldsymbol{\theta}\). Por ejemplo, es posible evaluar \(p(\boldsymbol{\theta} \mid \mathbf{y})\) y su media, mediana, varianza y algunas otras cantidades como cuantiles para tener estimaciones puntuales e intervalos. Además, la distribución a posteriori con frecuencia no tiene forma cerrada, por lo que va depender de los métodos computacionales para poder simular la distribución.

Para ajustar una regresión lineal bayesiana usaremos la función stan_glm de la librería rstanarm. Esta función, como la función lm anterior, requiere proporcionar la fórmula y los datos que se utilizarán, y dejar todos los siguientes argumentos con sus valores predeterminados. A continuación se procede a explicar los principales argumento de la función stan_glm.

family: por defecto esta función usa la distribución gaussiana como lo hacemos con la función clásica glm.
prior: La distribución a priori para los coeficientes de regresión. Por defecto, se usa la distribución normal. Hay un subconjunto de funciones que se utilizan para el a priori proporcionado por la libreria rstanarm como, la familia t-student, familia laplace, etc. Para obtener la lista completa con todos los detalles, ejecute este comando ?priors. Si queremos una distribución a priori uniforme, lo configuramos en NULL.
prior_intercept: priori para el intercepto, puede ser normal, t-student o cauchy. Si queremos una distribución a priori uniforme, lo configuramos en NULL.
prior_aux: priori para parámetros auxiliares como la desviación estándar del error para la familia de gaussiana.
algorithm: El método de estimación a utilizar. El valor predeterminado es el Markov chain Monte Carlo (MCMC).
QR: escalar lógico por defecto a FALSE, pero si es TRUE aplica una descomposición QR escalada a la matriz de diseño.
iter: es el número de iteraciones si se utiliza el método MCMC, el valor predeterminado es 2000.
chains: el número de cadenas de Markov, el valor predeterminado es 4.
warmup: también conocido como burnin, el número de iteraciones que se utilizan para la adaptación, y no debe utilizarse para inferencias. Por defecto es la mitad de las iteraciones.

library(rstanarm)
model_bayes <- stan_glm(medv ~., data=bost, seed=111)

## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.102 seconds (Warm-up)
## Chain 1:                0.098 seconds (Sampling)
## Chain 1:                0.2 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.058 seconds (Warm-up)
## Chain 2:                0.115 seconds (Sampling)
## Chain 2:                0.173 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.058 seconds (Warm-up)
## Chain 3:                0.092 seconds (Sampling)
## Chain 3:                0.15 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.065 seconds (Warm-up)
## Chain 4:                0.089 seconds (Sampling)
## Chain 4:                0.154 seconds (Total)
## Chain 4:

print(model_bayes, digits = 3)

## stan_glm
##  family:       gaussian [identity]
##  formula:      medv ~ .
##  observations: 506
##  predictors:   4
## ------
##             Median MAD_SD
## (Intercept) 32.834  2.285
## age         -0.143  0.020
## dis         -0.258  0.257
## chas1        7.543  1.432
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 8.324  0.260 
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

La estimación de la mediana es la calculada a partir de la simulación MCMC y MAD_SD es la desviación absoluta mediana calculada a partir de la misma simulación. Para comprender bien cómo obtener estos resultados, tracemos la simulación MCMC de cada predictor usando bayesplot.

library(bayesplot)
mcmc_dens(model_bayes, pars = c("age"))+
  vline_at(-0.143, col="red")

Como se puede observar, la estimación puntual de la variable age cae en la mediana de esta distribución (línea roja). Lo mismo ocurre con las variables dis y shas.

mcmc_dens(model_bayes, pars=c("chas1"))+
  vline_at(7.496, col="red")

mcmc_dens(model_bayes, pars=c("dis"))+
  vline_at(-0.244, col="red")

Ahora, ¿cómo podemos evaluar los parámetros del modelo?. La respuesta es analizando las distribuciones a posteriori y utilizando algunas estadísticas específicas. Para obtener estas estadísticas usaremos la función describe_posterior de la librería bayestestR.

library(bayestestR)
describe_posterior(model_bayes)

## # Description of Posterior Distributions
## 
## Parameter   | Median |           89% CI |      pd |        89% ROPE | % in ROPE |  Rhat |      ESS
## --------------------------------------------------------------------------------------------------
## (Intercept) | 32.834 | [29.218, 36.295] | 100.00% | [-0.920, 0.920] |         0 | 1.002 | 2029.279
## age         | -0.143 | [-0.175, -0.112] | 100.00% | [-0.920, 0.920] |       100 | 1.001 | 2052.155
## dis         | -0.258 | [-0.667,  0.179] |  81.92% | [-0.920, 0.920] |       100 | 1.002 | 2115.192
## chas1       |  7.543 | [ 5.159,  9.813] | 100.00% | [-0.920, 0.920] |         0 | 1.000 | 3744.403

7. Criterios para el análisis bayesiano

CI:Intervalo Credibilidad, se utiliza para cuantificar la incertidumbre sobre los coeficientes de regresión. Hay métodos de arrastre para calcular el CI, el intervalo de densidad más alto HDI que es el predeterminado y el intervalo inferior ETI. con un 89% de probabilidad (dados los datos) de que un coeficiente se encuentre por encima del valor CI_low y por debajo del valor CI_high. Esta interpretación probabilística estricta es completamente diferente del intervalo de confianza utilizado en la regresión lineal clásica donde el coeficiente cae dentro de este intervalo de confianza (si elegimos el 95% de confianza) 95 veces si repetimos el estudio 100 veces.
PD: Probabilidad de Dirección, que es la probabilidad de que el efecto vaya en dirección positiva o negativa, y se considera como el mejor equivalente del p-valor.
ROPE_CI: Región de Equivalencia Práctica, dado que el método de Bayes trata con probabilidades verdaderas, no tiene sentido calcular la probabilidad de obtener el efecto igual a cero (la hipótesis nula) como un punto (probabilidad de que un punto en intervalos continuos sea igual a cero). Por lo tanto, definimos en cambio un pequeño rango alrededor de cero que se puede considerar prácticamente lo mismo que ningún efecto (cero), por lo que este rango se llama ROPE. Por defecto (según Cohen, 1988) el rope de los coeficientes estandarizados es [-0.1,0.1].
Rhat: Factor de Reducción de Escala (\(\hat R\)), se calcula para cada cantidad escalar de interés, como la desviación estándar de esa cantidad de todas las cadenas incluidas juntas, dividida por la raíz cuadrada media de las desviaciones estándar dentro de la cadena separadas. Cuando este valor se acerca a 1 no tenemos ningún problema de convergencia mediante MCMC.
ESS: tamaño de muestra efectivo, captura cuántos extractos independientes contienen la misma cantidad de información que la muestra dependiente obtenida por el algoritmo MCMC, cuanto mayor sea el ESS, mejor. El umbral utilizado en la práctica es 400.

Alternativamente, podemos obtener las estimaciones de coeficientes (que son las medianas por defecto) por separado utilizando el paquete insight.

library(insight)
post <- get_parameters(model_bayes)
print(purrr::map_dbl(post,median),digits = 3)

## (Intercept)         age         dis       chas1 
##      32.834      -0.143      -0.258       7.543

También podemos calcular el Máximo a Posteriori (MAP) y la media de la siguiente manera

print(purrr::map_dbl(post, map_estimate),digits = 3)

## (Intercept)         age         dis       chas1 
##      33.025      -0.145      -0.295       7.573

print(purrr::map_dbl(post, mean),digits = 3)

## (Intercept)         age         dis       chas1 
##      32.761      -0.143      -0.248       7.523

Como vemos, los valores están más cerca entre sí debido a la misma normalidad de la distribución de las distribuciones a posteriori, donde todas las estadísticas centrales (media, mediana, moda) están más cerca entre sí. Usando la siguiente gráfica para visualizar el coeficiente de la variable age observamos que:

mcmc_dens(model_bayes, pars=c("age"))+
  vline_at(median(post$age), col="red")+
  vline_at(mean(post$age), col="yellow")+
  vline_at(map_estimate(post$age), col="green")

Como era de esperar, están aproximadamente uno encima del otro.

8. ¿Como probar la significancia de los coeficientes?

Como hacemos con la regresión lineal clásica (frecuentista), podemos probar la significancia de los coeficientes de regresión bayesiana comprobando si el intervalo de credibilidad correspondiente contiene a cero o no, si no, entonces este coeficiente es significativo. Regresemos a nuestro modelo y verifiquemos la importancia de cada coeficiente (usando credibilidad basada en el HDI predeterminado).

hdi(model_bayes)

## # Highest Density Interval
## 
## Parameter   |        89% HDI
## ----------------------------
## (Intercept) | [29.22, 36.29]
## age         | [-0.18, -0.11]
## dis         | [-0.67,  0.18]
## chas1       | [ 5.16,  9.81]

eti(model_bayes)

## # Equal-Tailed Interval
## 
## Parameter   |        89% ETI
## ----------------------------
## (Intercept) | [29.20, 36.28]
## age         | [-0.17, -0.11]
## dis         | [-0.67,  0.18]
## chas1       | [ 5.17,  9.83]

Usando ambos métodos, el único coeficiente no significativo es la variable dis, que está coherente con lo obtenido al usar la regresión lineal clásica.

Nota: este resultado similar entre la regresión lineal clásica y bayesiana puede deberse al supuesto de normalidad para el primero que está bien satisfecho, lo que da resultados satisfactorios. Sin embargo, en el mundo real es menos frecuente estar seguro del supuesto de normalidad que puede dar conclusiones contradictorias entre los dos enfoques.

Otra forma de probar la significancia al verificar la parte del intervalo de credibilidad que cae dentro del intervalo ROPE. Podemos obtener esto llamando a la función rope del paquete bayestestR.

rope(post$age)

## # Proportion of samples inside the ROPE [-0.10, 0.10]:
## 
## inside ROPE
## -----------
## 0.00 %

Para la variable age, casi todo el intervalo (IDH) está fuera del rango de ROPE, lo que significa que el coeficiente es muy significativo.

rope(post$chas1)

## # Proportion of samples inside the ROPE [-0.10, 0.10]:
## 
## inside ROPE
## -----------
## 0.00 %

rope(post$`(Intercept)`)

## # Proportion of samples inside the ROPE [-0.10, 0.10]:
## 
## inside ROPE
## -----------
## 0.00 %

Lo mismo es cierto para la variable chas y el intercepto.

rope(post$dis)

## # Proportion of samples inside the ROPE [-0.10, 0.10]:
## 
## inside ROPE
## -----------
## 20.02 %

En contraste, casi la quinta parte del intervalo de credibilidad de la variable dis está dentro del intervalo ROPE. En otras palabras, la probabilidad de que este coeficiente sea cero es del 20,02%.

rope_range(model_bayes)

## [1] -0.9197104  0.9197104

9. ¿Existe el P-valor en el enfoque Bayesiano?

A veces solo nos interesa comprobar la dirección del coeficiente (positivo o negativo). Este es el papel del estadístico PD en la tabla anterior, un valor alto significa que el efecto asociado se concentra en el mismo lado que la mediana. Para nuestro modelo, dado que PD es igual a 1, casi todas las distribuciones a posteriori de las variables age, chas1 y el intercepto están en el mismo lado (si la mediana es negativa, todos los demás valores son negativos). Sin embargo, cabe señalar que esta estadística no evalúa la importancia del efecto. Algo más importante de mencionar es que existe una fuerte relación entre esta probabilidad y el p-valor que es aproximado de la siguiente manera: \(p-value=1-pd\). comprobemos esto con nuestras variables en estudio.

library(broom)
df1 <-dplyr::select(tidy(model_freq), c(term,p.value))
df1$p.value <- round(df1$p.value, digits = 3)
df2 <- 1- purrr::map_dbl(post, p_direction)
df <- cbind(df1,df2)
df

##                    term p.value     df2
## (Intercept) (Intercept)   0.000 0.00000
## age                 age   0.000 0.00000
## dis                 dis   0.354 0.18075
## chas1             chas1   0.000 0.00000

10. Conclusión

En los último años, especialmente en algunos campos como la medicina y la psicología, se están volcando hacia el análisis bayesiano, ya que casi todo puede interpretarse directamente de manera probabilística. Sin embargo, el análisis bayesiano también tiene algunos inconvenientes, como la forma subjetiva de definir las distribuciones a prioris (que juegan un papel importante para calcular la distribución a posteriori), o para problemas que no tienen un priori conjugado, no siempre el algoritmo MCMC converge fácilmente a los valores correctos (especialmente con datos complejos).