Análisis de correspondencias con R: aplicación a datos de encuestas

II Seminario Análisis de datos avanzados en Ciencias de la Salud, Facultad de Ciencias de la Salud de la Universidad Rey Juan Carlos, Campus de Alcorcón.

Sesión impartida por Emilio López Cano, 24 de mayo de 2018

Las diapositivas del seminario se encuentran en: http://emilio.lcano.com/p/seminariourjc18/

Más ejemplos y explicaciones en mi libro de apuntes Análisis de datos con R (licencia CC)

Preparando el entorno

Descarga e instala R y RStudio en tu sistema, y en ese orden. Puedes encontrar las instrucciones y los archivos de instalación en las siguientes direcciones:

Después de instalar R y RStudio, abre RStudio y ejecuta el siguiente código en la consola de RStudio:

install.packages(c("usethis", "readxl", "FactoMineR", "factoextra", 
                   "dplyr", "gplots", "corrplot", "knitr"))
usethis::use_course("https://goo.gl/p3DU1Y")

En la consola se mostrarán mensajes de confirmación para descargar los materiales en un fichero zip. Tras confirmar, el fichero se descarga y descomprime automáticamente en la carpeta de usuario y se abre el proyecto RStudio con el que trabajaremos, incluido código y datos, que hay que descomprimir:

unzip("datos.zip")

Descripción e importación de datos

Vamos a trabajar con la última Encuesta europea de salud en España (EESE) que elabora el Instituto Nacional de Estadística. Se realiza con periodicidad quinquenal, y la más reciente disponible es la de 2014. Esta encuesta contiene bastantes variables cualitativas adecuadas para realizar análisis de correspondencias, y obviamente de interés para la investigación en Ciencias de la Salud. En la web del INE podemos encontrar y descargar los documentos metodológicos:

El INE publica los resultados de la encuesta después de su tratamiento. Pero además, el pone a disposición los microdatos de la encuesta, lo que es realmente interesante para la investigación académica. Los microdatos contienen las respuestas de cada cuestionario, con un formato determinado que se detalla en los documentos que acompañan a los datos. En el caso de la EESE, disponemos de ficheros en dos formatos:

  • Texto plano TXT con las columnas de ancho fijo

  • Formato Excel, una columna para cada variable

La descripción de las variables y posición en el fichero txt se encuentran en el fichero Excel de diseño de registro. Para este seminario vamos a importar los datos del fichero excel, cuya importación en R es trivial con el paquete readxl:

library(readxl)
eese <- read_excel("datos/MICRODAT.ADULTO.xlsx")

Tras la importación, tenemos un data frame en el espacio de trabajo de R con el que podemos trabajar. Contiene 22842 observaciones (encuestas) de 429 variables (respuestas a preguntas del cuestionario):

dim(eese)
## [1] 22842   429

Para este seminario vamos a seleccionar algunas variables cualitativas como ejemplo. Siguiendo las mismas pautas que se dan a continuación se pueden analizar cualesquiera otras variables de interés. Del mismo modo, el tratamiento y análisis de otras encuestas, ya sean microdatos publicados o de encuestas propias realizadas en las consultas médicas u hospitales, se realizaría de forma similar. Lo más eficiente es guardar los datos originales de la encuesta en un fichero excel con una fila por encuesta y una columna por variable e importar los datos a un data frame de R.

Seleccionamos las variables:

  • E4: Convivencia en pareja
  • G21: Estado de salud percibido en los últimos 12 meses
  • T112: Frecuencia con la que realiza alguna actividad física en su tiempo libre
  • V121: ¿Fuma actualmente?
  • W127: Frecuencia de consumo de alcohol en los últimos 12 meses

La codificación de las preguntas se deduce del propio cuestionario: letra del apartado y nº de pregunta, que se puede comprobar también en el fichero de diseño de registro. En este punto es importante señalar que los microdatos de encuestas suelen incluir una variable de homogeneización, llamada “factor de elevación”, que viene a representar la ponderación que se le asigna a cada individuo que contesta la encuesta. Para facilitar la exposición no tendremos en cuenta este factor de elevación en los ejemplos, que sí habría que tener en cuenta, junto con el resto de la metodología aplicada, en una investigación rigurosa.

Preparación de los datos

Vamos a trabajar con un data frame que contiene solo las variables de interés:

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
datos <- eese %>% 
  select(convivencia = E4, 
         estado.salud = G21, 
         actividad = T112,
         fuma = V121,
         alcohol = W127)

Que tiene la siguiente estructura:

str(datos)
## Classes 'tbl_df', 'tbl' and 'data.frame':    22842 obs. of  5 variables:
##  $ convivencia : chr  "3" "1" "1" "1" ...
##  $ estado.salud: chr  "4" "1" "2" "2" ...
##  $ actividad   : chr  "1" "2" "2" "1" ...
##  $ fuma        : chr  "4" "4" "3" "1" ...
##  $ alcohol     : chr  "01" "04" "04" "01" ...

Es decir, 22842 observaciones de 5 variables, todas ellas tipo carácter (chr). Esta sería una muestra de las primeras filas:

knitr::kable(head(datos, 10))
convivencia estado.salud actividad fuma alcohol
3 4 1 4 01
1 1 2 4 04
1 2 2 3 04
1 2 1 1 01
1 1 3 3 04
1 2 2 3 01
3 1 3 4 09
3 2 1 1 06
3 3 4 4 04
1 2 1 1 03

Vemos que los valores de las variables son códigos, cuyo significado aparece tanto en la encuesta como en el archivo de diseño de registro. Las variables cualitativas en R son de tipo factor, en vez de carácter, por lo que vamos primero a convertir estas variables en factores:

datos <- datos %>% transmute_all(as.factor)

Ahora cada variable es un factor con varios niveles:

str(datos)
## Classes 'tbl_df', 'tbl' and 'data.frame':    22842 obs. of  5 variables:
##  $ convivencia : Factor w/ 5 levels "1","2","3","8",..: 3 1 1 1 1 1 3 3 3 1 ...
##  $ estado.salud: Factor w/ 5 levels "1","2","3","4",..: 4 1 2 2 1 2 1 2 3 2 ...
##  $ actividad   : Factor w/ 6 levels "1","2","3","4",..: 1 2 2 1 3 2 3 1 4 1 ...
##  $ fuma        : Factor w/ 6 levels "1","2","3","4",..: 4 4 3 1 3 3 4 1 4 1 ...
##  $ alcohol     : Factor w/ 11 levels "01","02","03",..: 1 4 4 1 4 1 9 6 4 3 ...

Pero los códigos de los niveles no nos dicen mucho. Vamos a etiquetar estos niveles de forma más descriptiva para que las salidas y gráficos sean interpretables fácilmente. Vamos a asignar como valores perdidos (NA) los valores “No sabe” o “No contesta” en las encuestas. De nuevo, esto lo hacemos para centrarnos en la explicación del análisis de correspondencias. En una investigación rigurosa habría que analizar el efecto de estos valores faltantes.

levels(datos$convivencia) <- c("cónyuge", "pareja", "no", NA, NA)
levels(datos$estado.salud) <- c("Muy bueno", "Bueno", "Regular", "Malo", "Muy malo")
levels(datos$actividad) <- c("ninguna", "ocasional", "mensual", "semanal", NA, NA)
levels(datos$fuma) <- c("mucho", "poco", "ex", "nunca", NA, NA)
levels(datos$alcohol) <- c("diario", "semanal", "semanal", "semanal", "mensual", 
                           "mensual", "anual", "ex", "nunca",  NA, NA)

Quedando nuestro data.frame final en:

str(datos)
## Classes 'tbl_df', 'tbl' and 'data.frame':    22842 obs. of  5 variables:
##  $ convivencia : Factor w/ 3 levels "cónyuge","pareja",..: 3 1 1 1 1 1 3 3 3 1 ...
##  $ estado.salud: Factor w/ 5 levels "Muy bueno","Bueno",..: 4 1 2 2 1 2 1 2 3 2 ...
##  $ actividad   : Factor w/ 4 levels "ninguna","ocasional",..: 1 2 2 1 3 2 3 1 4 1 ...
##  $ fuma        : Factor w/ 4 levels "mucho","poco",..: 4 4 3 1 3 3 4 1 4 1 ...
##  $ alcohol     : Factor w/ 6 levels "diario","semanal",..: 1 2 2 1 2 1 6 3 2 2 ...

Asociación de variables

Vamos a estudiar el factor estado de salud con el resto de atributos, para ver en cuál o cuáles hay alguna relación. Tomemos a modo de ejemplo el consumo de alcohol, veamos en primer lugar la tabla de frecuencias (la guardamos primero en un objeto para utilizarlo después):

freqs <- table(datos$estado.salud, datos$alcohol)
freqs
##            
##             diario semanal mensual anual   ex nunca
##   Muy bueno    535    1091     952   502  331   796
##   Bueno       2062    2625    2061  1365 1084  1954
##   Regular      878     728     695   623  889  1364
##   Malo         249     135     172   157  453   508
##   Muy malo      60      29      49    48  177   222

El paquete gplots permite visualizar gráficamente la tabla de frecuencias:

library(gplots)
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
balloonplot(freqs, label = FALSE, show.margins = FALSE,
            main = "Consumo de alcohol vs. Estado de salud")

Para contrastar la relación entre ambos atributos, realizamos el test de la chi-cuadrado. De nuevo guardamos el objeto para acceder posteriormente a él:

ct1 <- chisq.test(freqs)
ct1
## 
##  Pearson's Chi-squared test
## 
## data:  freqs
## X-squared = 1640.9, df = 20, p-value < 2.2e-16

Para que el contraste sea potente debemos tener un número alto de observaciones, típicamente más de 30, y que los valores esperados sean más de cinco en cada celda, lo que se cumple sobradamente:

ct1$expected
##            
##                 diario   semanal   mensual      anual         ex     nunca
##   Muy bueno  698.39817  850.4807  725.1603  497.40568  541.51698  894.0383
##   Bueno     1851.16188 2254.2690 1922.0970 1318.41471 1435.33535 2369.7220
##   Regular    859.42652 1046.5744  892.3591  612.09156  666.37352 1100.1750
##   Malo       277.89839  338.4133  288.5472  197.92182  215.47407  355.7452
##   Muy malo    97.11503  118.2627  100.8364   69.16623   75.30008  124.3196

De igual modo podemos repetir el análisis para el resto de variables cualitativas que teníamos, o realizar un bucle para simplemente ver el p-valor y comprobar en cuáles hay relación o no:

for (atributo in c("convivencia", "actividad", "fuma", "alcohol")){}
sapply(datos[, c(1, 3:5)], function(x){
  atributo <- factor(x)
  chisq.test(table(datos$estado.salud, atributo))$p.value
})
##  convivencia    actividad         fuma      alcohol 
## 1.122581e-25 0.000000e+00 9.470356e-33 0.000000e+00

Lo que nos indica que todos están relacionados con el estado de salud. Cuando el tamaño de la muestra es tan grande, suele ser habitual que siempre se encuentre asociación.

Volver a la presentación

Análisis de correspondencias simple

La primera opción que tenemos para realizar el análisis de correspondencias con R es la función ca del paquete ca. No obstante, ya que el paquete FactoMineR tiene más funcionalidad, vamos a utilizar las funciones de este paquete directamente. Cargamos el paquete y realizamos el análisis con la función CA:

library(FactoMineR)
corres1 <- CA(freqs)

Automáticamente obtenemos el gráfico con la representación de las dos primeras dimensiones. Vemos que la primera dimensión ya explica casi toda la relación entre los atributos. A primera vista, la buena percepción del propio estado de salud se relaciona con algún consumo de alcohol, regular con los que nunca bebieron, y malo con los ex bebedores.

Como hemos guardado el objeto, podemos obtener los resultados numéricos con la función genérica summary:

summary(corres1)
## 
## Call:
## CA(X = freqs) 
## 
## The chi square of independence between the two variables is equal to 1640.89 (p-value =  0 ).
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4
## Variance               0.066   0.005   0.001   0.000
## % of var.             92.158   6.662   1.104   0.076
## Cumulative % of var.  92.158  98.820  99.924 100.000
## 
## Rows
##             Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## Muy bueno |    11.838 | -0.218 13.205  0.740 | -0.129 64.034  0.259 |
## Bueno     |    11.213 | -0.144 15.214  0.900 |  0.046 21.907  0.094 |
## Regular   |    12.234 |  0.227 17.571  0.953 |  0.026  3.136  0.012 |
## Malo      |    22.278 |  0.547 33.096  0.986 | -0.015  0.343  0.001 |
## Muy malo  |    14.424 |  0.735 20.914  0.962 | -0.141 10.581  0.035 |
##            Dim.3    ctr   cos2  
## Muy bueno -0.005  0.567  0.000 |
## Bueno      0.012  8.299  0.006 |
## Regular   -0.043 53.450  0.035 |
## Malo       0.063 36.827  0.013 |
## Muy malo   0.016  0.857  0.000 |
## 
## Columns
##             Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## diario    |     3.502 | -0.034  0.297  0.056 |  0.141 68.919  0.944 |
## semanal   |    18.233 | -0.299 27.212  0.990 | -0.020  1.720  0.005 |
## mensual   |     8.703 | -0.213 11.828  0.902 | -0.069 17.283  0.095 |
## anual     |     0.738 | -0.062  0.689  0.619 |  0.029  2.091  0.136 |
## ex        |    28.139 |  0.465 41.955  0.989 | -0.006  0.102  0.000 |
## nunca     |    12.673 |  0.237 18.019  0.943 | -0.047  9.885  0.037 |
##            Dim.3    ctr   cos2  
## diario     0.000  0.001  0.000 |
## semanal    0.020 10.169  0.004 |
## mensual    0.007  1.202  0.001 |
## anual     -0.038 21.394  0.230 |
## ex         0.048 37.572  0.011 |
## nunca     -0.033 29.662  0.019 |

Vemos que automáticamente nos ha proporcionado el resultado del contraste chi-cuadrado, por lo que podemos evitar realizarlo previamente. La primera tabla nos porporciona los autovalores de cada dimensión. El número de dimensiones es el menor número de categorías menos uno (en este caso hay 5 categorías fila y 6 categorías columna, por lo que el número de dimensiones es 5-1=4). Aquí nos fijamos en el porcentaje de varianza acumulado, y vemos que incluso con una dimensión sería suficiente. Después tenemos, para cada atributo (fila y columna) las inercias de cada nivel de atributo, y para cada dimensión su coordenada y dos medidas más: la contribución a esa coordenada concreta y la medida de calidad cos2.

Vamos a ver algunas visualizaciones interesantes para interpretar los resultados numéricos. Utilizamos algunas funciones del paquete factoextra. La primera de ellas nos sirve para seleccionar el número de dimensiones que necesitamos para explicar la mayor parte posible de la variabilidad. Es un gráfico de sedimentación, al que en este caso añadimos un criterio utilizado para esta selección (máximo entre los inversos del número de categorías menos uno, en este caso 25%).

max.porc <- max(1/(dim(freqs)-1)*100)
library(factoextra)
## Loading required package: ggplot2
## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
fviz_screeplot(corres1) +
 geom_hline(yintercept = max.porc, linetype = 2, color = "red") + 
  labs(title = "Gráfico de sedimentación", x = "Dimensiones", 
       y = "Porcentaje de variabilidad explicada")

Si necesitáramos más de dos dimensiones, deberiamos plantearnos si no podemos reagrupar los factores, para lo cual podemos visualizarlos con alguna de las medidas disponibles, y así detectar afinidades

fviz_ca_row(corres1, col.row = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE)

fviz_ca_col(corres1, col.col = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE)

En nuestro caso, como ya sabíamos, no es necesario. Ahora podemos ver para cada atributo cuál es la contribución de sus categorías a las dos dimensiones principales:

fviz_contrib(corres1, choice = "row", axes = 1:2)