Saltar al contenido principal

# 17.8: Regresión Bayesiana

$$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$

$$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$

$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$

( \newcommand{\kernel}{\mathrm{null}\,}\) $$\newcommand{\range}{\mathrm{range}\,}$$

$$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$

$$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}[1]{\| #1 \|}$$

$$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$$

$$\newcommand{\Span}{\mathrm{span}}$$

$$\newcommand{\id}{\mathrm{id}}$$

$$\newcommand{\Span}{\mathrm{span}}$$

$$\newcommand{\kernel}{\mathrm{null}\,}$$

$$\newcommand{\range}{\mathrm{range}\,}$$

$$\newcommand{\RealPart}{\mathrm{Re}}$$

$$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$

$$\newcommand{\Argument}{\mathrm{Arg}}$$

$$\newcommand{\norm}[1]{\| #1 \|}$$

$$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$$

$$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\AA}{\unicode[.8,0]{x212B}}$$

$$\newcommand{\vectorA}[1]{\vec{#1}} % arrow$$

$$\newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow$$

$$\newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$

$$\newcommand{\vectorC}[1]{\textbf{#1}}$$

$$\newcommand{\vectorD}[1]{\overrightarrow{#1}}$$

$$\newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}}$$

$$\newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}}$$

$$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$

$$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$

$$\newcommand{\avec}{\mathbf a}$$ $$\newcommand{\bvec}{\mathbf b}$$ $$\newcommand{\cvec}{\mathbf c}$$ $$\newcommand{\dvec}{\mathbf d}$$ $$\newcommand{\dtil}{\widetilde{\mathbf d}}$$ $$\newcommand{\evec}{\mathbf e}$$ $$\newcommand{\fvec}{\mathbf f}$$ $$\newcommand{\nvec}{\mathbf n}$$ $$\newcommand{\pvec}{\mathbf p}$$ $$\newcommand{\qvec}{\mathbf q}$$ $$\newcommand{\svec}{\mathbf s}$$ $$\newcommand{\tvec}{\mathbf t}$$ $$\newcommand{\uvec}{\mathbf u}$$ $$\newcommand{\vvec}{\mathbf v}$$ $$\newcommand{\wvec}{\mathbf w}$$ $$\newcommand{\xvec}{\mathbf x}$$ $$\newcommand{\yvec}{\mathbf y}$$ $$\newcommand{\zvec}{\mathbf z}$$ $$\newcommand{\rvec}{\mathbf r}$$ $$\newcommand{\mvec}{\mathbf m}$$ $$\newcommand{\zerovec}{\mathbf 0}$$ $$\newcommand{\onevec}{\mathbf 1}$$ $$\newcommand{\real}{\mathbb R}$$ $$\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}$$ $$\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}$$ $$\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}$$ $$\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}$$ $$\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$$ $$\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$$ $$\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$$ $$\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$$ $$\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}$$ $$\newcommand{\laspan}[1]{\text{Span}\{#1\}}$$ $$\newcommand{\bcal}{\cal B}$$ $$\newcommand{\ccal}{\cal C}$$ $$\newcommand{\scal}{\cal S}$$ $$\newcommand{\wcal}{\cal W}$$ $$\newcommand{\ecal}{\cal E}$$ $$\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}$$ $$\newcommand{\gray}[1]{\color{gray}{#1}}$$ $$\newcommand{\lgray}[1]{\color{lightgray}{#1}}$$ $$\newcommand{\rank}{\operatorname{rank}}$$ $$\newcommand{\row}{\text{Row}}$$ $$\newcommand{\col}{\text{Col}}$$ $$\renewcommand{\row}{\text{Row}}$$ $$\newcommand{\nul}{\text{Nul}}$$ $$\newcommand{\var}{\text{Var}}$$ $$\newcommand{\corr}{\text{corr}}$$ $$\newcommand{\len}[1]{\left|#1\right|}$$ $$\newcommand{\bbar}{\overline{\bvec}}$$ $$\newcommand{\bhat}{\widehat{\bvec}}$$ $$\newcommand{\bperp}{\bvec^\perp}$$ $$\newcommand{\xhat}{\widehat{\xvec}}$$ $$\newcommand{\vhat}{\widehat{\vvec}}$$ $$\newcommand{\uhat}{\widehat{\uvec}}$$ $$\newcommand{\what}{\widehat{\wvec}}$$ $$\newcommand{\Sighat}{\widehat{\Sigma}}$$ $$\newcommand{\lt}{<}$$ $$\newcommand{\gt}{>}$$ $$\newcommand{\amp}{&}$$ $$\definecolor{fillinmathshade}{gray}{0.9}$$

Bien, entonces ahora hemos visto equivalentes bayesianos a pruebas ortodoxas de chi-cuadrado y pruebas t. ¿Cuál es el siguiente? Si tuviera que seguir la misma progresión que usé al desarrollar las pruebas ortodoxas esperarías ver ANOVA a continuación, pero creo que es un poco más claro si empezamos con la regresión.

## actualización rápida

En el capítulo 15 utilicé los datos de paternidad para ilustrar las ideas básicas detrás de la regresión. Para recordarte cómo es ese conjunto de datos, aquí tienes las primeras seis observaciones:

load("./rbook-master/data/parenthood.Rdata")
head(parenthood)
##   dan.sleep baby.sleep dan.grump day
## 1      7.59      10.18        56   1
## 2      7.91      11.66        60   2
## 3      5.14       7.92        82   3
## 4      7.71       9.61        55   4
## 5      6.68       9.75        67   5
## 6      5.99       5.04        72   6

De vuelta en el Capítulo 15 propuse una teoría en la que mi maldad (dan.gruñón) en un día cualquiera está relacionado con la cantidad de sueño que tuve la noche anterior (dan.sleep), y posiblemente con la cantidad de sueño que nuestro bebé consiguió (bebé.dormir), aunque probablemente no al día en que tomamos la medida. Se probó esto usando un modelo de regresión. Para estimar el modelo de regresión se utilizó la función lm (), así:

model <- lm(
formula = dan.grump ~ dan.sleep + day + baby.sleep,
data = parenthood
)

Las pruebas de hipótesis para cada uno de los términos del modelo de regresión se extrajeron utilizando la función summary () como se muestra a continuación:

summary(model)
##
## Call:
## lm(formula = dan.grump ~ dan.sleep + day + baby.sleep, data = parenthood)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -10.906  -2.284  -0.295   2.652  11.880
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) 126.278707   3.242492  38.945   <2e-16 ***
## dan.sleep    -8.969319   0.560007 -16.016   <2e-16 ***
## day          -0.004403   0.015262  -0.288    0.774
## baby.sleep    0.015747   0.272955   0.058    0.954
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.375 on 96 degrees of freedom
## Multiple R-squared:  0.8163, Adjusted R-squared:  0.8105
## F-statistic: 142.2 on 3 and 96 DF,  p-value: < 2.2e-16

Al interpretar los resultados, cada fila de esta tabla corresponde a uno de los posibles predictores. El término (Intercept) no suele ser interesante, aunque es muy significativo. Lo importante para nuestros propósitos es el hecho de que dan.sleep es significativo a p<.001 y ninguna de las otras variables lo son.

## Versión bayesiana

Bien, entonces, ¿cómo hacemos lo mismo usando el paquete BayesFactor? La forma más fácil es usar la función regressionBF () en lugar de lm (). Como antes, usamos la fórmula para indicar cómo se ve el modelo de regresión completo, y el argumento data para especificar el marco de datos. Entonces el comando es:

regressionBF(
formula = dan.grump ~ dan.sleep + day + baby.sleep,
data = parenthood
)
## Bayes factor analysis
## --------------
## [1] dan.sleep                    : 1.622545e+34 ±0.01%
## [2] day                          : 0.2724027    ±0%
## [3] baby.sleep                   : 10018411     ±0%
## [4] dan.sleep + day              : 1.016576e+33 ±0%
## [5] dan.sleep + baby.sleep       : 9.77022e+32  ±0%
## [6] day + baby.sleep             : 2340755      ±0%
## [7] dan.sleep + day + baby.sleep : 7.835625e+31 ±0%
##
## Against denominator:
##   Intercept only
## ---
## Bayes factor type: BFlinearModel, JZS

Entonces eso es bastante sencillo: es exactamente lo que hemos estado haciendo a lo largo del libro. La salida, sin embargo, es un poco diferente de lo que obtienes de lm (). El formato de esto es bastante familiar. En la parte inferior tenemos alguna basura tecnológica, y en la parte superior tenemos alguna información sobre los factores Bayes. Lo nuevo es el hecho de que parece que tenemos muchos factores Bayes aquí. ¿De qué se trata todo esto?

El truco para entender esta salida es reconocer que si estamos interesados en averiguar cuáles de las 3 variables predictoras están relacionadas con dan.grump, en realidad hay 8 posibles modelos de regresión que podrían considerarse. Una posibilidad es el modelo de intercepción única, en el que ninguna de las tres variables tiene efecto. En el otro extremo del espectro se encuentra el modelo completo en el que importan las tres variables. Entonces, lo que hace regressionBF () es tratar el modelo de intercepción únicamente como la hipótesis nula, e imprimir los factores Bayes para todos los demás modelos cuando se comparan con ese nulo. Por ejemplo, si miramos la línea 4 de la tabla, vemos que la evidencia es de aproximadamente 1033 a 1 a favor de la afirmación de que un modelo que incluye tanto dan.sleep como día es mejor que el modelo de intercepción única. O si miramos la línea 1, podemos ver que las probabilidades son de aproximadamente 1.6×1034 de que un modelo que contenga la variable dan.sleep (pero no otras) sea mejor que el modelo intercept only.

## Encontrar el mejor modelo

En la práctica, esto no es muy útil. En la mayoría de las situaciones el modelo de intercepción única es uno que realmente no te importa en absoluto. Lo que me parece útil es comenzar trabajando en qué modelo es el mejor, y luego ver qué tan bien se comparan todas las alternativas con él. Así es como lo haces. En este caso, es bastante fácil ver que el mejor modelo es en realidad el que contiene dan.sleep solamente (línea 1), porque tiene el factor Bayes más grande. Sin embargo, si tienes muchos modelos posibles en la salida, es útil saber que puedes usar la función head () para elegir los mejores modelos. Primero, tenemos que regresar y guardar la información del factor Bayes en una variable:

models <- regressionBF(
formula = dan.grump ~ dan.sleep + day + baby.sleep,
data = parenthood
)

Digamos que quiero ver los tres mejores modelos. Para ello, utilizo la función head () especificando n=3, y esto es lo que obtengo como resultado:

head( models, n = 3)

## Bayes factor analysis
## --------------
## [1] dan.sleep              : 1.622545e+34 ±0.01%
## [2] dan.sleep + day        : 1.016576e+33 ±0%
## [3] dan.sleep + baby.sleep : 9.77022e+32  ±0%
##
## Against denominator:
##   Intercept only
## ---
## Bayes factor type: BFlinearModel, JZS

Esto nos está diciendo que el modelo en la línea 1 (es decir, dan.grump ~ dan.sleep) es el mejor. Eso es casi lo que estoy buscando, pero sigue comparando todos los modelos con el modelo solo de interceptación. Eso parece una tontería. Lo que me gustaría saber es qué tan grande es la diferencia entre el mejor modelo y los otros buenos modelos. Para eso, está este truco:

head( models/max(models), n = 3)
## Bayes factor analysis
## --------------
## [1] dan.sleep              : 1         ±0%
## [2] dan.sleep + day        : 0.0626532 ±0.01%
## [3] dan.sleep + baby.sleep : 0.0602154 ±0.01%
##
## Against denominator:
##   dan.grump ~ dan.sleep
## ---
## Bayes factor type: BFlinearModel, JZS

Observe el bit en la parte inferior mostrando que el “denominador” ha cambiado. Lo que eso significa es que los factores Bayes están comparando ahora cada uno de esos 3 modelos listados con el modelo dan.grump ~ dan.sleep. Obviamente, el factor Bayes en la primera línea es exactamente 1, ya que eso es solo comparar el mejor modelo consigo mismo. Más al grano, los otros dos factores Bayes son ambos menores que 1, lo que indica que todos son peores que ese modelo. Los factores Bayes de 0.06 a 1 implican que las probabilidades para el mejor modelo sobre el segundo mejor modelo son alrededor de 16:1. Esto se puede resolver por simple aritmética (es decir, 0.06/1≈16), pero la otra forma de hacerlo es comparando directamente los modelos. Para ver a lo que me refiero, aquí está la salida original:

models
## Bayes factor analysis
## --------------
## [1] dan.sleep                    : 1.622545e+34 ±0.01%
## [2] day                          : 0.2724027    ±0%
## [3] baby.sleep                   : 10018411     ±0%
## [4] dan.sleep + day              : 1.016576e+33 ±0%
## [5] dan.sleep + baby.sleep       : 9.77022e+32  ±0%
## [6] day + baby.sleep             : 2340755      ±0%
## [7] dan.sleep + day + baby.sleep : 7.835625e+31 ±0%
##
## Against denominator:
##   Intercept only
## ---
## Bayes factor type: BFlinearModel, JZS

El mejor modelo corresponde a la fila 1 de esta tabla, y el segundo mejor modelo corresponde a la fila 4. Todo lo que tienes que hacer para comparar estos dos modelos es esto:

models[1] / models[4]
## Bayes factor analysis
## --------------
## [1] dan.sleep : 15.96088 ±0.01%
##
## Against denominator:
##   dan.grump ~ dan.sleep + day
## ---
## Bayes factor type: BFlinearModel, JZS

Y ahí lo tienes. Encontraste el modelo de regresión con el factor Bayes más alto (es decir, dan.grump ~ dan.sleep), y sabes que la evidencia de ese modelo sobre la siguiente mejor alternativa (es decir, dan.grump ~ dan.sleep + día) es de aproximadamente 16:1.

## Extracción de factores Bayes para todos los términos incluidos

Bien, digamos que te has decidido por un modelo de regresión específico. ¿Qué factores Bayes debes reportar? En este ejemplo, voy a fingir que decidiste que dan.gruñón ~ dan.sleep + baby.sleep es el modelo que crees que es mejor. A veces es sensato hacer esto, incluso cuando no es el que tiene el factor Bayes más alto. Normalmente esto sucede porque tienes una razón teórica sustantiva para preferir un modelo sobre el otro. Sin embargo, en este caso lo estoy haciendo porque quiero usar un modelo con más de un predictor como mi ejemplo!

Habiendo averiguado qué modelo prefieres, puede ser muy útil llamar a la función regressionBF () y especificar WhichModels="top”. Utiliza tu modelo “preferido” como argumento de fórmula, y luego la salida te mostrará los factores Bayes que resultan cuando intentas soltar predictores de este modelo:

regressionBF(
formula = dan.grump ~ dan.sleep + baby.sleep,
data = parenthood,
whichModels = "top"
)
## Bayes factor top-down analysis
## --------------
## When effect is omitted from dan.sleep + baby.sleep , BF is...
## [1] Omit baby.sleep : 16.60705     ±0.01%
## [2] Omit dan.sleep  : 1.025403e-26 ±0.01%
##
## Against denominator:
##   dan.grump ~ dan.sleep + baby.sleep
## ---
## Bayes factor type: BFlinearModel, JZS

Bien, entonces ahora puedes ver los resultados un poco más claramente. El factor Bayes cuando intentas dejar caer el predictor dan.sleep es de aproximadamente 10−26, lo que es una evidencia muy fuerte de que no debes dejarlo caer. Por otro lado, el factor Bayes en realidad sube a 17 si dejas caer bebé.dormir, así que normalmente dirías que esa es una evidencia bastante sólida para dejar caer esa.

This page titled 17.8: Regresión Bayesiana is shared under a CC BY-SA 4.0 license and was authored, remixed, and/or curated by Danielle Navarro via source content that was edited to the style and standards of the LibreTexts platform.