1 INTRODUCCION

De acuerdo con lo expuesto en el apartado anterior, para conocer la relación entre dos variables cualitativas (nominales u ordinales) se pueden utilizar:

2 TABLAS DE CONTINGENCIA

Supongamos que queremos analizar las preferencias políticas en una región en la que existen 3 opciones políticas diferentes. Sospechamos que existe una diferencia acusada entre el ámbito urbano y el ámbito rural. El objetivo es determinar si esa hipótesis es verdad. Para ello, primero crearemos los datos de partida.

datos <- matrix(c(1200, 900, 400, 680, 250, 650), ncol=3, byrow=TRUE)
colnames(datos) <- c("P1","P2","P3")
rownames(datos) <- c("Rural","Urbano")

El análisis de dos variables cualitativas se realiza sobre una tabla de frecuencias, que es creada a continuación. Otra opción sería utilizar la función table() aplicada a un dataframe en el que contamos con dos variables cualitativas, como se realizó en el tema 2.

tabla1 <- as.table(datos)

Como no volveremos a utlizar los datos originales, los borramos del Global Environment.

rm(datos)

Figura 1: Ejemplo de tabla de contingencia

Las frecuencias relativas, en términos de proporciones (tantos por uno), son:

prop.table(tabla1)

Las cifras en la columna de la derecha y en la fila inferior (“Total”) reciben el nombre de Frecuencias marginales; la cifra situada en la esquina inferior derecha es el Gran total. Para añadir las frecuencias marginales a la tabla de contingencia, se añade una nueva fila con las sumas de cada columna y una nueva columna con las sumas de cada fila. Con R, esto se lleva a cabo con la función addmargins() cuya sintaxis básica es:

addmargins(tabla, margin=…, FUN=…)

donde:

El resultado es otro objeto de la clase table al que se le han añadido una o varias filas o columnas, que contienen las frecuencias marginales, tanto absolutas como relativas

addmargins(tabla1)                                          # Frecuencias marginales absolutas

addmargins(prop.table(tabla1))                              # Frecuencias marginales relativas

En caso de que necesitásemos únicamente la fila o la columna de frecuencias marginales, las funciones colSums() y rowSums() crean tablas sólo con la suma por columnas o por filas, respectivamente.

colSums(tabla1)
rowSums(tabla1)

3 GRÁFICOS

Para la visualización de las tablas de contingencia es posible utilizar diferentes tipos de gráficos.

3.1 Gráfico sencillo

La función balloonplot()(paquete gplots) dibuja una matriz gráfica donde cada celda contiene un punto cuyo tamaño refleja la magnitud relativa del componente correspondiente.

library("gplots")
balloonplot(tabla1, 
            main ="Relación entre género y uso de móvil",                       # Título del gráfico
            xlab ="género",                                                     # Etiqueta del eje X
            ylab="procedencia",                                                 # Etiqueta del eje Y
            label = TRUE,                                                       # Dibuja el valor numérico de los componentes de la tabla.
            show.margins = TRUE)                                                # Dibuja las frecuencias marginales.

3.2 Gráfico de barras apiladas/superpuestas

Como ya sabemos, se obtiene cuando se aplica la función barplot() a una tabla con dos variables (bidimensional), con dos posibilidades dependiendo del argumento beside:

  • si es “FALSE” barras apiladas.

  • si “TRUE” barras yuxtapuestas.

OPCION 1: barras apiladas (por defecto)

par(mar = c(5, 5, 4, 11))
barplot(tabla1,
        main = "Gráfico de barras apiladas",
        sub = "Subtítulo",
        xlab = "Etiqueta eje X",
        ylab = "Etiqueta eje Y",
        col = c("red", "blue"),
        axes = TRUE, 
        ylim = c(0, 2000),
        legend.text = rownames(tabla1),                                         # Texto de la leyenda
        args.legend = list(x = "topright", inset = c(-0.30, 0)))                # Posición de la leyenda

par(mar = c(5, 4, 4, 2))

OPCION 2: barras yuxtapuestas

barplot(tabla1,
      main = "Gráfico de barras yuxtapuestas",
      sub = "Subtítulo",
      ylim = c(0,1500),
      xlab = "Etiqueta eje X",
      ylab = "Etiqueta eje Y",
      axes = TRUE, 
      beside = TRUE)

Desafío

A partir del dataframe zonas_verdes muestra la relación entre las variables obras y barrio a partir de

  • Un gráfico creado con la función ballplot.

  • Un gráfico de barras yuxtapuestas que incorpore un título, una leyenda y las etiquetas de los ejes X e Y.

balloonplot(tabla2, 
            main ="Relación entre género y uso de móvil",
            xlab ="obras",
            ylab="barrio",
            label = FALSE, 
            show.margins = FALSE)

barplot(tabla2,
      main = "Relación entre género y uso de móvil",
      sub = "Subtítulo",
      ylim = c(0,20),
      xlab = "Barrios",
      ylab = "Número de parques",
      axes = TRUE, 
      beside = TRUE,
      legend.text = rownames(tabla2),
      args.legend = list(x = "topright"))

3.3 Diagrama de mosaico

Es un diagrama que facilita la representación visual de la asociación entre dos o más variables categóricas. En este caso, cada entrada de la tabla de frecuencias corresponde a una región rectangular de área proporcional a su valor, utilizando para ello la función mosaicplot() del paquete graphics.

Como argumentos adicionales se incluye:

  • shade usado para colorear el gráfico.

    • Color azul: indica que los valores observados son más elevados que los esperados si estos fueran aleatorios

    • Color rojo: indica que los valores observados son más bajos que los esperados.

  • las = 2 dispone verticalmente las etiquetas.

library("graphics")
mosaicplot(tabla1,
           shade = TRUE, las = 1)

Existe otro paquete, denominado vcd, que también puede ser usado para elaborar un gráfico de mosaico con la función mosaic()

library("vcd")

Para la elaboración del diagrama de mosaico hay que señalar que la primera variable mencionada en la sintaxis corresponde a las columnas (vertical), mientras que la segunda variable corresponde a las filas (horizontal).

mosaic(genero ~ uso_movil_let, 
       data = df, rot_labels=c(0,90,0,0))

La información proporcionada por este diagrama puede mejorarse.

mosaic(~genero + uso_movil_let, 
       data = df,
       main = "Título", 
       shade = TRUE,                                  # Colorea la figura
       legend = TRUE)                                 # Muestra una leyenda con los resíduos de Pearson.

La leyenda representada en el gráfico de mosaico muestra cuáles son las celdas más relevantes para los resultados de la prueba de chi-cuadrado. Contiene las mismas gamas de colores que en el caso del paquete graphics.

  • Color azul: significa que hay más observaciones en esa celda de lo que se esperaría bajo el modelo nulo (independencia).

  • Color rojo (no aparece en este ejemplo) significa que hay menos observaciones de las esperadas.

3.4 Spine plot

Los spineplots son una generalización de los gráficos de barras apilados. Un gráfico de mosaico es un spineplot para dos o más variables.

ATENCIÓN: se requiere una tabla de frecuencias

spineplot(tabla1)

Este gráfico puede ser personalizado, por ejemplo, trasponiendo las variables

spineplot(t(tabla1))                                            # Trasposición de variables

También puede modificarse el color de los cuadros y el borde de éstos

spineplot(tabla1, 
          col = c("red", "yellow", "blue"),
          border = "darkblue")                                      # Color del borde de las barras

Por defecto, la función spineplot crea un spineplot, pero el argumento off = 0 crea un espinograma, eliminando la distancia entre las barras.

spineplot(tabla1, off = 0)

4 PRUEBAS DE ASOCIACIÓN

Una vez tabulados los datos y visualizadas las posibles diferencias entre variables, el tercer paso es verificar estadísticamente esas posibles diferentes. Para ello recurrimos a varias pruebas.

4.1 La prueba de contraste de independencia con \(\chi^2\).

Una vez elaborada la tabla, es posible evaluar la significación estadística de la diferencia entre ellas con la prueba \(\chi^2\) de Pearson. Esta prueba tiene varios usos:

  • Opción 1: Evaluar si hay una asociación significativa entre dos o más variables cualitativas en una tabla de frecuencias.

  • Opción 2: determinar si un variable categórica sigue una distribución hipotética de valores esperados obtenidos de acuerdo a una hipótesis de proporcionalidad (bondad de ajuste).

En el caso que nos ocupa se utiliza la primera opción, siendo denominada prueba de “contraste de independencia” para tablas de contingencia bidimensionales. Consiste en decidir si las dos variables de la tabla tienen distribuciones independientes, es decir, si la distribución de probabilidades conjunta es igual al producto de las probabilidades marginales. Si la proporción de individuos en cada columna varía entre las diversas filas y viceversa, se dice que existe asociación entre las dos variables. Si no existe asociación se dice que ambas variables son independientes.

En nuestro ejemplo, queremos decidir si podemos aceptar que las variables garaje y terraza son independientes o si por el contrario hay evidencia de que la posibilidad de disponer de un garaje depende de la presencia o no de terraza (y esto podría atribuirse a la calidad de la vivienda).

La prueba de \(\chi^2\) presenta algunas características singulares:

  • Test no paramétrico: no requiere el conocimiento de ciertas características de la población de la que se han extraído las muestras.

  • Las categorías deben ser por lo menos dos y excluyentes (ningún individuo puede estar en más de una categoría a la vez.

Además, requiere varios pasos.

  • PRIMER PASO: definir las hipótesis nula y alternativa

    • Hipótesis nula (H0): las dos variables son independientes (no hay relación entre ellas).

    • Hipótesis alternativa (H1): las dos variables están relacionadas (no son independientes).

Además, debe decidirse el valor alfa. Esto implica decidir qué riesgo se asume para llegar a una conclusión errónea. Normalmente α=0,05 para las pruebas de independencia, es decir, existe un riesgo del 5 % de concluir que las dos variables son independientes, cuando en realidad no lo son. En este caso:

  • Si el valor-p α>0,05 se acepta la Hipótesis nula (H0): las dos variables son independientes (no hay relación entre ellas).

  • Si el valor-p α<0,05 se rechaza la hipótesis nula y se acepta Hipótesis alternativa (H1): las dos variables están relacionadas (no son independientes).

  • SEGUNDO PASO: cálculo de los valores esperados de cada celda de la tabla bajo la hipótesis nula, de acuerdo con la fórmula

Figura 2: Fórmula valores_esperados

Para ello, en R podemos realizar las siguientes operaciones. En primer lugar, extraemos las frecuencias marginales:

n <- length(tabla1)
            
frec.mar.fil <- rowSums(tabla1)                           # Frecuencias marginales de las filas
frec.mar.col <- colSums(tabla1)                           # Frecuencias marginales de las columnas

El test de independencia usa las frecuencias absolutas esperadas bajo la hipótesis nula de independencia, que se obtienen, para cada celda (i,j), multiplicando la frecuencia marginal de la fila i por la de la columna j y dividiendo por el tamaño de la muestra. En nuestro ejemplo estas frecuencias esperadas son

Figura 3: Cálculo de los valores esperados

Por lo tanto, con R obtenemos esta tabla de frecuencias esperadas de la manera siguiente:

frec.esperadas <- frec.mar.fil%*%t(frec.mar.col)/n
frec.esperadas
  • TERCER PASO: cálculo del estadístico comparando los valores observados (O) y esperados (E)

Figura 4: Fórmula chi-cuadrado

Una precaución antes de efectuar el test de independencia es necesario comprobar que las frecuencias esperadas son mayores o iguales que 5. Si se cumple esta condición se puede calcular de forma directa o aplicando la función chisq.test a la tabla de contingencia con las frecuencias absolutas.

La sintaxis de forma directa es:

chi2.estadistico <- sum((tabla1-frec.esperadas)^2/frec.esperadas)
chi2.estadistico

El p-valor del contraste se calcula así

p.valor <- 1-pchisq(chi2.estadistico,df=(dim(tabla1)[1]-1)*(dim(tabla1)[2]-1))
p.valor

Mediante la función de R

chisq.test(tabla1)

Como el p-valor es inferior a 0.05, se rechaza la hipótesis nula y se acepta la alternativa: las variables objeto de estudio están relacionadas.

La función chisq.test() ofrece una serie de outputs adicionales, como son los siguientes

chisq.test(tabla1)$observed         # Valores observados
chisq.test(tabla1)$expected         # Valores esperados bajo la hipótesis nula
chisq.test(tabla1)$parameter        # Grados de libertad
chisq.test(tabla1)$p.value          # Probabilidad

qchisq(0.95,                        # Valor critico de chi-cuadrado, al 95% con 1 grado de libertad
       1)                           # Grados de libertad, calculados a partir de n-1 (numero de categorias menos uno).

En el caso de tablas formadas por 2 filas y 2 columnas, se debe aplicar una corrección de continuidad.

chisq.test(tabla1, 
           correct=FALSE)                # Indicador lógico para aplicar una corrección de continuidad

Si algunas frecuencias absolutas esperadas fueran inferiores a 5, la aproximación del p-valor por una distribución \(\chi^2\) podría no ser adecuada. En este caso, al ser las variables cualitativas, no podemos recurrir al agrupamiento de valores consecutivos, puesto que no tienen orden. Si se da esta situación, lo mejor es recurrir a simular el p-valor usando el argumento simulate.p.value=TRUE.

set.seed(NULL)
chisq.test(tabla1,
           correct=FALSE,
           simulate.p.value=TRUE,           # Calcula los valores de probabilidad mediante una simulación Monte Carlo
           B=5000)$p.value                  # Número de réplicas simuladas y probabilidad

Otra opción en caso de frecuencias menores a 5 en tablas de contingencia formadas por 2 filas y 2 columnas es el test de Fisher.

?fisher.test

La prueba exacta de Fisher utiliza las siguientes hipótesis nula y alternativa:

  • Hipótesis nula H0: las dos variables son independientes.

  • Hipótesis alternativa H1: las dos variables no son independientes (existe relación entre ellas).

fisher.test(tabla1)

En la prueba exacta de Fisher, la hipótesis nula es que las dos columnas son independientes (o, de manera equivalente, que la odds ratio -razón de posibilidades- es igual a 1).

Para determinar si las dos columnas son independientes, se utiliza el p-value de la prueba. Si el valor es inferior a 0.05, no rechazamos la hipótesis nula: no tenemos evidencia suficiente para decir que existe una asociación estadísticamente significativa entre ambas variables.

fisher.test(tabla1,
            simulate.p.value=TRUE, 
            B=2^3)

Si se desea conocer qué celdas contribuyen más a los valores de la prueba, se calculan los resíduos de Pearson correspondientes a cada celda (residuos estandarizados)

round(chisq.test(tabla1)$residuals, 3)

Finalmente, es posible visualizar estos resíduos: el tamaño del círculo es proporcional a la contribución de cada celda.

library(corrplot)
corrplot(chisq.test(tabla1)$residuals,is.cor = FALSE)

4.2 Otros coeficientes

El grado de asociación entre dos variables se puede cuantificar adicionalmente empleando otros coeficientes:

  • El coeficiente phi es el valor de la chi cuadrado entre el número de observaciones. En general, cuanto más lejos de 0, más fuerte será la relación entre las dos variables.

    • Un valor próximo a 0 indica independencia entre los factores.

    • Valores próximos a -1 indica una relación perfectamente negativa entre las dos variables.

    • Valores próximos a 1 indica una relación perfectamente positiva entre las dos variables.

  • El coeficiente de contingencia también es una medida de la intensidad de la relación basado en la chi cuadrado y toma valores entre 0 (independencia) y 1 (dependencia).

  • La V de Cramer es menos susceptible a valores muestrales.

    • 0 implica independencia

    • 1 una relación perfecta entre los factores. Habitualmente valores superiores a 0,30 ya nos están indicando que hay una posible relación entre las variables; un valor de 0.5 o superior ya indicaría una gran relación entre ambas variables.. En este caso tenemos valores muy próximos a 0 en todos los estadísticos que nos ofrece, ambos factores son independientes. Saludos.

library(vcd)
assocstats(tabla1)

Desafío

PRIMER EJERCICIO

  • Importa el fichero alumnos_UC.RData. Señala qué variables son cuantativas y qué variables son cualitativas. A partir de ese dataframe:
load("D:/G14_EE_CC_SS_2023/TEMA_2_Estadistica_Descriptiva/alumnos_UC.Rdata")

df$genero <- as.factor(df$genero)
  • Crea una tabla que relacione el género y el uso del móvil.
tabla2 <- table(df$genero, df$uso_movil_let)
  • Calcula las frecuencias relativas en tantos por ciento.
prop.table(tabla2) * 100
  • Calcula por separado las frecuencias marginales por filas y por columnas.
addmargins(tabla_nueva, margin = 1) 
addmargins(tabla_nueva, margin = 2)
  • Elabora un gráfico que represente las dos variables.
par(mar = c(5, 5, 4, 11))
barplot(tabla_nueva,
        main = "Gráfico de barras apiladas",
        sub = "Subtítulo",
        xlab = "Etiqueta eje X",
        ylab = "Etiqueta eje Y",
        col = c("red", "blue"),
        axes = TRUE, 
        ylim = c(0, 60),
        legend.text = rownames(tabla_nueva),                                         # Texto de la leyenda
        args.legend = list(x = "topright", inset = c(-0.30, 0)))                # Posición de la leyenda
par(mar = c(5, 4, 4, 2))
  • Analiza si existe relación entre el género y el uso del móvil aplicando la prueba de Chi cuadrado.
chisq.test(tabla_nueva)

SEGUNDO EJERCICIO

  • Importa el fichero zonas_verdes.Rdata. A partir de ese dataframe:
load("D:/G14_EE_CC_SS_2023/TEMA_2_Estadistica_Descriptiva/zonas_verdes.Rdata")
str(zonas_verdes)
  • Crea una nueva variable recodificando la variable fecha_inauguración en dos categorías: los parques inaugurados después del 1 de enero de 1990 serán clasificados como “Nuevos” mientras que los parques inaugurados con posterioridad serán clasificados como “Antiguos”. Esta nueva variable, que pueden denominar periodo deberá ser convertida en factor para su uso posterior.
zonas_verdes <- within(zonas_verdes,{
  periodo <- NA                                                 # Crea la variable grupo, pero vacía  
  periodo[fecha_inauguracion > "1990-01-01"] <- "Nuevos"        # Los restantes comandos rellenan esa variable de manera ordenada
  periodo[fecha_inauguracion <= "1990-01-01"] <- "Antiguos"})

zonas_verdes$periodo <- as.ordered(zonas_verdes$periodo)
zonas_verdes$obras <- factor(zonas_verdes$obras,
                             levels = c(0,1),
                             labels = c("Sin obras", "Con obras"))
  • Crea una tabla que relacione el periodo y las obras, incorporando las frecuencias marginales.
tabla_nueva <- addmargins(table(zonas_verdes$obras, zonas_verdes$periodo), margin = c(1,2)) 
tabla_nueva
  • Elabora el gráfico que consideres más conveniente.
mosaicplot(tabla_nueva,
           shade = TRUE, 
           las = 1)
  • Calcula si existe relación entre la antigüedad del parque y el que ahora mismo se encuentre en obras mediante la prueba del chi-cuadrado y con la V de Cramer.
chisq.test(tabla_nueva)
assocstats(tabla_nueva)

ACTIVIDADES DE EVALUACIÓN CONTINUA

El fichero word 03_Actividad_evaluacion_continua_1_2023 contiene una serie de actividades encaminadas a poner en práctica los conocimientos adquiridos. Cada alumnos deberá elaborar un script que contenga los procedimientos necesarios para responder a las cuestiones planteadas, y enviarlo posteriormente al profesor a través de correo electrónico.

Los resultados de esta actividad pueden consultarse en 03_Actividad_evaluacion_continua_1_2023_resuelta