Análisis multivariante. Ordenación II

Análisis de correspondencias (CA)

Técnica estadística que, al igual que el PCA, trata de resumir gran cantidad de información en un número reducido de dimensiones.

Está diseñada para analizar tablas de dos vías (CA) o de múltiples vías (MCA), y nos permite explorar la estructura de variables categóricas. Concretamente permite visualizar y analizar patrones de asociación entre modalidades de dos o más variables categóricas.

Un ejemplo sencillo

Hacer un AC con los datos de smoke e interpretar los resultados

El fichero smoke.dat es una matriz de contingencia correspondiente a 193 individuos pertenecientes a las categorias profesionales: Junior Manager (JM), Senior Manager (SM), Junior Employed (JE), Senior Employed (SE) y Secretaries (SC) y su grado de tabaquismo.

Se busca saber si existe alguna relación entre el puesto que desempeñan y el grado de tabaquismo

smk <- read.table( "http://ares.inf.um.es/00Rteam/pub/clas/data/smoke.dat" )
smk

   none light medium heavy
SM    4     2      3     2
JM    4     3      7     4
SE   25    10     12     4
JE   18    24     33    13
SC   10     6      7     2

Lo primero es entender la tabla de contingencia

¿Cuántas variables tenemos?
¿Cuántas observaciones?

Aunque la tabla de contigencia pueda confundirnos un poco si no estamos acostumbrados a ella, lo que nos está indicando es el número de casos (la frecuencia de ocurrencia) de pares de categorías.

Así para responder a la pregunta de cuántas variables tenemos, la respuesta será 2 variables: categoría profecional (con 5 categorías) y la varible tabaquismo (con 4 categorías).

Respecto al número de observaciones, obviamente será la suma total de las frecuencias en la tabla de contingencia, en este ejemplo 193.

Tabla con frecuencias marginales

dfsmk <- as.data.frame(smk)
rows <- apply(smk, 1, sum)
cols <- apply(smk, 2, sum)
dfsmk$TotRows <- rows
dfsmk <- rbind( dfsmk, Totcols = c( cols, sum( smk ) ) )
knitr::kable( round( dfsmk, 2 ), "html" ) %>% 
    kable_styling(full_width = FALSE ) %>% 
    row_spec( 6, color = "black", background = "#6a96b5", bold = TRUE ) %>% 
    column_spec( 6, color = "black", background = "#6a96b5", bold = TRUE ) %>% 
    column_spec(1, color = "black", bold = FALSE)

	none	light	medium	heavy	TotRows
SM	4	2	3	2	11
JM	4	3	7	4	18
SE	25	10	12	4	51
JE	18	24	33	13	88
SC	10	6	7	2	25
Totcols	61	45	62	25	193

Tabla de distribuciones condicionadas por filas

dfsmk <- as.matrix( smk )
dfsmk <- prop.table( dfsmk, 1 )
sumsrows <- apply( dfsmk, 1, sum )
dfsmk <- as.data.frame( dfsmk )
dfsmk$TotFreq <- sumsrows
knitr::kable( round( dfsmk, 2 ), "html" ) %>% 
    kable_styling(full_width = FALSE ) %>% 
    column_spec( 6, color = "black", background = "#6a96b5", bold = TRUE )

	none	light	medium	heavy	TotFreq
SM	0.36	0.18	0.27	0.18	1
JM	0.22	0.17	0.39	0.22	1
SE	0.49	0.20	0.24	0.08	1
JE	0.20	0.27	0.38	0.15	1
SC	0.40	0.24	0.28	0.08	1

Tabla de distribuciones condicionadas por columnas

dfsmk <- as.matrix( smk )
dfsmk <- prop.table( dfsmk, 2 )
sumscols <- apply( dfsmk, 2, sum )
dfsmk <- as.data.frame( dfsmk )
dfsmk <- rbind( dfsmk, TotFreq = c( sumscols ) )
kable( round( dfsmk, 2 ), "html" ) %>% 
    kable_styling(full_width = FALSE ) %>% 
    row_spec( 6, color = "black", background = "#6a96b5", bold = TRUE )

	none	light	medium	heavy
SM	0.07	0.04	0.05	0.08
JM	0.07	0.07	0.11	0.16
SE	0.41	0.22	0.19	0.16
JE	0.30	0.53	0.53	0.52
SC	0.16	0.13	0.11	0.08
TotFreq	1.00	1.00	1.00	1.00

smkCA <- CA( smk, graph = FALSE)
summary( smkCA )


Call:
CA(X = smk, graph = FALSE) 

The chi square of independence between the two variables is equal to 16.44164 (p-value =  0.1718348 ).

Eigenvalues
                       Dim.1   Dim.2   Dim.3
Variance               0.075   0.010   0.000
% of var.             87.756  11.759   0.485
Cumulative % of var.  87.756  99.515 100.000

Rows
         Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr
SM     |     2.673 | -0.066  0.330  0.092 |  0.194 21.356
JM     |    11.881 |  0.259  8.366  0.526 |  0.243 55.115
SE     |    38.314 | -0.381 51.201  0.999 |  0.011  0.300
JE     |    26.269 |  0.233 33.097  0.942 | -0.058 15.177
SC     |     6.053 | -0.201  7.006  0.865 | -0.079  8.052
         cos2    Dim.3    ctr   cos2  
SM      0.800 |  0.071 69.433  0.107 |
JM      0.465 | -0.034 25.619  0.009 |
SE      0.001 | -0.005  1.698  0.000 |
JE      0.058 |  0.003  1.205  0.000 |
SC      0.133 | -0.008  2.045  0.001 |

Columns
         Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr
none   |    49.186 | -0.393 65.400  0.994 |  0.030  2.934
light  |     7.059 |  0.099  3.085  0.327 | -0.141 46.317
medium |    12.610 |  0.196 16.562  0.982 | -0.007  0.174
heavy  |    16.335 |  0.294 14.954  0.684 |  0.198 50.575
         cos2    Dim.3    ctr   cos2  
none    0.006 | -0.001  0.061  0.000 |
light   0.657 |  0.022 27.282  0.016 |
medium  0.001 | -0.026 51.140  0.017 |
heavy   0.310 |  0.026 21.517  0.005 |

¿Qué nos dice la Chi cuadrado sobre las variables?
¿Son las variables independientes? ¿Cuándo lo son?
¿Cuánta varianza queda explicada por cada dimensión?

Al realizar un análisis de correspondencias (CA) se realiza un test de \(\tilde{\chi}^2\) cuya hipótesis nula (\(H_0\)) es que las variables son independientes, por tanto no están relacionadas. Sería algo equivalente a la comprobación de independencia de variables que hacíamos en PCA mediante correlaciones.

Si obtenemos un p-valor por encima de nuestro \(\alpha\) fijado (p.e \(\alpha=0.05\)) entonces no podemos rechazar la hipótesis nula y concluiríamos la existencia de independencia entre las variables. En nuestro ejemplo equivaldría a decir que no existe relación entre la categoría profesional y el grado de tabaquismo.

Sin embargo con un p valor por debajo del \(\alpha\) fijado, rechazaríamos la hipótesis nula, concluyendo, por tanto, la existencia de dependencia entre las variables.

Al igual que ocurría en PCA, aquí también obtenemos una tabla con la varianza asociada a cada dimensión (recordemos que dimensión, componente, factor son equivalentes en este contexto), varianza acumulada, etc.

Podemos visualizar las categorías de cada variable de forma independiente, o las dos variables de manera conjunta, de forma que resulte más sencilla su interpretación.

fviz_ca_col(smkCA, alpha.col = "contrib" )

fviz_ca_row(smkCA, alpha.row = "contrib" )

fviz_ca_biplot( smkCA )

Otros parámetros de interés:

Las inercias que sería el equivalente a las varianzas de las variables en un PCA. Nos indica cuánta información tiene cada valor de cada variable.
Contribución relativa (ctr), nos indica en % cuánto contribuye cada variable en cada dimensión.
\(Coseno^2\) (cos2) es la distancia de la proyección de cada variable original sobre los nuevos ejes. Es otra forma de expresar la distancia de cada variable sobre cada nueva dimensión.
Las coordenadas de las variables sobre los nuevos ejes

Interpretación:

La primera dimensión parece ser un gradiente de tabaquismo, principalmente distingue entre no fumadores (parte negativa del eje de abscisas) y fumadores en sus diferentes grados (parte positiva).
La segunda dimensión (eje de ordenadas) parece distinguir entre los diferentes puestos de trabajo, situando en la parte superior los jefes y en la parte media-inferior empleados y secretariado.

Se puede observar que estos “gradientes” de tabaquismo o puesto de trabajo, no son lineales, sino que forman una especie de parábola, que es típica de los gráficos de CA en los que existe una clara dependencia entre variables.

Un ejemplo con variables ilustrativas

Los datos corresponden a una encuesta realizada a 1724 mujeres francesas en 1974.

Las preguntas:

¿Cómo piensa vd. que es la familia perfecta?
- Mujer y hombre trabajan
- Hombre trabaja más que mujer
- Sólo hombre trabaja
- ¿Qué actividad es la mejor para una madre cuando los hijo están en el colegio?
  - Estar en casa
  - trabajar a tiempo parcial
  - Trabajar a tiempo completo
- ¿Qué piensa de la siguente afirmación?: Las mujeres que no trabajan se sienten excluidas del mundo
  - Totalmente de acuerdo
  - Bastante de acuerdo
  - Bastante en desacuerdo
  - Totalmente en desacuerdo

par( cex = 0.7)
women_work <- read.table( "http://ares.inf.um.es/00Rteam/pub/clas/data/women_work.dat", header =TRUE, 
                         row.names = 1, sep = "" )

# View( women_work )

res.2var <- CA( women_work[,1:3], graph = FALSE )
summary( res.2var )


Call:
CA(X = women_work[, 1:3], graph = FALSE) 

The chi square of independence between the two variables is equal to 233.4304 (p-value =  2.410248e-49 ).

Eigenvalues
                       Dim.1   Dim.2
Variance               0.117   0.019
% of var.             86.292  13.708
Cumulative % of var.  86.292 100.000

Rows
                          Iner*1000    Dim.1    ctr   cos2
both.man.and.woman.work |    55.487 | -0.559 40.432  0.851
man.morks.more          |    28.675 | -0.244 16.371  0.667
only.man.works          |    51.239 |  0.310 43.197  0.985
                           Dim.2    ctr   cos2  
both.man.and.woman.work |  0.233 44.429  0.149 |
man.morks.more          | -0.172 51.436  0.333 |
only.man.works          |  0.038  4.135  0.015 |

Columns
                          Iner*1000    Dim.1    ctr   cos2
stay.at.home            |    68.489 |  0.618 53.913  0.920
part.time.work          |     6.478 | -0.004  0.007  0.001
full.time.work          |    60.434 | -0.541 46.079  0.891
                           Dim.2    ctr   cos2  
stay.at.home            |  0.183 29.613  0.080 |
part.time.work          | -0.100 34.853  0.999 |
full.time.work          |  0.189 35.533  0.109 |

fviz_ca( res.2var )

Con variable ilustrativa

# Con variable ilustrativa
res.total <- CA( women_work, col.sup = 4:ncol( women_work ), graph = FALSE )
fviz_ca_biplot( res.total )

Interpretación:

La primera dimensión estaría relacionada principalmente con la disposición de la mujer a estar en su casa (izquierda trabajar full time, derecha en casa full time, en el centro parcialmente)
La segunda dimensión que sólo retiene un 13.7% de la información podría relacionarse con la postura más o menos radical de las opiniones. Situándose en la parte positiva las posturas más definidas (estar en casa - trabajar a tiempo completo, trabajar hombres y mujeres por igual - trabajar sólo los hombres) frente a la parte negativa del eje con posturas intermedias (probablemente relacionado con mujeres que les gustaría conciliar la vida familiar y laboral).

Análisis factorial exploratorio (FA)

En disciplinas como la psicología o sociología, en ocasiones no es posible medir directamente variables de interés, sino que es necesario recoger medidas indirectas que estén relacionadas con esas variables “latentes” de interés.

El AF parte de la idea de que las fuertes correlaciones entre un conjunto de variables observadas es debida a la existencia de unos pocos factores “latentes” de naturaleza más abstracta.

Desde un punto de vista conceptual, el PCA trata de encontrar factores que expliquen la mayor parte de la varianza total. Mientras que el AF busca factores que expliquen la mayor parte de la varianza común o covarianza.

Tiene en común con PCA que en ambos casos se busca nuevas variables —factores– que expliquen los datos. Y al igual que en el caso de PCA las variables no pueden estar incorreladas, porque si así fuese, no habría covarianza que explicar.

Algunos conceptos:

Varianza común:: parte de la varianza compartida por todas las variables.
Varianza única:: parte de la varianza de una variable que es propia de esa variable.
Carga factorial:: correlación entre un factor \(F_j\) y la una variable \(x_j\)

La función factanal() ejecuta un FA por el método de máxima verosimilitud sobre una matriz de covariazas o matriz de datos directamente.

El fichero mundodes.dat contiene las observaciones correspondientes a 91 países. Las variables medidas son tasa de natalidad por cada 1000 habitantes (TN1000), tasa de mortalidad por cada 1000 habitantes (TM1000), mortalidad infantil por debajo de un año (MortInfan), esperanza de vida de hombre (EspHom), esperanza de vida de mujeres (EspMuj), y producto nacional bruto per cápita (PNBpCapita).

mundodes <-  mundodes <- read.csv2("http://gauss.inf.um.es/datos/mundodes.csv", row.names=1, header = TRUE, sep = "\t")
mundoFA <- factanal( mundodes, 2, scores = "regression" )
round( mundoFA$loadings, 2 )


Loadings:
           Factor1 Factor2
TN1000     -0.87    0.28  
TM1000     -0.27    0.96  
MortInf    -0.83    0.47  
EspHom      0.82   -0.55  
EspMuj      0.86   -0.50  
PNBpCapita  0.68   -0.12  

               Factor1 Factor2
SS loadings      3.393   1.788
Proportion Var   0.566   0.298
Cumulative Var   0.566   0.863

Datos del análisis

Unicidad (uniquenesses): Porcentaje de varianza de esa variable que no ha sido explicada por los factores.
Dicho de otra manera, es la proporción de información que puede ser explicada únicamente por la propia variable.
Cargas factoriales (loadings): Cargas factoriales, son las correlaciones entre variables y factores
SS loadings: Es la suma de cuadrados de las cargas factoriales
—suma de los cuadrados de las columnas de la matriz de cargas factoriales —

round( apply( mundoFA$loadings, 2, function( x ){ sum( x^2 ) } ), 2 )

Factor1 Factor2 
   3.42    1.78

Comunalidades: Es el porcentaje de la varianza de la variable que es explicada por los factores.
Se puede calcular como \(1 - unicidades\) o como la suma de los cuadrados de las filas de la matriz de cargas factoriales.

round( ( 1 - mundoFA$uniquenesses ) , 2 )

    TN1000     TM1000    MortInf     EspHom     EspMuj 
      0.84       1.00       0.91       0.97       1.00 
PNBpCapita 
      0.48

round( apply( mundoFA$loadings, 1, function( x ){ sum( x^2 ) } ), 2 )

    TN1000     TM1000    MortInf     EspHom     EspMuj 
      0.84       1.00       0.91       0.97       1.00 
PNBpCapita 
      0.48

Una cuestión a tener en cuenta en el FA es el número de factores necesarios para explicar suficientemente bien nuestros datos. Aunque hay muchos métodos, vamos a ver uno muy sencillo que consiste simplemente en probar con distinto número de factores y factanal nos indica el p-valor de la Chi-cuadrado aplicada bajo la hipótesis de que el número de factores aplicados es suficiente. Vamos a verlo.

fa <- factanal( mundodes, factors = 1, scores = "regression" ) 
fa


Call:
factanal(x = mundodes, factors = 1, scores = "regression")

Uniquenesses:
    TN1000     TM1000    MortInf     EspHom     EspMuj 
     0.201      0.483      0.087      0.031      0.005 
PNBpCapita 
     0.577 

Loadings:
           Factor1
TN1000     -0.894 
TM1000     -0.719 
MortInf    -0.956 
EspHom      0.984 
EspMuj      0.998 
PNBpCapita  0.651 

               Factor1
SS loadings      4.617
Proportion Var   0.769

Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 45.61 on 9 degrees of freedom.
The p-value is 7.12e-07

Vemos que con p-valor de 7.12e-07 \(< 0.05\), se rechazaría la hipótesis nula de que 1 factor es suficiente.

Probaremos ahora con 2 factores

fa <- factanal( mundodes, factors = 2, scores = "regression" ) 
fa


Call:
factanal(x = mundodes, factors = 2, scores = "regression")

Uniquenesses:
    TN1000     TM1000    MortInf     EspHom     EspMuj 
     0.158      0.005      0.087      0.026      0.005 
PNBpCapita 
     0.517 

Loadings:
           Factor1 Factor2
TN1000     -0.875   0.277 
TM1000     -0.274   0.959 
MortInf    -0.833   0.468 
EspHom      0.818  -0.552 
EspMuj      0.864  -0.498 
PNBpCapita  0.684  -0.120 

               Factor1 Factor2
SS loadings      3.419   1.783
Proportion Var   0.570   0.297
Cumulative Var   0.570   0.867

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 3.99 on 4 degrees of freedom.
The p-value is 0.408

Ahora el p-valor es de 4.08e-01 \(> 0.05\), por lo que no podemos rechazar la hipótesis nula de que dos factores son suficientes. Ese sería el número mínimo de factores, para este ejemplo, necesarios para explicar nuestros datos.

Representación gráfica e interpretación

Podemos generar un gráfico biplot similar al que teníamos en PCA, donde las relaciones de las variables con los factores son representadas como vectores sobre el plano factorial y las observaciones como puntos.

plot( mundoFA$loadings[,1:2], type = "n", xlim = c( -2, 3 ), ylim = c( -2, 3 ) )
abline( h = 0, v = 0 )
text( mundoFA$loadings[,1:2], labels = names( mundodes ), 
      cex = 1, col = "blue", pos = c( 3,4,3,2,3,3 ) )
arrows( 0, 0, mundoFA$loadings[,1], mundoFA$loadings[,2], 
        col = "red", angle = 7, lwd=2 )
par( new = TRUE )
plot( mundoFA$scores[,1:2], type = "p", xlab = "", 
     ylab = "", xaxt = "n", yaxt = "n", cex = 1, pch=20)

El primer eje (Factor 1) podría interpretarse como una medida del grado de desarrollo de un país. Con peso positivo hay variables como PNB, Esp. de vida de mujeres y de hombres (valores altos en paises desarrollados), mientras que con peso negativo están principalmente la Mort.Infantil y la Tasa de Natalidad, que son bajos en estos mismos paises. Por otro lado, el segundo eje, está fuertemente asociado a la Tasa de Mortalidad, separa a los paises con alta tasa de los de baja tasa.

Comprobando la interpretación

abline( h = 0, v = 0 )
text( mundoFA$loadings[,1:2], labels = names( mundodes ), 
      cex = 1, col = "blue", pos = c( 3,4,3,2,3,3 ) )
arrows( 0, 0, mundoFA$loadings[,1], mundoFA$loadings[,2], 
        col = "red", angle = 7, lwd=2 )
par( new = TRUE )
plot( mundoFA$scores[,1:2], type = "p", xlab = "", 
     ylab = "", xaxt = "n", yaxt = "n", cex = 1, pch=20)
Cng <- mundoFA$scores[ rownames( mundoFA$scores ) == "Congo", c( 1, 2 ) ]
Swd <- mundoFA$scores[ rownames( mundoFA$scores ) == "Sweden", c( 1, 2 ) ]
puntos <- as.matrix( rbind( Cng, Swd ) )
points( puntos[ , 1 ], puntos[ , 2 ], col = c("red", "green" ) , pch = c( 20, 20 ) , cex = 2 )
text( puntos[ , 1], puntos[ , 2] , labels = c( "Cng", "Swd" ), pos = 4 )

mundodes[ ( rownames( mundodes ) == "Congo" | rownames( mundodes ) == "Sweden" ),  ]

       TN1000 TM1000 MortInf EspHom EspMuj PNBpCapita
Sweden   14.5   11.1     5.6   74.2   80.0      23660
Congo    46.1   14.6    73.0   50.1   55.3       1010

¿ Dónde estaríamos nosotros?

abline( h = 0, v = 0 )
text( mundoFA$loadings[,1:2], labels = names( mundodes ), 
      cex = 1, col = "blue", pos = c( 3,4,3,2,3,3 ) )
arrows( 0, 0, mundoFA$loadings[,1], mundoFA$loadings[,2], 
        col = "red", angle = 7, lwd=2 )
par( new = TRUE )
plot( mundoFA$scores[,1:2], type = "p", xlab = "", 
     ylab = "", xaxt = "n", yaxt = "n", cex = 1, pch=20)
Spn <- mundoFA$scores[ rownames( mundoFA$scores ) == "Spain", c( 1, 2 ) ]
points( Spn[ 1 ], Spn[ 2 ], col = c("blue" ) , pch = c( 20 ) , cex = 3 )

par( cex = 0.8 )
biplot( jitter(mundoFA$scores, amount = 0.5), mundoFA$loadings )
abline( h = 0, v = 0 )

Biplot coloreados

Colorear los gráficos según valores, suele ayudar en muchos casos a mejorar la interpretación o a corroborar que nuestra interpretación es correcta.

En este gráfico los paises han sido coloreados según su PIB, desde el rojo —PIB más alto— al azul —más bajo—

mundoOrder <- mundodes[ order( -mundodes$PNBpCapita ),  ]
redblue <- colorRampPalette( c( "red", "orange", "yellow", "green", "blue" ), interpolate = "linear" )
colors <- redblue( nrow(mundoOrder) )
mundoFA <- factanal( mundoOrder, 2, scores = "regression" )
plot( mundoFA$loadings[,1:2], type = "n", xlim = c( -2, 3 ), ylim = c( -2, 3 ) )
abline( h = 0, v = 0 )
text( mundoFA$loadings[,1:2], labels = names( mundodes ), 
      cex = 1, col = "blue", pos = c( 3,4,3,2,3,3 ) )
arrows( 0, 0, mundoFA$loadings[,1], mundoFA$loadings[,2], 
        col = "red", angle = 7, lwd=2 )
par( new = TRUE )
plot( mundoFA$scores[,1:2], type = "p", xlab = "", 
     ylab = "", xaxt = "n", yaxt = "n", cex = 1, pch=20, col = colors)

En este caso se se ha coloreado por la tasa de mortalidad —TM1000—

mundoOrder <- mundodes[ order( mundodes$TM1000 ),  ]
redblue <- colorRampPalette( c( "red", "orange", "yellow", "green", "blue" ), interpolate = "linear" )
colors <- redblue( nrow(mundoOrder) )
mundoFA <- factanal( mundoOrder, 2, scores = "regression" )
plot( mundoFA$loadings[,1:2], type = "n", xlim = c( -2, 3 ), ylim = c( -2, 3 ) )
abline( h = 0, v = 0 )
text( mundoFA$loadings[,1:2], labels = names( mundodes ), 
      cex = 1, col = "blue", pos = c( 3,4,3,2,3,3 ) )
arrows( 0, 0, mundoFA$loadings[,1], mundoFA$loadings[,2], 
        col = "red", angle = 7, lwd=2 )
par( new = TRUE )
plot( mundoFA$scores[,1:2], type = "p", xlab = "", 
     ylab = "", xaxt = "n", yaxt = "n", cex = 1, pch=20, col = colors)

Manos a la obra

El reto de este módulo es hacer y comparar varios análisis de correspondencias muy sencillos para ir practicando, asi que… manos a la obra.

Referencias

Nuevamente aparece como referencia el archiconocido paquete FactoMineR.
El mismo libro de módulo anterior, en su capítulo 2 trata el análisis de correspondencias. (visitar).
Breve tutorial introductorio al análisis factorial exploratorio (visitar).
Completo tutorial sobre correspondencias, para extraer todo su jugo (visitar).