Técnicas de ordenación II: CA, FA
Curso 013: Introducción al análisis de datos multivariantes
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?
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?
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-inferor empleados y secretariado.
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)