Tablas de contingencia y contraste de independencia (test chi cuadrado)

Contraste de independencia

El Barómetro del CIS  (Centro de Investigaciones Sociológicas) permite obtener datos sobre las creencias religiosas de la población en España. Una pregunta que puede interesarnos es ¿hay alguna diferencia al respecto entre hombres y mujeres? Vamos a utilizar los datos del Barómetro para intentar contestar esa pregunta.

Por ejemplo, en el mes de enero de 2013 el Barómetro} recoge las respuestas de n=2452 personas sobre sus creencias religiosas. (En realidad son 2483, pero para simplificar vamos a eliminar de nuestra consideración a las 19 mujeres y a los 12 hombres que decidieron no contestar). Observa que, como de costumbre, vamos a usar n para el número total de personas encuestadas. Agrupamos a todos los creyentes de distintas religiones por un lado y a los que se declaran no creyentes o ateos por otro. Y así tenemos una tabla de doble entrada:

Hombres Mujeres Total
Creyentes ?? ?? 1864
No creyentes ?? ?? 588
Total 1205 1247 2452

En la que hemos dejado pendientes los valores interiores de la tabla, y sólo se muestran los valores marginales. Si las creencias religiosas no dependieran del género, ¿cuáles son los valores que esperaríamos ver en esta tabla? Los que se deducen de los marginales, claro. Por ejemplo, el número de hombres creyentes se puede calcular multiplicando el número de hombres por la proporción de creyentes en la población. Esto es:

1205 \cdot \frac{1864}{2452}\approx 916

Con R podemos obtener estos valores esperados fácilmente. Como es la primera vez que hacemos algo así en el blog, vamos a dar los detalles. Empezamos con dos vectores, los de los valores marginales:

creencia=c(1864,588)
genero=c(1205,1247)

Y ahora los multiplicamos como matrices, trasponiendo el vector creencia, que va en vertical en la tabla, y que por tanto hace el papel de vector columna. Para trasponer un vector basta con usar la función t(). El producto matricial se representa mediante el símbolo %*%, así que la tabla (matriz producto) es (mostramos la salida):

> creencia %*% t(genero)
        [,1]    [,2]
[1,] 2246120 2324408
[2,]  708540  733236

Naturalmente, aún tenemos que dividir por el total de personas (el número n). El resultado completo es esta tabla de valores esperados:

> (tablaEsperada=creencia %*% t(genero) / sum(genero))
         [,1]     [,2]
[1,] 916.0359 947.9641
[2,] 288.9641 299.0359

Está bien, pero sería más fácil de interpretar si añadimos nombres a las filas y columnas, y además incluimos los valores marginales con los que habíamos empezado. Hacemos esto con los comandos:

colnames(tablaEsperada)=c("H","M")
rownames(tablaEsperada)=c("CREY","NO_CREY")
(addmargins(tablaEsperada))

y el resultado es mucho más agradable de ver:

                H         M  Sum
CREY     916.0359  947.9641 1864
NO_CREY  288.9641  299.0359  588
Sum     1205.0000 1247.0000 2452

De acuerdo, eso es lo que esperaríamos, si las creencias y el género fueran independientes. ¿Cuáles son los datos reales que se obtuvieron en el Barómetro? Aquí está la tabla de valores observados:

Hombres Mujeres Total
Creyentes 849 1015 1864
No creyentes 356 232 588
Total 1205 1247 2452

Introducimos esta tabla en R directamente usando cbind, y le añadimos la información necesaria para hacerla más útil:

(tablaObservada=rbind(c(849,1015),c(356,232)))
colnames(tablaObservada)=c("H","M")
rownames(tablaObservada)=c("CREY","NO_CREY")
(addmargins(tablaObservada))

El resultado es:

           H    M  Sum
CREY     849 1015 1864
NO_CREY  356  232  588
Sum     1205 1247 2452

Comparando las dos tablas, resulta evidente que los valores observados no coinciden con los  esperados. De hecho, el número de hombres no creyentes es más alto de lo que habíamos estimado a partir de la población en conjunto (y,  lógicamente, el número de mujeres no creyentes es más bajo que la estimación). Pero ese número de hombres no creyentes, ¿es significativamente más alto?

Para contrastar la hipótesis de independencia entre género y creencias, usaremos la función chisq.test de R (se muestra la salida):

> (chisq.Observada=chisq.test(tablaObservada,correct=FALSE))

	Pearson's Chi-squared test

data:  tablaObservada
X-squared = 40.2253, df = 1, p-value = 2.263e-10

Como en otros casos, hemos incluido la opción correct=FALSE para que R (lo haga un poco peor de lo que sabe y) nos devuelva el mismo valor que aparecería en un curso de Introducción a la Estadística.  Se recomienda leer la ayuda de chisq.test para ver lo que implica esa opción.

En este caso el p-valor es muy pequeño, así que con esos datos deberíamos rechazar la hipótesis nula, y concluir que el género y las creencias no son independientes.

El test \chi^2 que estamos usando consiste, para este ejemplo, en calcular el siguiente estadístico:

X^2=\dfrac{(o_{11}-e_{11})^2}{e_{11}}+\dfrac{(o_{12}-e_{12})^2}{e_{12}}+\dfrac{(o_{21}-e_{21})^2}{e_{21}}+\dfrac{(o_{22}-e_{22})^2}{e_{22}}.

donde los valores o_{ij} son los de la tabla observada, y los e_{ij} son los valores esperados. Cada término de la forma

\frac{o_{ij}-e_{ij}}{e_{ij}}

es un residuo, y el estadístico es la suma de los residuos al cuadrado. En R,  al aplicar el test, el programa devuelve,  dentro de la variable chisq.Observada, una lista con todo lo  necesario. Accedemos a las componentes de esa lista con el operador $. Por
ejemplo, las tablas de valores esperados y observados se obtienen así (se muestra la salida):

> (E=chisq.Observada$expected)
               H        M
CREY    916.0359 947.9641
NO_CREY 288.9641 299.0359

> (O=chisq.Observada$observed)
          H    M
CREY    849 1015
NO_CREY 356  232

Y a partir de aquí el cálculo de los residuos es inmediato:

> (residuos=(O-E)/sqrt(E))
                H         M
CREY    -2.214885  2.177266
NO_CREY  3.943532 -3.876553

Lo hemos hecho de esta manera para que quede claro cómo se definen, pero los residuos también se pueden obtener directamente con:

chisq.Observada$res

De cualquier manera, para calcular a mano el estadístico a partir de los residuos, podemos ahora elevar los residuos al cuadrado y sumarlos usando addmargins, de esta manera:

> (tablaChi2=residuos^2)
                H         M
CREY     4.905714  4.740486
NO_CREY 15.551448 15.027663

> addmargins(tablaChi2)
                H         M      Sum
CREY     4.905714  4.740486  9.64620
NO_CREY 15.551448 15.027663 30.57911
Sum     20.457163 19.768148 40.22531

El estadístico es el valor de la esquina inferior derecha de esta tabla:

> (estadistico=addmargins(tablaChi2)[3,3])
[1] 40.22531

Y podemos ver que coincide con lo que obtuvimos antes usando chisq.test. El p-valor se obtiene a partir del estadístico usando pchisq como en cualquier contraste (con un grado de libertad, porque es una tabla 2×2). Usamos la cola derecha, porque los valores del estadístico que contradicen a la hipótesis nula se sitúan ahí:

> (pValor=1-pchisq(estadistico,df=1))
[1] 2.262972e-10

Y como antes, el p-valor coincide con lo que obtuvimos en chisq.test.

Una posible representación gráfica de este contraste es mediante un gráfico de columnas apiladas (la altura indica el porcentaje), con una columna por  género. Para obtener esa representación empezamos por fabricar una tabla de proporciones (porcentajes), en la que los porcentajes son sobre el total de la muestra. Añadimos las sumas marginales para que sea más fácil entender esta tabla (recordamos que se muestra la salida):

> tablaProporciones1=prop.table(tablaObservada)
> addmargins(tablaProporciones1)
                H          M       Sum
CREY    0.3462480 0.41394780 0.7601958
NO_CREY 0.1451876 0.09461664 0.2398042
Sum     0.4914356 0.50856444 1.0000000

Para que cada columna de la tabla nos de las alturas de las columnas del gráfico, tenemos que dividir la tabla por las sumas por columnas, calculadas con colSums. Además, vamos a redondearla para mejorar la presentación.

> (tablaProporciones2=round(tablaObservada/colSums(tablaObservada),digits=1))
          H   M
CREY    0.7 0.8
NO_CREY 0.3 0.2

Ya podemos usar barplot para dibujar el grafico, con colores y etiquetas en los ejes. Aunque con esto ya se obtiene un gráfico adecuado, la interpretación de un gráfico como este es mucho más fácil, si se rotulan las columnas con el porcentaje de una de las categoría en cada columna. Primero fabricamos los rótulos, y luego los colocamos en sus posiciones con text. Estos son los comandos:

(bp_t2=barplot(tablaProporciones2,col=c("lightseagreen","darkgoldenrod1"),xlab="Género",ylab="Creyente/No creyente",main="Creencias religiosas por género",font=2))

rotulos=paste(format(100*signif(tablaProporciones2[1,],digits=2)),"% creyentes",sep="")

text(bp_t2, tablaProporciones2[1,]+0.05,rotulos , xpd = TRUE, font=2)

Y este es el resultado:

20130514-CreenciasReligiosasPorGenero

Un problema como este, con una tabla de contingencia 2×2, también se puede abordar con un contraste de igualdad de proporciones para dos muestras, como los que vimos en esta entrada. Vamos a hacerlo aquí (de nuevo, con correct=FALSE), para que se vea que R devuelve el mismo valor del estadístico, y el mismo p-valor.

> (creyentesObs=tablaObservada[1,])
   H    M
 849 1015
> prop.test(creyentesObs,genero,correct=F)

	2-sample test for equality of proportions without continuity correction

data:  creyentesObs out of genero
X-squared = 40.2253, df = 1, p-value = 2.263e-10
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.14300581 -0.07577254
sample estimates:
   prop 1    prop 2
0.7045643 0.8139535

Por último, si tenemos dudas sobre las hipótesis de normalidad subyacentes al test \chi^2, podemos hacer un test no paramétrico para confirmar lo anterior.

> fisher.test(tablaObservada)

	Fisher's Exact Test for Count Data

data:  tablaObservada
p-value = 2.745e-10
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.4490685 0.6611693
sample estimates:
odds ratio
 0.5452314

que, en este ejemplo, confirma que el resultado del test es significativo.

Código en R y observaciones finales

En este enlace hay un fichero con el código R que hemos usado en este ejemplo. Naturalmente, el caso 2×2 no agota las posibilidades del test \chi^2 de independencia. Estoy terminando de escribir el correspondiente capítulo de los apuntes de teoría (avisaré en un post cuando acabe), que incluye ejemplos y código para el caso general, con tablas de contingencia con un número mayor de filas y columnas.

Aparte de esto, y como se explica en ese capítulo, existe otra variante del test, que se aplica a problemas de bondad de ajuste, es decir a la pregunta de si los datos observados se corresponden con una distribución teórica. Concretamente, ejemplos de este tipo de preguntas son: ¿cómo podemos saber si un dado está cargado? O ¿corresponden estos datos a una población normal?.  Próximamente le dedicaremos una entrada del blog a este tema. Y, puesto que aquí ha surgido, siquiera sea de pasada, el tema de las operaciones con matrices, también le dedicaremos una entrada a ese tema.

Gracias por la atención.

Y una última novedad (si P. Krugman se atreve, no voy a  ser yo quien se quede atrás): Recomendación musical de esta entrada Headlock de Imogen Heap,

Anuncios

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s