17.9: ANOVA Bayesiano
- Page ID
- 151340
\( \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}\)Como puedes ver, el paquete BayesFactor
es bastante flexible, y puede hacer versiones bayesianas de casi todo en este libro. De hecho, puede hacer algunas otras cosas pulcras que no he cubierto en absoluto en el libro. No obstante, tengo que detenerme en alguna parte, así que solo hay otro tema que quiero cubrir: ANOVA bayesiano.
actualización rápida
Al igual que con los otros ejemplos, creo que es útil comenzar con un recordatorio de cómo discutí ANOVA anteriormente en el libro. Primero, recordemos cuáles fueron los datos. El ejemplo que utilicé originalmente es el marco de datos clin.trial
, que se ve así
load("./rbook-master/data/clinicaltrial.Rdata")
head(clin.trial)
## drug therapy mood.gain
## 1 placebo no.therapy 0.5
## 2 placebo no.therapy 0.3
## 3 placebo no.therapy 0.1
## 4 anxifree no.therapy 0.6
## 5 anxifree no.therapy 0.4
## 6 anxifree no.therapy 0.2
Para ejecutar nuestro análisis ortodoxo en capítulos anteriores utilizamos la función aov ()
para hacer todo el trabajo pesado. En el Capítulo 16 recomendé usar la función Anova ()
del paquete de autos
para producir la tabla ANOVA, ya que usa pruebas de Tipo II por defecto. Si has olvidado lo que son las “pruebas Tipo II”, podría ser una buena idea volver a leer la Sección 16.10, porque volverá a ser relevante en un momento. En cualquier caso, así es como se veía nuestro análisis:
library(car)
## Loading required package: carData
model <- aov( mood.gain ~ drug * therapy, data = clin.trial )
Anova(model)
## Anova Table (Type II tests)
##
## Response: mood.gain
## Sum Sq Df F value Pr(>F)
## drug 3.4533 2 31.7143 1.621e-05 ***
## therapy 0.4672 1 8.5816 0.01262 *
## drug:therapy 0.2711 2 2.4898 0.12460
## Residuals 0.6533 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Eso nos muestra claramente evidencias de un efecto principal del fármaco
a p<.001, un efecto de la terapia
a p<.05 y sin interacción.
Versión bayesiana
¿Cómo hacemos lo mismo usando métodos bayesianos? El paquete BayesFactor
contiene una función llamada AnovaBF ()
que hace esto por usted. Utiliza una fórmula
bastante estándar y una estructura de datos
, por lo que el comando debería parecer muy familiar. Al igual que hicimos con la regresión, será útil guardar la salida en una variable:
models <- anovaBF(
formula = mood.gain ~ drug * therapy,
data = clin.trial
)
El resultado es bastante diferente al ANOVA tradicional, pero no es tan malo una vez que entiendes lo que buscas. Echemos un vistazo:
models
Esto se ve muy similar a la salida que obtuvimos de la función regressionBF ()
, y con buena razón. Recuerden lo que dije en la Sección 16.6: bajo el capó, el ANOVA no es diferente a la regresión, y ambos son solo ejemplos diferentes de un modelo lineal. Debido a esto, el anovaBF ()
reporta la salida de la misma manera. Por ejemplo, si queremos identificar el mejor modelo podríamos usar los mismos comandos que usamos en la última sección. Una variante que me parece bastante útil es esta:
models/max(models)
## Bayes factor analysis
## --------------
## [1] drug : 0.3521042 ±0.94%
## [2] therapy : 0.001047568 ±0.94%
## [3] drug + therapy : 1 ±0%
## [4] drug + therapy + drug:therapy : 0.978514 ±1.29%
##
## Against denominator:
## mood.gain ~ drug + therapy
## ---
## Bayes factor type: BFlinearModel, JZS
Al “dividir” la salida de los modelos
por el mejor modelo (es decir, max (modelos)
), lo que R está haciendo es usar el mejor modelo (que en este caso es drogas + terapia
) como denominador, lo que te da una idea bastante clara de lo cerca que están los competidores. Por ejemplo, el modelo que contiene el término de interacción es casi tan bueno como el modelo sin la interacción, ya que el factor Bayes es de 0.98. Es decir, los datos no indican claramente si existe o no una interacción.
Construcción de pruebas bayesianas tipo II
Bien, eso está muy bien, podrías estar pensando, pero ¿qué informo como alternativa al valor p? En la tabla ANOVA clásica, se obtiene un único valor p por cada predictor del modelo, así se puede hablar sobre la significancia de cada efecto. ¿Cuál es el análogo bayesiano de esto?
Es una buena pregunta, pero la respuesta es complicada. Recuerda lo que dije en la Sección 16.10 sobre el hecho de que ANOVA sea complicado. Incluso en la versión clásica del ANOVA hay varias “cosas” diferentes a las que el ANOVA podría corresponder. Específicamente, comenté cómo se obtienen diferentes valores p dependiendo de si usa pruebas de Tipo I, Pruebas de Tipo II o pruebas de Tipo III. Para averiguar qué factor Bayes es análogo al “el” valor p en un ANOVA clásico, necesitas averiguar para qué versión de ANOVA quieres un análogo. Para los efectos de esta sección, asumiré que quiere pruebas de Tipo II, porque esas son las que creo que son más sensatas en general. Como ya comenté en la Sección 16.10, las pruebas Tipo II para un ANOVA bidireccional son razonablemente sencillas, pero si has olvidado esa sección no sería mala idea volver a leerla antes de continuar.
Asumiendo que has tenido un repaso en las pruebas Tipo II, echemos un vistazo a cómo sacarlas de la tabla de factores Bayes. Supongamos que queremos probar el efecto principal de la droga
. La hipótesis nula para esta prueba corresponde a un modelo que incluye un efecto de terapia
, pero ningún efecto de fármaco
. La hipótesis alternativa es el modelo que incluye ambos. Es decir, lo que queremos es el factor Bayes correspondiente a esta comparación:
knitr::kable(tibble::tribble(
~V1, ~V2,
"Null model:", "`mood.gain ~ therapy`",
"Alternative model:", "`mood.gain ~ therapy + drug`"
), col.names = c("", ""))
Modelo nulo: | estado de ánimo.ganancia ~ terapia |
Modelo alternativo: | estado de ánimo.ganancia ~ terapia + fármaco |
Como sucede, podemos leer la respuesta a esto directamente de la tabla porque corresponde a una comparación entre el modelo en la línea 2 de la tabla y el modelo en la línea 3: el factor Bayes en este caso representa evidencia para el nulo de 0.001 a 1. O, más amablemente, las probabilidades son de aproximadamente 1000 a 1 contra el nulo.
El efecto principal de la terapia
se puede calcular de la misma manera. En este caso, el modelo nulo es el que contiene sólo un efecto de fármaco, y la alternativa es el modelo que contiene ambos. Por lo que la comparación relevante es entre las líneas 2 y 1 de la tabla. Las probabilidades a favor del nulo aquí son solo de 0.35 a 1. Nuevamente, me parece útil enmarcar las cosas al revés, así que me referiría a esto como evidencia de aproximadamente 3 a 1 a favor de un efecto de la terapia
.
Finalmente, para probar un efecto de interacción, el modelo nulo aquí es uno que contiene ambos efectos principales pero ninguna interacción. El modelo alternativo agrega la interacción. Es decir:
knitr::kable(tibble::tribble(
~V1, ~V2,
"Null model:", "`mood.gain ~ drug + therapy`",
"Alternative model:", "`mood.gain ~ drug + therapy + drug:therapy`"
), col.names = c("", ""))
Modelo nulo: | estado de ánimo.ganancia ~ fármaco + terapia |
Modelo alternativo: | estado de ánimo.ganancia ~ fármaco + terapia + medicamento:terapia |
Si miramos esos dos modelos arriba en la tabla, vemos que esta comparación es entre los modelos de las líneas 3 y 4 de la tabla. Las probabilidades de 0.98 a 1 implican que estos dos modelos están bastante igualados.
Podrías estar pensando que todo esto es bastante laborioso, y admito que es verdad. En algún momento podría considerar agregar una función al paquete lsr
que automatizaría este proceso y construiría algo así como una “tabla ANOVA Bayesiana Tipo II” a partir de la salida de la función anovaBF ()
. Sin embargo, todavía no he tenido tiempo de hacer esto, ni me he tomado una decisión sobre si realmente es una buena idea hacer esto. Mientras tanto, pensé que debería mostrarte el truco de cómo hago esto en la práctica. El comando que utilizo cuando quiero agarrar los factores Bayes correctos para un ANOVA Tipo II es este:
max(models)/models
## denominator
## numerator drug therapy drug + therapy
## drug + therapy 2.840068 954.5918 1
## denominator
## numerator drug + therapy + drug:therapy
## drug + therapy 1.021958
El resultado no es tan bonito como el último, pero lo bueno es que puedes leer todo lo que necesites. El mejor modelo es la farmacoterapia
, por lo que todos los demás modelos se están comparando con eso. ¿Cuál es el factor Bayes para el principal efecto de la droga
? La hipótesis nula relevante es la que contiene sólo terapia
, y el factor Bayes en cuestión es 954:1. El efecto principal de la terapia
es más débil, y la evidencia aquí es de solo 2. 8:1. Finalmente, la evidencia contra una interacción es muy débil, a 1. 01:1.
Leer los resultados de esta tabla es algo contrario a la intuición, porque hay que leer las respuestas de la parte “incorrecta” de la tabla. Por ejemplo, la evidencia de un efecto de la droga
se puede leer de la columna etiquetada como terapia
, lo cual es bastante extraño. Para ser justos con los autores del paquete, no creo que alguna vez pensaron que la función AnovaBF ()
se usara de esta manera. Mi entendimiento 273 es que su punto de vista es simplemente que deberías encontrar el mejor modelo y reportar ese modelo: no hay ninguna razón inherente por la que un ANOVA bayesiano deba tratar de seguir exactamente el mismo diseño que un ANOVA ortodoxo. 274
En cualquier caso, si sabes lo que estás buscando, puedes mirar esta tabla y luego reportar los resultados del análisis bayesiano de una manera bastante análoga a cómo reportarías un ANOVA tipo II regular. Como mencioné anteriormente, todavía no hay convención sobre cómo hacerlo, pero suelo ir por algo como esto:
Un ANOVA Bayesiano Tipo II encontró evidencia de los principales efectos del fármaco (factor Bayes: 954:1) y la terapia (factor Bayes: 3:1), pero no hay evidencia clara a favor o en contra de una interacción (factor Bayes: 1:1).