Análisis de correspondencias

Definición general

El Análisis de correspondencias (CA) es un método utilizado para establecer posibles asociaciones entre variables categóricas: entre las clases de cada una de ellas. El fin es establecer patrones o estructuras en los datos, mediante su observación.

Es un método de análisis más bien exploratorio, que se complementa con análisis de tipo inferencial, tales como el Analisis de Modelos Log-lineales o el Análisis de Regresión Logísitca.

Después de rechazar la hipótesis \(H_0\) de independéncia de carácteres con el test de \(\chi^2\), este test consiste en identificar, visualizar y analizar las magnitudes obtenidas entre variables más próximas, por lo tanto más relacionadas.

Consiste en descomponer el estadístico \(\chi^2\) de la tabla de doble entrada

\[ \begin{array}{cc|cccc|c} & B & 1 & 2 & ... & b & \\ A \\ \hline 1 & & n_{11} & n_{12} & ... & n_{1b} & n_1.\\ 2 & & n_{21} & n_{22} & ... & n_{2b} & n_2.\\ ...& & ...&...&....&.... &...\\ a & & n_{a1} & n_{a2} &...&n_{ab} & n_a.\\ \hline & & n._1 & n._2 &... & n._b & n \\ \end{array} \]

en distintas componentes, correspondientes a diferentes dimensiones, de la heterogenidad entre filas y columnas asignándoles una escala a las filas y otra a las columnas que maximice la correlación entre los pares de variables resultantes.

Definición EAAR: En otras palabra, el \(CA\) representa datos cualitativos multivariantes, determinando coordenadas que permitan representar las clases (o “valores”, según EAAR) de las variables, de manera que las clases (“valores”) más cercanas sean las más relacionadas.

Análisis directo de la contribución a \(\lambda\) de cada combinación

Una primera aproximación para sería utilizar la contribución al estadístico \(\lambda\)

\[ \lambda=\sum_{i=1}^a\sum_{j=1}^b\frac{(n_{ij}-n_i. n._j/n)^2}{n_i.n._j/n} \]

utilizando la función \(\texttt{chisq.test()}\), como en el siguiente ejemplo:

Ejemplo 3.2 EAAR

Tabla de doble entrada de Nanny Wermuth (1973) de \(n=6851\) nacimientos, con dos variables:

• Características de la madre:
  • Madre joven no fumadora (\(\texttt{jnf}\))
  • Madre joven fumadora (\(\texttt{jf}\))
  • Madre mayor no fumadora (\(\texttt{mnf}\))
  • Madre mayor fumadora (\(\texttt{mf}\))
• Estado del recien nacido:
  • Prematuro, murió (\(\texttt{pm}\))
  • Prematuro, sobrevivió (\(\texttt{pv}\))
  • Gestación completa, murió (\(\texttt{gcm}\))
  • Gestación completa, sobrevivió (\(\texttt{gcv}\))

X <- matrix(c(50,9,41,4,315,40,147,11,24,6,14,1,4012,459,1594,124),ncol = 4)
colnames(X) <- c('pm','pv', 'gcm','gcv')
rownames(X) <- c('jnf','jf','mnf','mf')

-	pm	pv	gcm	gcv	ni.
jnf	50	315	24	4,012	4,401
jf	9	40	6	459	514
mnf	41	147	14	1,594	1,796
mf	4	11	1	124	140
n.j	104	513	45	6,189	6,851

Si ejecutamos el test de \(\chi^2\),

chisq.test(X)

## 
##  Pearson's Chi-squared test
## 
## data:  X
## X-squared = 19.109, df = 9, p-value = 0.02428

Veremos que, con \(\lambda=\) 19.109 y un \(\text{p-value}=\) 0.024, podemos inferir que ambas variables están relcionadas.

Siguiendo utilizando la función \(\texttt{chisq.test()}\), podemos aprovechar de varios elementos de clase \(\texttt{htest}\) para poder obtener una table de valores observados:

chisq.test(X)$observed

##     pm  pv gcm  gcv
## jnf 50 315  24 4012
## jf   9  40   6  459
## mnf 41 147  14 1594
## mf   4  11   1  124

de valores esperados:

chisq.test(X)$expected

##            pm        pv        gcm       gcv
## jnf 66.808349 329.54503 28.9074588 3975.7392
## jf   7.802657  38.48810  3.3761495  464.3331
## mnf 27.263757 134.48373 11.7968180 1622.4557
## mf   2.125237  10.48314  0.9195738  126.4720

Matriz de discrepancias

y podemos descomponer el estadístico \(\lambda\) de Pearson en nuestra tabla, directamente con el componente \(\texttt{residuals}\) obteniendo directamente la matriz de discrepacias

chisq.test(X)$residuals

##             pm         pv         gcm        gcv
## jnf -2.0564099 -0.8012301 -0.91274971  0.5750808
## jf   0.4286447  0.2437018  1.42800017 -0.2474937
## mnf  2.6307229  1.0792952  0.64145758 -0.7064523
## mf   1.2860049  0.1596343  0.08386956 -0.2198162

o con el cálculo manual

(chisq.test(X)$observed - chisq.test(X)$expected) / sqrt(chisq.test(X)$expected)

##             pm         pv         gcm        gcv
## jnf -2.0564099 -0.8012301 -0.91274971  0.5750808
## jf   0.4286447  0.2437018  1.42800017 -0.2474937
## mnf  2.6307229  1.0792952  0.64145758 -0.7064523
## mf   1.2860049  0.1596343  0.08386956 -0.2198162

Los valores de la tabla, elevados al cuadrado y sumados darán el valor del estadístico:

sum(((chisq.test(X)$observed - chisq.test(X)$expected) / sqrt(chisq.test(X)$expected))^2)

## [1] 19.10895

Análisis de correspondencias bidimensional

Siguiendo con el objetivo de determinar las coordenadas de cada clase (“valores”) de la tabla, valga la advertencia inicial que el número \(a\) de filas tiene que ser \(\leq\) al número de columnas, \(b\). Si no fuera así, habrá que invertirlas (sin embargo, ejecutando en R, no haría falta).

Con este método, podemos considerar nuestra matriz de discrepancias \(\textbf{C}\) como el producto de tres matrices: \[ \textbf{C=UDV}^t \]

Descomposición en valores singulares

Donde:

• Las columnas de \(\textbf{U}\) representan los \(a\) autovectores del producto \(\textbf{CC}^t\).
• Las columnas de \(\textbf{V}\) son los primeros \(a\) autovectores de \(\textbf{C}^t\textbf{C}\).
• \(\textbf{D}\) es una matriz diagonal, \(a\times a\), donde los elementos de la diagonal son las raices cuadradas de los autovalores de \(\textbf{CC}^t\).

Esta es la Descomposición en valores singulares. En R se obtiene fácilmente con la función \(\texttt{svd()}\), si la matriz es cuadrada. Si no, hay que hacer los productos de las matrices.

Volviendo a la matriz cuadrada del ejemplo,

C <- chisq.test(X)$residuals
C

##             pm         pv         gcm        gcv
## jnf -2.0564099 -0.8012301 -0.91274971  0.5750808
## jf   0.4286447  0.2437018  1.42800017 -0.2474937
## mnf  2.6307229  1.0792952  0.64145758 -0.7064523
## mf   1.2860049  0.1596343  0.08386956 -0.2198162

y, habiendo comprobado que es cuadrada, podremos determinar los elementos que las descomponen, con \(\texttt{svd()}\) y sin:

U

U con svd()

U <- svd(C)$u
U

##            [,1]        [,2]        [,3]      [,4]
## [1,] -0.5890434  0.09691339  0.03532335 0.8014911
## [2,]  0.2323953 -0.91664751  0.17528323 0.2739079
## [3,]  0.7150964  0.25779684 -0.40003306 0.5120073
## [4,]  0.2960701  0.28966520  0.89888909 0.1429509

U calculando el eigenvector de \(\textbf{CC}^t\)

eigen(C %*% t(C))$vectors

##            [,1]        [,2]        [,3]      [,4]
## [1,]  0.5890434 -0.09691339  0.03532335 0.8014911
## [2,] -0.2323953  0.91664751  0.17528323 0.2739079
## [3,] -0.7150964 -0.25779684 -0.40003306 0.5120073
## [4,] -0.2960701 -0.28966520  0.89888909 0.1429509

V

V con svd()

V <- svd(C)$v
V

##            [,1]        [,2]       [,3]        [,4]
## [1,]  0.8578520  0.35470453  0.3508481 -0.12320822
## [2,]  0.3235723  0.01813457 -0.9055896 -0.27364133
## [3,]  0.3248658 -0.93436886  0.1218550 -0.08104557
## [4,] -0.2320628  0.02847212  0.2048522 -0.95045872

V calculando el eigenvector de \(\textbf{C}^t\textbf{C}\)

eigen(t(C) %*% C)$vectors

##            [,1]        [,2]       [,3]       [,4]
## [1,]  0.8578520  0.35470453  0.3508481 0.12320822
## [2,]  0.3235723  0.01813457 -0.9055896 0.27364133
## [3,]  0.3248658 -0.93436886  0.1218550 0.08104557
## [4,] -0.2320628  0.02847212  0.2048522 0.95045872

D

D con svd()

D <- svd(C)$d
D

## [1] 4.164935e+00 1.292606e+00 3.023947e-01 5.329125e-16

Recordar que \(D\) es una matriz diagonal:

D <- matrix(c(D[1],0,0,0,0,D[2],0,0,0,0,D[3],0,0,0,0,D[4]), ncol = 4)
D

##          [,1]     [,2]      [,3]         [,4]
## [1,] 4.164935 0.000000 0.0000000 0.000000e+00
## [2,] 0.000000 1.292606 0.0000000 0.000000e+00
## [3,] 0.000000 0.000000 0.3023947 0.000000e+00
## [4,] 0.000000 0.000000 0.0000000 5.329125e-16

D calculando el eigenvalue de \(\textbf{CC}^t\)

sqrt(pmax(0, eigen(C %*% t(C))$values)) # pmax() para que no de NaN en el 4º elemento

## [1] 4.1649348 1.2926058 0.3023947 0.0000000

Si recuperamos los tres elementos \(\textbf{U,D y V}\) y los multiplicamos, como arriba (\(\textbf{C=UDV}^t\)), obtendremos

U %*% D %*% t(V)

##            [,1]       [,2]        [,3]       [,4]
## [1,] -2.0564099 -0.8012301 -0.91274971  0.5750808
## [2,]  0.4286447  0.2437018  1.42800017 -0.2474937
## [3,]  2.6307229  1.0792952  0.64145758 -0.7064523
## [4,]  1.2860049  0.1596343  0.08386956 -0.2198162

Que es exactamente nuestra matriz de discrepancias.

Obtención de la coordendas

\[ \begin{align} & X_{1i}=U_{1i}·\sqrt{\frac{n}{n_{i·}}}\\ & X_{2i}=U_{2i}·\sqrt{\frac{n}{n_{i·}}}\\ & Y_{j1}=V_{j1}·\sqrt{\frac{n}{n_{·j}}}\\ & Y_{j2}=V_{j2}·\sqrt{\frac{n}{n_{·j}}} \end{align} \]

Volviendo al ejemplo de arriba: podemos generar unos loops para calcular las primera y las segundas coordenadas de \(X\) y de \(Y\) de las matrices \(\textbf{U}\) y \(\textbf{V}\):

Coordenadas de X

cX <- matrix(rep(NA,8), ncol = 2) 

for (i in 1:4) {
  for (j in 1:2) {
    cX[i,j] <- U[i,j] * sqrt(sum(X)/sum(X[i,]))
  }
}
# U[2,1] * sqrt(sum(X)/sum(X[2,])) # Ejemplo

cX # Primeras y segundas coordenadas de X

##            [,1]       [,2]
## [1,] -0.7349344  0.1209164
## [2,]  0.8484433 -3.3465536
## [3,]  1.3966526  0.5035022
## [4,]  2.0711306  2.0263259

Coordenadas de Y

cY <- matrix(rep(NA,8), ncol = 2) 

for (i in 1:4) {
  for (j in 1:2) {
    cY[i,j] <- V[i,j] * sqrt(sum(X)/sum(X[,i]))
  }
}

# V[2,1] * sqrt(sum(X)/sum(X[,2])) # Ejemplo

cY # Primeras y segundas coordenadas de Y

##            [,1]         [,2]
## [1,]  6.9626198   2.87890308
## [2,]  1.1824685   0.06627130
## [3,]  4.0084341 -11.52893220
## [4,] -0.2441587   0.02995619

Con \(a\neq b\), ya no vale la igualdad \[ \textbf{C=UDV}^t \]

Sin embargo, las coordenadas se calculan exactamente igual.

Ejemplo 3.3 EAAR

Los famosos datos del color del pelo y de los ojos de los 5387 habitantes del pueblo escocés de Caithness. Datos contenidos en la

#df <- read.delim('datosr/fisher_eyes.txt')[1:4,1:6] 

#rownames(df) <- df[,1]

#m <- matrix(c(326,688,343,98,38,116,84,48,241,584,909,403,110,188,412,681,3,4,26,85),ncol = 5)
#colnames(m) <- c("Rubio", "Rojo", "Cast","Moreno", "Negro")
#rownames(m) <- c("Azules", "Claros", "Castaño", "Obscuro")

m <- MASS::caith


m %>%
  tibble::rownames_to_column("Clases") %>% 
  flextable::flextable() %>% 
  flextable::set_caption("Colores de ojos y del pelo en Caithness (Fisher, 1940)")

Colores de ojos y del pelo en Caithness (Fisher, 1940)

Clases	fair	red	medium	dark	black
blue	326	38	241	110	3
light	688	116	584	188	4
medium	343	84	909	412	26
dark	98	48	403	681	85

Si realizamos un test de independencia donde la \(H_0\) es la independiencia de las dos variables (Color del pelo, Color de los ojos),

chisq.test(m)

## 
##  Pearson's Chi-squared test
## 
## data:  m
## X-squared = 1240, df = 12, p-value < 2.2e-16

Obtenemos un \(\texttt{p-value}\) infinitesimal que exluye la independencia de las dos variables.

Mediante el Análisis de correspondencias podremos verificar si hay relación entre ciertas clases de una variable con alguna clase de la otra.

Valores observados

chisq.test(m)$observed

##        fair red medium dark black
## blue    326  38    241  110     3
## light   688 116    584  188     4
## medium  343  84    909  412    26
## dark     98  48    403  681    85

Valores esperados

chisq.test(m)$expected

##            fair      red   medium     dark    black
## blue   193.9280 38.11918 284.8275 185.3978 15.72749
## light  426.7496 83.88342 626.7793 407.9785 34.60924
## medium 479.1479 94.18303 703.7383 458.0720 38.85873
## dark   355.1745 69.81437 521.6549 339.5517 28.80453

Matriz de discrepancias

C <- chisq.test(m)$residuals
C

##              fair         red    medium       dark     black
## blue     9.483979 -0.01930262 -2.596906  -5.537407 -3.209321
## light   12.646503  3.50663997 -1.708741 -10.890844 -5.203033
## medium  -6.219798 -1.04927862  7.737531  -2.152635 -2.062785
## dark   -13.646052 -2.61077971 -5.195102  18.529854 10.470584

Lambda

sum((chisq.test(m)$residuals)^2)

## [1] 1240.039

En este caso, siendo la matriz \(4\times 5\), no podemos usar la función \(\texttt{svd()}\).

Obtendremos, por lo tanto \(\textbf{U}\) y \(\textbf{V}\), recordando que

1. Las columnas de \(\textbf{U}\) representan los \(a\) autovectores del producto \(\textbf{CC}^t\).

U <- eigen(C %*% t(C))$vectors

#rownames(U) <- rownames(C)
#colnames(U) <- rownames(C[4:1,])
U

##             [,1]       [,2]        [,3]      [,4]
## [1,] -0.32740153  0.3481491  0.79894718 0.3650806
## [2,] -0.53470247  0.2762034 -0.58694655 0.5415706
## [3,]  0.04321499 -0.8105596  0.10869354 0.5738565
## [4,]  0.77783929  0.3814407 -0.07323162 0.4940710

2. Las columnas de \(\textbf{V}\) son los primeros \(a\) autovectores de \(\textbf{C}^t\textbf{C}\).

V <- eigen(t(C) %*% C)$vectors

#rownames(U) <- rownames(C)
#colnames(U) <- rownames(C[4:1,])
V

##             [,1]       [,2]        [,3]        [,4]       [,5]
## [1,] -0.63337328  0.5208721  0.22198126  0.00000000 -0.5274987
## [2,] -0.12040878  0.0641327 -0.92784507 -0.29524104 -0.1825513
## [3,] -0.05929716 -0.7563782  0.06953154  0.04105527 -0.6464176
## [4,]  0.67018848  0.3045289  0.17534533 -0.49222197 -0.4302105
## [5,]  0.36286535  0.2444040 -0.23291035  0.81784150 -0.2923755

Coordenadas de X

cX <- matrix(rep(NA,8), ncol = 2)
rownames(cX) <- rownames(m)
colnames(cX) <- c("X1","X2")

for (i in 1:4) {
  for (j in 1:2) {
    cX[i,j] <- U[i,j] * sqrt(sum(m)/sum(m[i,]))
  }
}
# U[2,1] * sqrt(sum(X)/sum(X[2,])) # Ejemplo

cX # Primeras y segundas coordenadas de X

##                 X1         X2
## blue   -0.89679252  0.9536227
## light  -0.98731818  0.5100045
## medium  0.07530627 -1.4124778
## dark    1.57434710  0.7720361

Coordenadas de Y

cY <- matrix(rep(NA,10), ncol = 2)
rownames(cY) <- colnames(m)
colnames(cY) <- c("Y1","Y2")

for (i in 1:5) {
  for (j in 1:2) {
    cY[i,j] <- V[i,j] * sqrt(sum(m)/sum(m[,i]))
  }
}

# V[2,1] * sqrt(sum(X)/sum(X[,2])) # Ejemplo

cY # Primeras y segundas coordenadas de Y

##                 Y1         Y2
## fair   -1.21871379  1.0022432
## red    -0.52257500  0.2783364
## medium -0.09414671 -1.2009094
## dark    1.31888486  0.5992920
## black   2.45176017  1.6513565

La distancia \(\chi^2\)

Volviendo a la ecuación de \(\chi^2\) \[ \chi^2=\sum_{i=1}^a\sum_{j=1}^b\frac{(n_{ij}-n_i. n._j/n)^2}{n_i.n._j/n}\\ \]

Dividiendo numerador y denominador por los totales de las filas podremos expresar la misma ecuación en terminos relativos de perfiles: \[ \text{perfil observado}:po=n_{ij}/n_j.\\ \text{perfil esperado}:pe=(n_i. n._j/n)/n_j.\\ \chi^2=\text{total de la fila}\times \frac{(\text{perfiles observados de la fila - perfiles esperados de la fila})^2}{\text{perfiles esperados de la fila}}=\\ =\sum_{i=1}^a\sum_{j=1}^b n_{·j}·\frac{(po-pe)^2}{pe} \]

Dividamos ahora los dos lados de la ecuación por el tamaño total de la muestra \(n\):

\[ \frac{\chi^2}{n}=\phi^2=\sum_{i=1}^a\sum_{j=1}^b \frac{n_{i·}}{n}·\frac{(po-pe)^2}{pe}=\\ \sum_{i=1}^a\sum_{j=1}^b mass_i·\frac{(po-pe)^2}{pe} \]

Dado que \[\sum \sqrt\frac{(po-pe)^2}{pe}\]

es equivalente a una medida de distancia euclidea respecto a un centro que podemos denominar \(Dist\chi^2\):

\[ \text{dist-}\chi^2=\sum \sqrt\frac{(po-pe)^2}{pe} \]

Se trata de una medida importante, que visualiza la distancia de cada valor de la tabla de un hipotético centro y que, como veremos quedará reflejado en una representació gráfica.

La inercia

Para cada fila, \(\phi\), la inercia, será: \[ \phi^2=mass_i·(\text{dist-}\chi^2)^2 \]

Dado que la suma de las masas es 1, podemos decir que la inercia es la media ponderada de los cuadrados de las distancias \(χ^2\) entre los perfiles fila y su perfil media. Por tanto, la inercia será alta cuando los perfiles fila presenten grandes desviaciones con relación a su media, y será baja cuando éstos se hallen cerca de la media.

Calcular \(mass\)

mass_x <- c()

for (i in 1:nrow(m)) {
  mass_x[i] <- sum(chisq.test(m)$observed[i,])/sum(chisq.test(m)$observed)
}

mass_x

## [1] 0.1332838 0.2932987 0.3293113 0.2441062

mass_y <- c()

for (i in 1:ncol(m)) {
  mass_y[i] <- sum(chisq.test(m)$observed[,i])/sum(chisq.test(m)$observed)
}

mass_y

## [1] 0.27009467 0.05309077 0.39669575 0.25821422 0.02190459

Calcular \(dist\chi^2\)

chiDist_x <- c()

for(i in 1:nrow(m)){
  
  pox <- chisq.test(m)$observed[i,]/sum(chisq.test(m)$observed[i,])
  
  pex <- chisq.test(m)$expected[i,]/sum(chisq.test(m)$observed[i,])
  
  chiDist_x[i] <- sqrt(sum((pox-pex)^2/pex))
}

#po1 <- chisq.test(m)$observed[1,]/sum(chisq.test(m)$observed[1,])
#pe1 <- chisq.test(m)$expected[1,]/sum(chisq.test(m)$observed[1,])

chiDist_x

## [1] 0.4378549 0.4506201 0.2473594 0.7153975

chiDist_y <- c()

for(i in 1:ncol(m)){
  
  poy <- chisq.test(m)$observed[,i]/sum(chisq.test(m)$observed[,i])
  
  pey <- chisq.test(m)$expected[,i]/sum(chisq.test(m)$observed[,i])
  
  chiDist_y[i] <- sqrt(sum((poy-pey)^2/pey))
}

#po1 <- chisq.test(m)$observed[1,]/sum(chisq.test(m)$observed[1,])
#pe1 <- chisq.test(m)$expected[1,]/sum(chisq.test(m)$observed[1,])

chiDist_y

## [1] 0.5712352 0.2658543 0.2125256 0.5979011 1.1321927

Calcular la inercia

inert_x <- mass_x * (chiDist_x)^2

inert_x

## [1] 0.02555277 0.05955678 0.02014947 0.12493199

inert_y <- mass_y * (chiDist_y)^2

inert_y

## [1] 0.088134492 0.003752377 0.017917615 0.092307908 0.028078616

La inercia total

sum(inert_x)

## [1] 0.230191

sum(inert_y)

## [1] 0.230191

Notar que las inercias totales son iguales.

La inercia principal

La inercia principal es representada por los valores propios (eigenvalues) de \(\textbf{U}=\textbf{CC}^t:eigenval(\textbf{CC}^t)\)

inPr <- eigen(C %*% t(C))$values/sum(m)
inPr

## [1]  1.992448e-01  3.008677e-02  8.594814e-04 -2.402361e-18

La suma de las inercia principal coincide con la inercia total

sum(inPr)

## [1] 0.230191

Resumen

Podemos por lo tanto reunir todos los elementos en unos cuadros

Inercia principal

inPr

## [1]  1.992448e-01  3.008677e-02  8.594814e-04 -2.402361e-18

Inercia principal

inPr

## [1]  1.992448e-01  3.008677e-02  8.594814e-04 -2.402361e-18

Rows

data.frame(cbind(mass_x,chiDist_x,inert_x,cX))

##           mass_x chiDist_x    inert_x          X1         X2
## blue   0.1332838 0.4378549 0.02555277 -0.89679252  0.9536227
## light  0.2932987 0.4506201 0.05955678 -0.98731818  0.5100045
## medium 0.3293113 0.2473594 0.02014947  0.07530627 -1.4124778
## dark   0.2441062 0.7153975 0.12493199  1.57434710  0.7720361

Columns

data.frame(cbind(mass_y,chiDist_y,inert_y,cY))

##            mass_y chiDist_y     inert_y          Y1         Y2
## fair   0.27009467 0.5712352 0.088134492 -1.21871379  1.0022432
## red    0.05309077 0.2658543 0.003752377 -0.52257500  0.2783364
## medium 0.39669575 0.2125256 0.017917615 -0.09414671 -1.2009094
## dark   0.25821422 0.5979011 0.092307908  1.31888486  0.5992920
## black  0.02190459 1.1321927 0.028078616  2.45176017  1.6513565

Escalado de coordenadas con el método “symmetric”

Las coordenadas obtenidas como se ha descrito arriba se pueden escalar. El método estándar de escalado es el simétrico: utilizando la inercia principal, las coordenadas se multiplican por las componentes de la inercia principal (la Dimensión 1 por la 1a Componente, la 2º con la 2º, etc.)

\[Coord_{escala}=Coord\times\sqrt{eigenvalue}\]

r <- matrix(rep(NA,nrow(m)*2),ncol = 2)

for (i in 1:2) {
  r[,i] <- cX[,i] * sqrt(inPr[i])
}

row.names(r) <- row.names(m)
r

##               [,1]        [,2]
## blue   -0.40029985  0.16541100
## light  -0.44070764  0.08846303
## medium  0.03361434 -0.24500190
## dark    0.70273880  0.13391383

col <- matrix(rep(NA,length(m)*2),ncol = 2)

for (i in 1:2) {
  col[,i] <- cY[,i] * sqrt(inPr[i])
}

rownames(col) <- colnames(m)



col

##               [,1]        [,2]
## fair   -0.54399533  0.17384449
## red    -0.23326097  0.04827895
## medium -0.04202412 -0.20830421
## dark    0.58870853  0.10395044
## black   1.09438828  0.28643670

Con la librería \(\texttt{CA}\)

Cálculo de parámetros

library(ca)

ca(m)

## 
##  Principal inertias (eigenvalues):
##            1        2        3       
## Value      0.199245 0.030087 0.000859
## Percentage 86.56%   13.07%   0.37%   
## 
## 
##  Rows:
##             blue    light    medium      dark
## Mass    0.133284 0.293299  0.329311  0.244106
## ChiDist 0.437855 0.450620  0.247359  0.715398
## Inertia 0.025553 0.059557  0.020149  0.124932
## Dim. 1  0.896793 0.987318 -0.075306 -1.574347
## Dim. 2  0.953623 0.510004 -1.412478  0.772036
## 
## 
##  Columns:
##             fair      red    medium      dark     black
## Mass    0.270095 0.053091  0.396696  0.258214  0.021905
## ChiDist 0.571235 0.265854  0.212526  0.597901  1.132193
## Inertia 0.088134 0.003752  0.017918  0.092308  0.028079
## Dim. 1  1.218714 0.522575  0.094147 -1.318885 -2.451760
## Dim. 2  1.002243 0.278336 -1.200909  0.599292  1.651357

Representación gráfica

La misma librería dispone de la función \(\texttt{plot.ca}\), para representar fácilmente el resultado de la \(CA\), directamente a partir de la tabla de doble entrada inicial.

plot(ca(m))

Representación con \(\texttt{ggplot2}\)

Con algunas fáciles manipulaciones, podemos representar el mismo gráfico en \(\texttt{ggplot2}\). Incluso lo podemos mejorar incluyendo la información de la inercia, que será proporcional al tamaño del del símbolo:

library(ggplot2)
library(ggrepel)


col <- as.data.frame(col) %>% bind_cols(data.frame(cbind(mass_y,chiDist_y,inert_y,cY))) %>% 
  tibble::rownames_to_column() %>% 
  mutate(Variable= 'Hair') %>% 
  rename(Dim.1=V1,
         Dim.2=V2,
         mass=mass_y,
         chiDist=chiDist_y,
         inertia=inert_y,
         C1=Y1,
         C2=Y2)

r <- as.data.frame(r) %>% bind_cols(data.frame(cbind(mass_x,chiDist_x,inert_x,cX))) %>% 
  tibble::rownames_to_column() %>% 
  mutate(Variable= 'Eyes') %>% 
  rename(Dim.1=V1,
         Dim.2=V2,
         mass=mass_x,
         chiDist=chiDist_x,
         inertia=inert_x,
         C1=X1,
         C2=X2) %>% 
  bind_rows(as.data.frame(col))

ggplot(r, aes(x=Dim.1, y=Dim.2, 
              color = Variable, 
              shape = Variable,
              size = inertia)) +
  geom_point()+
  geom_text_repel(aes(label = rowname)) +
  geom_vline(xintercept= 0)+
  geom_hline(yintercept = 0)+
  xlab(paste0("Dimensión 1: ",scales::percent(inPr[1]/sum(inPr),0.1))) +
  ylab(paste0("Dimensión 2: ", scales::percent(inPr[2]/sum(inPr),0.1)))

Como vemos, el gráfico es especular respecto al que produce por defecto la función \(\texttt{ca}\). Esto porque, como ya hemos visto para la \(PCA\), los signos en las matrices que generan las componentes principales se pueden invertir, pero no tiene importancia de cara al análisis.

Método mixto usando \(\texttt{ca}\) y \(\texttt{ggplot2}\)

ca() Permite evitar una grande cantidad de operaciones previas (construcción de matrices cálculo de parámetros, etc.). Podemos por lo tanto aprovechar de dicha ventaja sin perder la claredad y eficacia de representación de ggplot() de esta forma.

a <- ca(m)

# Inertia principal = eigenval(C %*% C^t) / n

inPr <- eigen(chisq.test(m)$residuals %*% t(chisq.test(m)$residuals))$values / sum(m)

# Construcción directamente desde de la library ca

r <- rbind(data.frame(Variable="Eyes",
                 mass=a$rowmass,distChi=a$rowdist,inertia=a$rowinertia,a$rowcoord[,1:2]),
data.frame(Variable="Hair",
           mass=a$colmass,distChi=a$coldist,inertia=a$colinertia,a$colcoord[,1:2])) %>% 
  mutate(Dim1 = Dim1 * sqrt(inPr[1]),
         Dim2 = Dim2 * sqrt(inPr[2])) %>% 
  tibble::rownames_to_column()

# Gráfico idéntico

ggplot(r, aes(x=Dim1, y=Dim2, 
              color = Variable, 
              shape = Variable,
              size = inertia)) +
  geom_point()+
  geom_text_repel(aes(label = rowname)) +
  geom_vline(xintercept= 0)+
  geom_hline(yintercept = 0)+
  xlab(paste0("Dimensión 1: ",scales::percent(inPr[1]/sum(inPr),0.1))) +
  ylab(paste0("Dimensión 2: ", scales::percent(inPr[2]/sum(inPr),0.1)))

Como se construye directamente desde \(\texttt{ca}\), ya no presenta especularidad respecto a plot(ca(m)).

Análisis de Correspondencias Múltiples

Matriz disyuntiva o matriz de indicadores

Es una matriz donde cada línea coresponde a un individuo. Las filas se rellenan de 0 y 1. Se puede utilizar para descomposición en valores singulares.

farms <- MASS::farms


# Para crear una matriz

c <- model.matrix(~ -1+Mois+Manag+Use+Manure,data= farms,
                  contrasts.arg = list(Manure = contrasts(farms$Manure, contrasts = FALSE),
                                       Use = contrasts(farms$Use, contrasts = FALSE),
                                       Manag = contrasts(farms$Manag, contrasts = FALSE)
                                       )
                  )
head(c)

##   MoisM1 MoisM2 MoisM4 MoisM5 ManagBF ManagHF ManagNM ManagSF UseU1 UseU2 UseU3
## 1      1      0      0      0       0       0       0       1     0     1     0
## 2      1      0      0      0       1       0       0       0     0     1     0
## 3      0      1      0      0       0       0       0       1     0     1     0
## 4      0      1      0      0       0       0       0       1     0     1     0
## 5      1      0      0      0       0       1       0       0     1     0     0
## 6      1      0      0      0       0       1       0       0     0     1     0
##   ManureC0 ManureC1 ManureC2 ManureC3 ManureC4
## 1        0        0        0        0        1
## 2        0        0        1        0        0
## 3        0        0        0        0        1
## 4        0        0        0        0        1
## 5        0        0        1        0        0
## 6        0        0        1        0        0

Matriz de Burt

La matriz de Burt se obtiene de la matriz disyuntiva: es el producto de su traspuesta por la misma:

burt <- t(c) %*% c
head(burt)

##         MoisM1 MoisM2 MoisM4 MoisM5 ManagBF ManagHF ManagNM ManagSF UseU1 UseU2
## MoisM1       7      0      0      0       2       3       1       1     2     3
## MoisM2       0      4      0      0       1       0       1       2     2     2
## MoisM4       0      0      2      0       0       1       0       1     1     1
## MoisM5       0      0      0      7       0       1       4       2     2     2
## ManagBF      2      1      0      0       3       0       0       0     1     1
## ManagHF      3      0      1      1       0       5       0       0     2     1
##         UseU3 ManureC0 ManureC1 ManureC2 ManureC3 ManureC4
## MoisM1      2        1        1        3        1        1
## MoisM2      0        1        1        0        0        2
## MoisM4      0        0        1        1        0        0
## MoisM5      3        4        0        0        3        0
## ManagBF     1        0        2        1        0        0
## ManagHF     2        0        1        2        2        0

svd(c)$d^2

##  [1] 2.302404e+01 1.512265e+01 1.082962e+01 9.992020e+00 6.207453e+00
##  [6] 4.688870e+00 3.545108e+00 2.067642e+00 1.896087e+00 1.503231e+00
## [11] 9.017280e-01 2.215565e-01 1.290742e-31 7.101473e-32 6.076003e-32
## [16] 1.031706e-32

eigen(burt,
      symmetric = TRUE,
      only.values = TRUE
      )$values

##  [1]  2.302404e+01  1.512265e+01  1.082962e+01  9.992020e+00  6.207453e+00
##  [6]  4.688870e+00  3.545108e+00  2.067642e+00  1.896087e+00  1.503231e+00
## [11]  9.017280e-01  2.215565e-01  2.468499e-15  9.363299e-16 -2.524068e-16
## [16] -7.134271e-16

Los últimos autovalores son son diferentes por redondeamiento comutacional.

\(\text{MCA}\) a través de la \(\text{CA}\) de la Matriz de Burt

caburt <- ca(burt)

Representación de la \(\text{CA}\) de la Matriz de Burt para clases: método estándar

plot.ca(caburt)

Representación de la \(\text{CA}\) por Clase de la Matriz de Burt: con \(\texttt{ggplot2}\)

aburt <- caburt

# Inertia principal = eigenval(C %*% C^t) / n

inPr <- aburt$sv^2

# Construcción directamente desde de la library ca

r <- data.frame(Class= aburt$rownames,
                mass=aburt$rowmass,
                distChi=aburt$rowdist,
                inertia=aburt$rowinertia,
                aburt$rowcoord[,1:2]) %>% 
  mutate(Variable= gsub("([A-Z][a-z].*?)[A-Z].*", "\\1", Class)) %>% 
  mutate(Dim1 = Dim1 * sqrt(inPr[1]), # Se realiza el mismo esclado que en la CA bidomensional
         Dim2 = Dim2 * sqrt(inPr[2]))

# Gráfico idéntico

ggplot(r, aes(x=Dim1, y=Dim2, 
              color = Variable
              )) +
  geom_point()+
  geom_text_repel(aes(label = Class)) +
  geom_vline(xintercept= 0)+
  geom_hline(yintercept = 0)+
  xlab(paste0("Dimensión 1: ",scales::percent(inPr[1]/sum(inPr),0.1))) +
  ylab(paste0("Dimensión 2: ", scales::percent(inPr[2]/sum(inPr),0.1)))

En este caso, no se visualiza según inercia, ya que parece que la matriz de Burt tiene tendencia en sobrestimarla.

Representación de los individuos

Si realizamos la trasformación de Burt al revés,

tburt <- c %*% t(c)

Obtenemos una matriz de Burt referente a los individuos que, tratada como la anterior, puede generar gráficos como los siguientes.

tcaburt <- ca(tburt)

plot.ca(tcaburt)

# Inertia principal = eigenval(C %*% C^t) / n

inPr2 <- tcaburt$sv^2

# Construcción directamente desde de la library ca

i <- data.frame(Ind= tcaburt$rownames,
                mass=tcaburt$rowmass,
                distChi=tcaburt$rowdist,
                inertia=tcaburt$rowinertia,
                tcaburt$rowcoord[,1:2]) %>% 
  #mutate(Variable= gsub("([A-Z][a-z].*?)[A-Z].*", "\\1", Class)) %>% 
  mutate(Dim1 = Dim1 * sqrt(inPr2[1]), # Se realiza el mismo esclado que en la CA bidomensional
         Dim2 = Dim2 * sqrt(inPr2[2]))

# Gráfico idéntico

ggplot(i, aes(x=Dim1, y=Dim2
              #color = Variable,
              #size = inertia,
              )) +
  geom_point()+
  geom_text_repel(aes(label = Ind)) +
  geom_vline(xintercept= 0)+
  geom_hline(yintercept = 0)+
  xlab(paste0("Dimensión 1: ",scales::percent(inPr2[1]/sum(inPr2),0.1))) +
  ylab(paste0("Dimensión 2: ", scales::percent(inPr2[2]/sum(inPr2),0.1)))

\(\text{MCA}\) Realizado con \(\texttt{MCA}\) de \(\texttt{FactoMineR}\)

library(FactoMineR)

Como vemos, la visualización obtenida arriba, aplicando \(\texttt{ca}\) a la matriz de Burt és muy similar a la que produce utilizando \(\texttt{MCA}\).

Representación de las Clases según primera y segunda dimensión

plot(mca,invisible= "ind")

Representación de los Individuos según primera y segunda dimensión

plot(mca,invisible= "var")

Biplot \(\text{MCA}\) Individuos y Clases

plot(mca)