17.6: Análisis Bayesiano de Tablas de Contingencia

Tiempo para cambiar de marcha. Hasta este punto he estado hablando sobre qué es la inferencia bayesiana y por qué podrías considerar usarla. Ahora quiero describir brevemente cómo hacer versiones bayesianas de diversas pruebas estadísticas. Las discusiones en las próximas secciones no son tan detalladas como me gustaría, pero espero que sean suficientes para ayudarte a empezar. Entonces comencemos.

El primer tipo de problema de inferencia estadística que discutí en este libro apareció en el Capítulo 12, en el que discutimos problemas categóricos de análisis de datos. En ese capítulo hablé de varios problemas estadísticos distintos que te podrían interesar, pero el que aparece con mayor frecuencia en la vida real es el análisis de las tablas de contingencia. En este tipo de situación de análisis de datos, tenemos una tabulación cruzada de una variable contra otra, y el objetivo es averiguar si existe alguna asociación entre estas variables. El conjunto de datos que utilicé para ilustrar este problema se encuentra en el archivo Chapek9.rData, y contiene una sola trama de datos chapek9

load("./rbook-master/data/chapek9.Rdata")
head(chapek9)

##   species choice
## 1   robot flower
## 2   human   data
## 3   human   data
## 4   human   data
## 5   robot   data
## 6   human flower

En este conjunto de datos, supuestamente muestreamos 180 seres y medimos dos cosas. Primero, verificamos si eran humanos o robots, según lo capturado por la variable especie. Segundo, les pedimos que nominaran si más preferían flores, cachorros o datos. Cuando producimos la tabulación cruzada, obtenemos esto como los resultados:

crosstab <- xtabs( ~ species + choice, chapek9 )
crosstab

##        choice
## species puppy flower data
##   robot    13     30   44
##   human    15     13   65

Sorprendentemente, los humanos parecían mostrar una preferencia mucho más fuerte por los datos que los robots. En ese momento especulamos que esto podría haber sido porque el interrogante era un robot grande que llevaba una pistola, y los humanos podrían haber tenido miedo.

texto ortodoxo

Solo para refrescar tu memoria, así es como analizamos estos datos en el Capítulo @refch:chisquare. Debido a que queremos determinar si existe alguna asociación entre especies y elección, utilizamos la función associationTest () en el paquete lsr para ejecutar una prueba de asociación de chi-cuadrado. Los resultados se veían así:

library(lsr)

## Warning: package 'lsr' was built under R version 3.5.2 

associationTest( ~species + choice, chapek9 )

## 
##      Chi-square test of categorical association
## 
## Variables:   species, choice 
## 
## Hypotheses: 
##    null:        variables are independent of one another
##    alternative: some contingency exists between variables
## 
## Observed contingency table:
##        choice
## species puppy flower data
##   robot    13     30   44
##   human    15     13   65
## 
## Expected contingency table under the null hypothesis:
##        choice
## species puppy flower data
##   robot  13.5   20.8 52.7
##   human  14.5   22.2 56.3
## 
## Test results: 
##    X-squared statistic:  10.722 
##    degrees of freedom:  2 
##    p-value:  0.005 
## 
## Other information: 
##    estimated effect size (Cramer's v):  0.244

Debido a que encontramos un pequeño valor de p (en este caso p<.01), concluimos que los datos son inconsistentes con la hipótesis nula de no asociación, y la rechazamos.

Prueba bayesiana

¿Cómo ejecutamos una prueba equivalente como bayesiano? Bueno, como cualquier otra cosa sangrienta en las estadísticas, hay muchas formas diferentes de hacerlo. No obstante, por el bien de la cordura de todos, a lo largo de este capítulo he decidido confiar en un paquete R para hacer el trabajo. Específicamente, voy a usar el paquete BayesFactor escrito por Jeff Rouder y Rich Morey, que a partir de este escrito está en la versión 0.9.10.

Para el análisis de tablas de contingencia, el paquete BayesFactor contiene una función llamada ContingencyTableBF (). Los datos que necesitas darle a esta función son la tabla de contingencia en sí (es decir, la variable de tabulación cruzada anterior), por lo que podría estar esperando usar un comando como este:

library( BayesFactor )           # ...because we have to load the package
contingencyTableBF( crosstab )   # ...because that makes sense, right?

Sin embargo, si intentas esto recibirás un mensaje de error. Esto se debe a que la función ContingencyTestBF () necesita otra información tuya: necesita saber qué plan de muestreo utilizó para ejecutar su experimento. Puede especificar el plan de muestreo utilizando el argumento SampleTipo de muestra. ¡Entonces probablemente debería decirte cuáles son tus opciones! La función ContingencyTableBF () distingue entre cuatro tipos diferentes de experimentos:

• Tamaño de muestra fijo. Supongamos que en nuestro ejemplo chapek9, nuestro experimento fue diseñado así: deliberadamente nos propusimos a probar a 180 personas, pero no intentamos controlar el número de humanos o robots, ni intentamos controlar las elecciones que tomaron. En este diseño, el número total de observaciones N es fijo, pero todo lo demás es aleatorio. Esto se conoce como muestreo “multinomial conjunto”, y si eso es lo que hiciste debes especificar sampleType = “JointMulti”. En el caso de los datos chapek9, eso es en realidad lo que tenía en mente cuando inventé el conjunto de datos.
• Totales fijos de fila (o columna). Un tipo diferente de diseño podría funcionar así. Decidimos con anticipación que queremos 180 personas, pero tratamos de ser un poco más sistemáticos al respecto. Específicamente, el experimentador lo constriñe para que obtengamos un número predeterminado de humanos y robots (por ejemplo, 90 de cada uno). En este diseño, ya sea los totales de fila o los totales de columna son fijos, pero no ambos. Esto se conoce como muestreo “multinomial independiente”, y si eso es lo que hiciste debes especificar sampleType = “indepMulti”.
• Tanto los totales de fila como de columna fijos. Otra posibilidad lógica es que diseñaste el experimento para que tanto los totales de fila como los totales de columna sean fijos. Esto no tiene ningún sentido en el ejemplo chapek9, pero hay otros deisgns que pueden funcionar de esta manera. Supongamos que te muestro una colección de 20 juguetes, para luego darles 10 pegatinas que dicen chico y otros 10 que dicen niña. Luego les doy 10 pegatinas azules y 10 pegatinas rosas. Entonces te pido que pongas las pegatinas en los 20 juguetes de tal manera que cada juguete tenga un color y cada juguete tenga un género. No importa cómo asignes las pegatinas, el número total de juguetes rosas y azules será de 10, al igual que el número de niños y niñas. En este diseño se fijan tanto las filas como las columnas de la tabla de contingencia. Esto se conoce como muestreo “hipergeométrico”, y si eso es lo que has hecho debes especificar sampleType = “hypergeom”.
• No se arregla nada. Por último, podría darse el caso de que no se arregle nada. Ni las columnas de fila, ni los totales de columna, ni tampoco el tamaño total de la muestra. Por ejemplo, en el escenario chapek9, supongamos que lo que había hecho es ejecutar el estudio por un período de tiempo fijo. Por casualidad, resultó que conseguí que 180 personas acudieran a estudiar, pero fácilmente podría haber sido otra cosa. Esto se conoce como muestreo de “Poisson”, y si eso es lo que has hecho debes especificar sampleType="Poisson”.

Bien, entonces ahora tenemos el conocimiento suficiente para hacer una prueba de hecho. Para los datos chapek9, insinué que diseñamos el estudio de tal manera que se fijara el tamaño total de la muestra N, por lo que deberíamos establecer sampleType = “JointMulti”. El comando que necesitamos es,

library( BayesFactor )

## Warning: package 'BayesFactor' was built under R version 3.5.2

## Loading required package: coda

## Warning: package 'coda' was built under R version 3.5.2

## Loading required package: Matrix

## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey ([e-mail]).
## 
## Type BFManual() to open the manual.
## ************ 

contingencyTableBF( crosstab, sampleType = "jointMulti" )

## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 15.92684 ±0%
## 
## Against denominator:
##   Null, independence, a = 1 
## ---
## Bayes factor type: BFcontingencyTable, joint multinomial

Al igual que con la mayoría de los comandos R, la salida inicialmente parece sospechosamente similar al galimatías absoluto. Afortunadamente, en realidad es bastante simple una vez que superas la impresión inicial. En primer lugar, tenga en cuenta que las cosas en la parte superior e inferior son pelusa irrelevante. Ya sabes que estás haciendo un análisis factorial de Bayes. Ya sabes que estás analizando una tabla de contingencia, y ya sabes que especificaste un plan de muestreo multinomial conjunto. Así que vamos a despojarnos de eso y echemos un vistazo a lo que queda:

[1] Non-indep. (a=1) : 15.92684 @plusorminus0%

Against denominator:
  Null, independence, a = 1 

Ignoremos también esos dos bits a=1, ya que son detalles técnicos que no necesitas conocer en esta etapa. ²⁶⁸ El resto de la salida es en realidad bastante sencillo. En la parte inferior, la salida define la hipótesis nula para ti: en este caso, la hipótesis nula es que no hay relación entre especies y elección. O, para decirlo de otra manera, la hipótesis nula es que estas dos variables son independientes. Ahora bien, si miras la línea por encima de ella, podrías (correctamente) adivinar que el no-INDEP. parte se refiere a la hipótesis alternativa. En este caso, la alternativa es que existe una relación entre las especies y la elección: es decir, no son independientes. Entonces lo único que queda en la salida es el bit que lee

15.92684 @plusorminus0%

La parte 15.9 es el factor Bayes, y te está diciendo que las probabilidades para la hipótesis alternativa contra el nulo son alrededor de 16:1. La parte ± 0% no es muy interesante: esencialmente, todo lo que te dice es que R ha calculado un factor Bayes exacto, por lo que la incertidumbre sobre el factor Bayes es del 0%. ²⁶⁹ En cualquier caso, los datos nos están diciendo que tenemos evidencia moderada para la hipótesis alternativa.

Redactar los resultados

Al escribir los resultados, mi experiencia ha sido que no hay tantas “reglas” sobre cómo “deberías” reportar las pruebas de hipótesis bayesianas. Eso podría cambiar en el futuro si los métodos bayesianos se vuelven estándar y algún grupo de trabajo comienza a escribir guías de estilo, pero mientras tanto sugeriría usar algún sentido común. Por ejemplo, evitaría escribir esto:

Una prueba bayesiana de asociación encontró un resultado significativo (BF=15.92)

En mi opinión, este escrito no está claro. Incluso asumiendo que ya has reportado las estadísticas descriptivas relevantes, hay una serie de cosas con las que no estoy contento. En primer lugar, el concepto de “significancia estadística” está bastante ligado a los valores p, por lo que se lee de manera ligeramente extraña. Segundo, la parte “BF=15.92” sólo tendrá sentido para las personas que ya entienden los métodos bayesianos, y no todos lo hacen. En tercer lugar, no está claro exactamente qué prueba se ejecutó y qué software se utilizó para hacerlo.

Por otro lado, a menos que la precisión sea sumamente importante, creo que esto está llevando las cosas un paso demasiado lejos:

Se realizó una prueba bayesiana de asociación usando la versión 0.9.10-1 del paquete BayesFactor usando priores por defecto y un plan de muestreo multinomial conjunto. El factor Bayes resultante de 15.92 a 1 a favor de la hipótesis alternativa indica que hay evidencia moderadamente fuerte para la no independencia de especies y elección.

Todo sobre ese pasaje es correcto, claro. Morey y Rouder (2015) construyeron sus pruebas bayesianas de asociación utilizando el papel de Gunel y Dickey (1974), la prueba específica que utilizamos asume que el experimento se basó en un plan de muestreo multinomial conjunto, y de hecho el factor Bayes de 15.92 es evidencia moderadamente fuerte. Es que es demasiado prolijo.

En la mayoría de las situaciones simplemente no necesitas tanta información. Mi preferencia suele ser ir por algo un poco más breve. Primero, si estás reportando múltiples análisis factoriales Bayes en tu escritura, entonces en algún lugar solo necesitas citar el software una vez, al comienzo de la sección de resultados. Entonces podrías tener una frase como esta:

Todos los análisis se realizaron usando el paquete BayesFactor en R y, a menos que se indique lo contrario, se utilizaron valores de parámetros predeterminados

¿Te das cuenta de que no me molesto en incluir el número de versión? Eso es porque la cita en sí incluye esa información (ve a revisar mi lista de referencias si no me crees). No hay necesidad de abarrotar sus resultados con información redundante que casi nadie necesitará realmente. Cuando llegues a la prueba real puedes salirte con la tuya con esto:

Una prueba de asociación produjo un factor Bayes de 16:1 a favor de una relación entre especies y elección.

Corto y dulce. He redondeado 15.92 a 16, porque no hay realmente ninguna diferencia importante entre 15. 92:1 y 16:1. Yo deletreé “factor Bayes” en lugar de truncarlo a “BF” porque no todos conocen la abreviatura. Indiqué exactamente cuál es el efecto (es decir, “una relación entre especies y elección”) y cuán fuerte fue la evidencia. No me molesté en indicar si se trataba de pruebas “moderadas” o “fuertes”, ¡porque las probabilidades mismas te lo dicen! No hay nada que te impida incluir esa información, y yo mismo lo he hecho en ocasiones, pero no la necesitas estrictamente. De igual manera, no me molesté en indicar que ejecuté el plan de muestreo “multinomial conjunto”, porque estoy asumiendo que la sección de método de mi escritura dejaría en claro cómo se diseñó el experimento. (Podría cambiar de opinión sobre eso si la sección del método fuera ambigua). Tampoco me molesté en indicar que se trataba de una prueba bayesiana de asociación: si tu lector no puede resolverlo por el hecho de que estás reportando un factor Bayes y el hecho de que estás citando el paquete BayesFactor para todos tus análisis, entonces no hay posibilidad de que ellos” Entenderé todo lo que hayas escrito. Además, si sigues escribiendo la palabra “Bayes” una y otra vez empieza a parecer estúpida. Bayes Bayes Bayes Bayes. ¿Ves?

Otros planes de muestreo

Hasta este punto todo lo que te he mostrado es cómo usar la función ContingencyTableBF () para el plan de muestreo multinomial conjunto (es decir, cuando el tamaño total de la muestra N es fijo, pero nada más lo es). Para el plan de muestreo de Poisson (es decir, nada fijo), el comando que necesita es idéntico excepto por el argumento SampleType:

contingencyTableBF(crosstab, sampleType = "poisson" )

## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 28.20757 ±0%
## 
## Against denominator:
##   Null, independence, a = 1 
## ---
## Bayes factor type: BFcontingencyTable, poisson

Observe que el factor Bayes de 28:1 aquí no es el idéntico al factor Bayes de 16:1 que obtuvimos de la última prueba. El plan de muestreo realmente importa.

¿Qué pasa con el diseño en el que se fijan las columnas de fila (o totales de columna)? Como mencioné anteriormente, esto corresponde al plan de muestreo “multinomial independiente”. Nuevamente, necesitas especificar el argumento SampleType, pero esta vez necesitas especificar si fijaste las filas o las columnas. Por ejemplo, supongamos que muestreé deliberadamente 87 humanos y 93 robots, entonces necesitaría indicar que el Margen Fijode la tabla de contingencia son las “filas”. Entonces el comando que usaría es:

contingencyTableBF(crosstab, sampleType = "indepMulti", fixedMargin="rows")

## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 8.605897 ±0%
## 
## Against denominator:
##   Null, independence, a = 1 
## ---
## Bayes factor type: BFcontingencyTable, independent multinomial

Nuevamente, el factor Bayes es diferente, con la evidencia de la alternativa cayendo a un mero 9:1. Como cabría esperar, las respuestas volverían a ser diferentes si fueran las columnas de la tabla de contingencia las que fijaba el diseño experimental.

Por último, si pasamos al muestreo hipergeométrico en el que todo está arreglado, obtenemos...

contingencyTableBF(crosstab, sampleType = "hypergeom")
#Error in contingencyHypergeometric(as.matrix(data2), a) : 
#  hypergeometric contingency tables restricted to 2 x 2 tables; see help for contingencyTableBF()

... un mensaje de error. Bien, algunas lecturas rápidas a través de los archivos de ayuda sugieren que se acerca el soporte para tablas de contingencia más grandes, pero aún no se ha implementado. Mientras tanto, imaginemos que tenemos datos del experimento de “etiquetado de juguetes” que describí anteriormente en esta sección. Específicamente, digamos que nuestros datos se ven así:

toys <- data.frame(stringsAsFactors=FALSE,
        gender = c("girl", "boy"),
        pink = c(8, 2),
        blue = c(2, 8)
        )

La prueba bayesiana con muestreo hipergeométrico nos da lo siguiente:

contingencyTableBF(toys, sampleType = "hypergeom")

#Bayes factor analysis
#--------------
#[1] Non-indep. (a=1) : 8.294321 @plusorminus0%
#
#Against denominator:
#  Null, independence, a = 1 
#---
#Bayes factor type: BFcontingencyTable, hypergeometric

El factor Bayes de 8:1 proporciona evidencia modesta de que las etiquetas estaban siendo asignadas de una manera que correlaciona el género con el color, pero no es concluyente.