Análisis multivariante. Ordenación I

Para ponernos en situación

Imaginemos que disponemos de los datos de altura de un conjunto de individuos. Si queremos representar esto para ver cómo se distribuyen, con un sólo eje sería suficiente:

Si además de la altura queremos representar el peso de estos mismos individuos, para ver la relación de ambas variables —por ejemplo ver si el peso aumenta con la altura—, necesitaríamos dos ejes para representar toda la información:

Os imagináis por donde vamos, ¿verdad?. Si queremos ver esas relaciones con 3, 4, … 100, … 500, variables, necesitaríamos otros tantos ejes y esto ya se escapa de nuestra capacidad de representación y casi que entendimiento.

Si encontramos que existe una cierta correlación entre dos variables, podemos encontrar una función que las relacione, de manera que podríamos construir una nueva variable que sea combinanción lineal de las anteriores y que aglutine prácticamente toda la información —puede haber una pequeña pérdida de la misma—. En este caso, los datos de la nueva variable quedarían representados así:

Como podemos, ver la nueva variable —que contiene la información combinada de peso y altura— por si sóla explica bastante bien la distribución de los individuos, solamente hay una mínima cantidad de información que no queda recogida por esta variable, pero que la recoge otra variable denominada en este gráfico nueva variable 2. Para conseguirlo, aparentemente lo único que hemos tenido que hacer es rotar los datos 45^o —realmente lo que hemos rotado han sido lo ejes—.

Las correlaciones entre las variables, por desgracia —o por suerte— no son perfectas, no tienen correlación 1 o -1, sino que el azar, aleatoriedad, diversidad, variabilidad, o como lo queramos llamar, hace que nos alejemos de esa correlación perfecta y por tanto necesitemos otra variable que recoja esas “migajas” de información que le ha dejado la primera nueva variable. Ahora la nueva variable tiene una información —mucha— y la nueva variable 2 tiene otra —muy poca—, pero lo importante es que estas dos variables no comparten nada de información entre ellas, son icorreladas.

En resumen, lo que queremos decir con esto, es que las técnicas multivariantes, son en general un conjunto de técnicas que intentan reducir las dimensiones de los datos originales, creando nuevas variables que nos permitan “ver el retrato” de nuestros datos con el mínimo número de variables posible y de esta manera podamos interpretar las relaciones existente entre los individuos basadas en todas estas variables.

Análisis de componentes principales

Es una técnica multivariante de reducción de dimensiones, cuya misión consiste en extraer la información importante de las variables originales, transformańdolas en nuevas variables llamadas componentes o factores, cumpliendo éstos que sean icorrelados entre si (ortogonales). Es un problema de álgrebra lineal, cada componente es resultado de la combinación lineal de las variables originales. De forma gráfica consistiría en una rotación de los datos sobre los ejes, para conseguir que la mayor parte de los datos caigan sobre ellos.

Cuanta mayor correlación haya entre las variables originales, mayor será el número de individuos explicados por la primera componente, y por tanto menos información quedará para ser explicada por la siguiente.

La interpretación de los factores es la clave del análisis, y es algo que no viene dado por el propio análisis, sino que hay que deducirlo después de observar la relación entre las variables y los factores.

Definiciones no formales

Matriz de correlaciones:: el PCA sólo tiene sentido si hay alta correlación entre variables, por tanto la matriz de correlaciones entre las variables puede ser un buen punto de partida.
Factores:: o componentes principales, son cada una de las nuevas variables, combinación lineal de las originales.
Matriz factorial:: matriz que contiene los coeficientes factoriales, es decir, las correlaciones entre componentes (columnas) y las variables originales (filas).
Autovalores:: o valores propios o eigenvalues, representan las varianzas de los factores. Ordenadas de mayor a menor.
Autovectores:: o vectores propios o eigenvectors, se corresponden con las columnas de la matriz factorial. Cada columna un autovector.

Instalación de las librerías

Para comenzar instalaremos los paquetes FactoMineR y factoextra y otros paquetes que iremos necesitando a lo largo del análisis

install.packages( "FactoMineR" )
install.packages( "factoextra" )
install.packages( "corrplot" )
install.packages( "psych" )
install.packages( "Factoshiny" )
install.packages( "rgl" )

y cargamos las las librerías

Cargar y preparar los datos

Vamos a trabajar con un conjunto de datos, correspondientes a las medidas biométricas de alumnos de la Universidad de Murcia. La tabla de datos consta de 98 observaciones de las que se han obtenido el peso, altura, talla pie, longitud de brazos, circunferencia de cadera, color de ojos, sexo y un valor numérico que representa un morfotipo autoasignado por los propios alumnos.

biom <- read.table( "http://ares.inf.um.es/mmcl/recursos/biom2003.dat" )
head( biom )

  Grupo Peso Altura Pie Hombros Brazos Caderas Sexo Ojos Tipo
1     1   60    163  37      41     68      95    1    1    2
2     1   52    166  37      37     70      87    1    2    1
3     1   61    172  39      39     69      91    1    2    1
4     1   73    181  43      50     78     101    2    2    1
5     1   53    172  39      39     72      89    1    1    1
6     1   63    169  40      37     66      96    1    2    1

Las variables grupo, sexo, ojos y tipo, a pesar de ser numéricas, no están representando una medida sino una categoría. Estas variables son conocidas en R como factores que no hay que confundir con las componentes pricipales que en ocasiones también se denominan factores.

Estadísticos descriptivos

De forma preliminar, podemos obtener una serie de estadísticos básicos que nos sirven para conocer un poco nuestros datos. Esto se podría hacer con la función summary

summary(biom[,2:7])

      Peso           Altura           Pie           Hombros          Brazos         Caderas      
 Min.   :40.00   Min.   :150.0   Min.   :35.00   Min.   :32.00   Min.   :63.00   Min.   : 73.00  
 1st Qu.:56.00   1st Qu.:163.0   1st Qu.:38.00   1st Qu.:38.00   1st Qu.:69.00   1st Qu.: 90.00  
 Median :61.50   Median :170.0   Median :39.00   Median :40.00   Median :71.00   Median : 95.00  
 Mean   :63.16   Mean   :169.7   Mean   :39.61   Mean   :40.86   Mean   :71.55   Mean   : 95.21  
 3rd Qu.:70.00   3rd Qu.:175.8   3rd Qu.:42.00   3rd Qu.:44.00   3rd Qu.:74.75   3rd Qu.:100.00  
 Max.   :95.00   Max.   :188.0   Max.   :47.00   Max.   :53.00   Max.   :86.00   Max.   :113.00

Aunque si queremos personalizar los estadísticos y sacarlos en forma de una tabla más manejable puede hacerse de la siguiente manera

bm <- biom[,2:7]
biom_summary <- data.frame( Min = apply( bm, 2, min ), 
                            Q1 = apply( bm, 2, quantile, 1/4 ),
                            Mediana = apply( bm, 2, median ),
                            Media = apply( bm, 2, mean ), 
                            Q3 = apply( bm, 2, quantile, 3/4 ), 
                            Max = apply( bm, 2, max ), 
                            Var = apply( bm, 2, var ) )
( biom_summary )

        Min  Q1 Mediana     Media     Q3 Max        Var
Peso     40  56    61.5  63.16327  70.00  95 104.385441
Altura  150 163   170.0 169.71429 175.75 188  76.206186
Pie      35  38    39.0  39.61224  42.00  47   7.621292
Hombros  32  38    40.0  40.85714  44.00  53  21.402062
Brazos   63  69    71.0  71.55102  74.75  86  18.765411
Caderas  73  90    95.0  95.21429 100.00 113  59.118557

Vimos en la introducción que lo adecuado era trabajar con los datos tipificados, y que la matriz de correlaciones era la matriz de covarianzas de datos tipificados. Vamos a calcular esa matriz de correlaciones y visualizar las relaciones entre las variables.

( biom_cor <- cor( bm ) )

             Peso    Altura       Pie   Hombros    Brazos   Caderas
Peso    1.0000000 0.7071186 0.7647123 0.6461115 0.5669994 0.6891848
Altura  0.7071186 1.0000000 0.8539068 0.5457743 0.7868916 0.4202303
Pie     0.7647123 0.8539068 1.0000000 0.6260485 0.6852809 0.4289266
Hombros 0.6461115 0.5457743 0.6260485 1.0000000 0.5446283 0.3785138
Brazos  0.5669994 0.7868916 0.6852809 0.5446283 1.0000000 0.2786997
Caderas 0.6891848 0.4202303 0.4289266 0.3785138 0.2786997 1.0000000

Podemos visualizar estas correlaciones en forma de gráficos de dispersión

pairs( bm )

y también podemos visualizar esas dependencias entre las varibles mediante un correlograma

corrplot( biom_cor, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

Esto nos permite ver que existen variables altamente correlacionadas como Altura y Peso o Altura y Pie y otras muy poco correlacionas como Brazos y Caderas o Pie y Caderas.

Numéricamente, hay distintas técnicas para determinar si existe una alta correlaciones en datos multivariantes de forma conjunta que justifiquen el hacer un PCA. Uno de ellos es el índice de Kaiser-Meyer-Olkin (KMO), que compara los valores de correlaciones entre pares de variables. El punto de partida de ese índice es una matriz de correlaciones entre las variables observadas. Si el índice es próximo a 1 el PCA se puede hacer, mientras que si es próximo a 0, el PCA no será relevante.

( biom_kmo <- psych::KMO( bm ) )

Kaiser-Meyer-Olkin factor adequacy
Call: psych::KMO(r = bm)
Overall MSA =  0.8
MSA for each item = 
   Peso  Altura     Pie Hombros  Brazos Caderas 
   0.80    0.78    0.82    0.88    0.83    0.72

También podemos ver las relaciones entre individuos y variables mediante gráficos de estrellas. Este gráfico nos da idea que cuánto pesa cada variable sobre cada individuo.

stars( bm , labels = abbreviate( rownames( bm ), 5 ), 
      nrow = 7, key.loc = c( 18, 19 ), full = TRUE )

También podemos colorear a los individuos en base a alguna categoría como por ejemplo según su sexo

color <- as.character( factor( biom$Sexo, labels = c( "pink", "skyblue" ) ) )

stars( bm, labels = abbreviate( rownames( bm ), 5 ), 
      nrow = 7, key.loc = c( 18, 19 ), full = TRUE, col.stars = color )

Aplicar componentes principales

Ya hemos explorado nuestros datos, los conocemos bastante bien y hemos visto que existen altas correlaciones entre las variables que justificaría poder hacer un PCA.

biom_pca <- PCA( bm, scale.unit = TRUE, ncp = 6, graph = FALSE )
print( biom_pca )

**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 98 individuals, described by 6 variables
*The results are available in the following objects:

   name               description                          
1  "$eig"             "eigenvalues"                        
2  "$var"             "results for the variables"          
3  "$var$coord"       "coord. for the variables"           
4  "$var$cor"         "correlations variables - dimensions"
5  "$var$cos2"        "cos2 for the variables"             
6  "$var$contrib"     "contributions of the variables"     
7  "$ind"             "results for the individuals"        
8  "$ind$coord"       "coord. for the individuals"         
9  "$ind$cos2"        "cos2 for the individuals"           
10 "$ind$contrib"     "contributions of the individuals"   
11 "$call"            "summary statistics"                 
12 "$call$centre"     "mean of the variables"              
13 "$call$ecart.type" "standard error of the variables"    
14 "$call$row.w"      "weights for the individuals"        
15 "$call$col.w"      "weights for the variables"

La función nos devuelve un objeto que contiene toda la información del análisis.

Autovalores

Si quisiésemos acceder a los autovalores para ver la cantidad de varianza acumulada por cada componente, tendríamos que hacer:

( biom_pca$eig )

       eigenvalue percentage of variance cumulative percentage of variance
comp 1  4.0267249              67.112082                          67.11208
comp 2  0.8656409              14.427349                          81.53943
comp 3  0.5177578               8.629296                          90.16873
comp 4  0.3117111               5.195184                          95.36391
comp 5  0.1664386               2.773977                          98.13789
comp 6  0.1117267               1.862112                         100.00000

Lo que nos indica que la primera componente acumula 4.0267249 de la varianza, lo que supone un 67.1120816% del total. La suma de todos los autovalores será igual al número de componentes:

( sum( biom_pca$eig[ ,1  ] ) )

[1] 6

Vemos que la mayor parte de la varianza queda explicada por la primera componente y que entre las dos primeras ya se explica un 81.5394307% de la varianza. Con lo que hemos pasado de tener 6 variables para explicar nuestros datos, a tener dos componentes que contienen prácticamente toda la información.

Un gráfico de sedimentación de los autovalores nos permite visualizar mejor la varianza acumulada por cada componente.

fviz_screeplot( biom_pca, ncp = 6 )

Variables

Otros parámetros de interés que nos devuelve el análisis sería:

Coord = cor: son las correlaciones entre las variables originales y los factores (nuevas variables). Son también las coordenadas de los vectores que representas las variables.

¿Por qué las correlaciones y las coordenadas son iguales? Esto ocurre porque le hemos pedido a la función que nos tipifique los datos ( scale.unit = TRUE ). Si los datos no se tipificasen, las coordenadas no estarían en el rango -1 a 1 como lo están siempre las correlaciones.

Las columnas de la matriz de correlaciones son los autovectores, es decir, los factores que se aplican sobre los datos para realizar la rotación y obtener las nuevas variables.
cos2: es igual al cuadrado de cor $cos2 = cor^2$. Mide la calidad de la representación de las variables en el plano factorial. La suma de los cos2 para una variable es 1.
```
( apply(biom_pca$var$cos2, 1, sum ) )
```
```
   Peso  Altura     Pie Hombros  Brazos Caderas 
      1       1       1       1       1       1 
```

contrib: Mide el porcentaje de peso de cada variable sobre cada factor. La suma por columnas es 100.

( apply(biom_pca$var$contrib, 2, sum ) )

Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6 
  100   100   100   100   100   100

Individuos

Los parámetros indicados anteriormente en los resultados también se obtienen en el caso de los individuos, lo que permite hacer una interpretación desde esa perspectiva. Es decir, podríamos ver qué individuo o individuos son los que más contribuyen o influyen sobre las componentes, ver las nuevas coordenadas de nuestras observaciones sobre el plano factorial.

Coord: coordenadas de las observaciones sobre el plano factorial. Son los valores que toman las nuevas variables sobre los individuos. Sería la nueva tabla de datos

head( biom_pca$ind$coord )

       Dim.1      Dim.2      Dim.3       Dim.4       Dim.5        Dim.6
1 -1.2376961  0.6087399  0.5472820  0.14791608 -0.04442171 -0.232230101
2 -1.9038719 -0.7226803 -0.2127355  0.15627710  0.04957810 -0.243800640
3 -0.6381043 -0.2495886 -0.1913736 -0.49914831  0.06857999 -0.470101232
4  3.1519143 -0.3780732  0.7518615  0.41096376  0.30339598 -0.004718868
5 -0.7921862 -1.0059793 -0.2826909  0.07894064  0.44568833 -0.167220217
6 -0.7792484  0.6756044 -0.5011026 -1.02488689  0.15276054 -0.101472595

Cos2: Proyección sobre los ejes. Mide la contribución de cada factor sobre un individuo concreto. La suma por filas será 1.

apply( biom_pca$ind$cos2, 1, sum )

 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1

Contrib: Mide el porcentaje de peso de cada individuo sobre un factor. La suma por columnas es 100.

apply( biom_pca$ind$contrib, 2, sum )

Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6 
  100   100   100   100   100   100

Dist: La distancia absoluta de cada individuo al origen de coordenadas del plano factorial. Se calcula, para un individuo concreto, con la fórmula de la distancia euclídea, como la raíz cuadrada de la suma de las coordenadas de ese individuo al cuadrado.
```
head( biom_pca$ind$dist, 5 )
```
```
        1         2         3         4         5 
1.5098863 2.0684706 0.9904256 3.3020844 1.3972472 
```
Si calculamos a mano la distancia de los individuos al origen de coordenadas (para los 5 primeros):
```
apply( biom_pca$ind$coord[ 1:5, ], 1, function( x ) sqrt( sum( x^2 ) ) )
```
```
        1         2         3         4         5 
1.5098863 2.0684706 0.9904256 3.3020844 1.3972472 
```

Representación gráfica de las variables y los individuos

Aunque las correlaciones entre variables y factores nos dan información sobre el peso de cada variable sobre cada factor, la representación gráfica de éstas sobre el plano factorial es fundamental para entender esas relaciones y poder interpretar mejor los resultados.

Variables con factoextra

fviz_pca_var( biom_pca, axes = c( 1, 2), col.var = "cos2", alpha.var = "contrib" ) + theme_grey()

Vimos anteriormente que la suma de cos2 de una variable determinada sobre cada factor es 1. Esto significa que cada vector debería estar tocando el perímetro de la circunferencia unidad, pero no lo está haciendo ninguna, ¿por qué?. Si observamos por ejemplo la variable Caderas, vemos que esta muy cerca de tocar dicho perímetro, su proyección sobre las dimensiones 1 y 2 (componentes) indica su contribución a éstas, pero aún le falta algo de contribución que debe estar repartida por otra u otras dimensiones. Si está variable solo tuviese peso sobre las dos primeras dimensiones estaría tocando la circunferencia.

Individuos con factoextra

fviz_pca_ind( biom_pca, col.ind = color, alpha.ind = "contrib" ) + theme_grey()

fviz_pca_ind( biom_pca, habillage = as.factor( biom$Sexo ), addEllipses = TRUE, ellipse.level = 0.99 )

Podemos colorear las observaciones según alguna variable, en este caso el sexo. Además podemos hacer que las variables que más contribuyen en este plano factorial, se resalten más que las que menos influencia tienen. También tenemos la posibilidad de dibujar elipses alrededor de cada grupo con un cierto nivel de confianza.

Todos juntos

Finalmente, igual que en el caso de la representación gráfica manual, podemos hace un gráfico conocido como biplot donde se representan variables e individuos simultáneamente.

fviz_pca_biplot( biom_pca, col.ind = color )

Si somos un poco observadores, veremos que este gráfico esta un poco “trucado”. Ya hemos comentado que la suma total de cos2 para cada variable es 1. Por tanto ningún vector de ninguna variable puede superar ese valor, sin embargo en esta gráfica, si nos fijamos en la escala vemos que esto si que sucede. Simplemente es una pequeña licencia para representar de forma conjunta observaciones y variable, dando más importancia a las coordenadas de las observaciones que a las de las variables.

Resultados interactivos

Factoshiny::PCAshiny(bm)

Esta función arranca una aplicación shiny en el navegador, donde se puede interactuar con los resultados y los gráficos.

Variables suplementarias

Se pueden añadir variables suplementarias al análisis que no se tienen en cuenta en éste, pero que pueden ayudar a interpretar los resultados. Estas variables también se conocen variables ilustrativas.

biom$Sexo <- as.factor( biom$Sexo )
biom_pca_sup <- PCA( biom[ ,2:8 ], quali.sup = 7  )

plot( biom_pca_sup, choix = "ind" )

Gráfico a mano

Una vez que hemos visto cómo representar los resultados fácilmente con las distintas funciones del paquete factoextra, vamos a construir este mismo gráfico a mano, para entender cómo se construye. Saber cómo están construidos los gráficos siempre facilita tanto la comprensión de la técnica, como la interpretación de los resultados.

Con este gráfico podríamos también hacer una predicción sobre dónde caería un nuevo individuo y nos serviría para comprobar si nuestra interpretación está siendo correcta.

par( new = TRUE )

Warning in par(new = TRUE): llamada par(new=TRUE) sin gráfico

#par( cex = 0.8 )
#mapa de variables
x <- biom_pca$var$coord[, 1]
y <- biom_pca$var$coord[, 2]
plot( x, y, type="n", xlim = c( -1.0, 1.0 ), 
     ylim = c( -1.0, 1.0 ), pch = 20, cex = 2.5, col = c( 1:6 ),
     xlab = "PC1", ylab = "PC2", asp = 1 )
abline( v = 0, h = 0 )
title( "Representación manual de componentes principales" )

#añadir direcciones
arrows( 0, 0, x, y, col= rgb( .5, .5, .5 ) , length = 0.15, angle = 10, lwd = 5 )
arrows( 0, 0, x, y, col = c( 1:6 ), length = 0.15, angle = 10, lwd = 2 )
text( x, y, labels=rownames(biom_pca$var$coord), cex = 0.8, col = c( 1:6 ), adj = c( -.3, 0 ) )

## Círculo unidad
symbols( 0, 0, circles = 1, add = TRUE, xlim = c( -0.8, 0.8 ), inches = FALSE )

#añadir observaciones
xc <- biom_pca$ind$coord[, 1]
yc <- biom_pca$ind$coord[, 2]
par( new = TRUE )
simbol <- ifelse( biom$Sexo == 1, 18, 20 )
plot( xc, yc, col = color, type = "p", pch = simbol, xlab = "",
      ylab = "", xaxt = "n", yaxt = "n", cex = 1.3 )

#predicción
ni <- as.data.frame( matrix( c( 77, 184, 44, 54, 81, 90 ), nrow = 1 ) )
colnames( ni ) <- colnames( bm )
nuevodato <- predict.PCA( biom_pca, newdata = ni, scale = TRUE )
ndNX<-  (nuevodato$coord[1] - mean( biom_pca$ind$coord[,1] ) ) / sd(biom_pca$ind$coord[, 1] )
ndNY<-  (nuevodato$coord[2] - mean( biom_pca$ind$coord[,2] ) ) / sd(biom_pca$ind$coord[, 2] )
points( ndNX, ndNY, pch = 20, cex = 2.5, col = 2 )

Con este gráfico podemos observar varias cosas. Por un lado ahora las variables están en su escala, ya que ninguna pasa de 1 y los individuos en la suya. Por otro lado podemos ver lo útil que puede ser introducir una observación conocidad en el gráfico (punto rojo) para comprobar si estamos realizando una correcta interpretación. Finalmente, hacer el gráfico “a mano”, aunque nos ha servido para ver cómo se generan, también no ha enseñado que utilizar las funciones existentes para ello es mucho mejor y menos trabajoso.

Interpretar los resultados

El eje de abscisas es una medida biométrica global, los individuos en general más grandes se sitúan en la parte positiva del eje.

En el eje de ordenadas, caderas y peso contribuyen positivamente a la segunda componente y el resto de variables negativamente, sobre todo longitud de brazos y altura, podríamos interpretar este eje como una medida de la anchura de los individuos relacionada con el peso. Los individuos situados en la parte superior serán de mayor peso y con caderas más anchas.

Se podría decir que el primer eje es una medida de tamaño y el segundo una medida de forma. Los individuos a la izquierda del eje son más pequeños y a la derecha más grandes. Mientras que los de la parte superior son más anchos en general y en la parte inferior más estilizados o “espigados”.

Si cada variable contribuyese en igual proporción, puesto que tenemos 6 variables, cada una contribuiría $1/6$, o lo que es lo mismo un $16.67\%$. Podemos representar la contribución de las variables a los factores mediante gráficos de barras.

miscolores <- colors()[c( 261, 300, 201, 230, 350 )]
miscolores <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2")
barplot( t( biom_pca$var$contrib ), beside = TRUE,
         border = 1, col = miscolores, 
         legend.text = colnames(biom_pca$var$contrib),
         args.legend = list(x = "top", ncol = 5), 
         ylim = c( 0, 90 ) )
abline( h = 16.67, col = "red", lwd = 1.5 )

Visualización de 3 componentes

Los gráficos en 3 dimensiones nos permite tener una visión un poco más clara y realista de la distribución de los individuos y las variables en este espacio multidimensional. Hay que tener en cuenta que con los datos de este ejemplo, teníamos alrededor del 90 % de la información representada en 3 ejes, lo que quiere decir que un gráfico 3D prácticamente no deja incógnitas sobre la distribución de observaciones y variables. Sin embargo, si la varianza explicada fuese bastante inferior, aún en un espacio tridimensional veríamos que los vectores no llegan a tocar la superficie de la esfera. Necesitaríamos un espacio n–dimensional para poder verlo.

Visualización 3D de los individuos

# individuos en 3D
plot3d( biom_pca$ind$coord[, 1:3 ], col = biom$Sexo )

You must enable Javascript to view this page properly.

Visualización 3D de las variables

# Variables en 3D
matrix3d <- matrix( c(  0.99924529, 0.01935360, -0.03367962, 0,
                       -0.01976939, 0.99973202, -0.01205689, 0,
                        0.03343736, 0.01271362,  0.99936002, 0,
                        0,          0,           0,          1), 
                    nrow = 4,  ncol = 4, byrow = TRUE)

coords <- NULL
for ( i in 1 : nrow( biom_pca$var$coord ) ) {
  coords <- rbind( coords, rbind( c( 0, 0 , 0 ), biom_pca$var$coord[ i, 1:3 ] ) )
}
plot3d( biom_pca$var$coord[, 1:3], type = "n", xlim = c(-1,1), ylim = c(-1,1), zlim = c(-1,1))

spheres3d( 0, 0, 0, radius = 1, alpha = c( "0.2" ), color = c( "gray" )  )

lines3d(coords, col = "red" , lwd = 4, add = TRUE )
text3d( biom_pca$var$coord[ , 1:3 ], texts = rownames( biom_pca$var$coord ), col = "black" )

par3d( userMatrix = matrix3d )

You must enable Javascript to view this page properly.

Resumiendo

A pesar de que parezca demasiado complejo y demasiada información, realmente es muy sencillo aplicar componentes principales en R. En apenas cuatro o cinco líneas de código estaría resuelto todo el análisis con gráficos incluidos.

Como breve guía, básicamente el análisis sería:

Importar los datos a R df <- read.table( mydata ).
Comprobar que existe correlación entre las variables pairs( df ), cor( df ).
Calcular componentes principales dfPCA <- PCA( df, scale.unit = TRUE ).
Visualizar resultados fviz_pca_ind( dfPCA ); fviz_pca_var( dfPCA ).
Interpretar resultados para esto no hay función :-(

Escalado métrico multidimensional (MDS)

El escalado multidimensional (MDS) o coordenadas principales es una aplicación del PCA cuando la información que de la disponemos para los individuos es una matriz de distancias o similitudes, cuyo objetivo de representar en un espacio de reducidas dimensiones una configuración de puntos que represente a los individuos conservando al máximo la información original.

Podría entenderse como una técnica complementaria al PCA, pues mientras el PCA explora las relaciones entre las variables a partir de una matriz de covarianzas o correlaciones, el MDS explora la estructura o configuración de los individuos a partir de una matriz de distancias o proximidades.

En el escalado multidimensional se distigue el escalado métrico multidimensional (MMDS) en el que asumimos que las relaciones entre distancias y disimmilaridades son lineales y se resuelve mediante una función de mínimos cuadrados, y el escalado no métrico multidimensional (NMDS), en el que no asumimos que exista esa relación lineal, resolviéndose mediante una regresión monótona de mínimos cuadrados.

No os preocupéis, no vamos a entrar en tanto detalle, únicamente lo mencionamos para que conozcáis de su existencia. Aquí veremos el caso más simple de MMDS.

Ejemplo MDS métrico

eurodist contiene las distancias en kilómetros por carretera entre diferentes ciudades europeas,

#Cargar datos
data( eurodist )
dst <- eurodist #Directamente objeto dist no matrix ni data.frame
str( dst )

 'dist' num [1:210] 3313 2963 3175 3339 2762 ...
 - attr(*, "Size")= num 21
 - attr(*, "Labels")= chr [1:21] "Athens" "Barcelona" "Brussels" "Calais" ...

El objeto generado es de tipo dist, puede ser visualizado como una matriz de datos, convirtiéndolo como tal mediante la función as.matrix, por ejemplo:

( as.matrix( dst )[ 1:4, 1:4 ] )

          Athens Barcelona Brussels Calais
Athens         0      3313     2963   3175
Barcelona   3313         0     1318   1326
Brussels    2963      1318        0    204
Calais      3175      1326      204      0

Para realizar el análisis:

par( cex = 0.7 )

# Escalado métrico
dstMDS <- cmdscale( dst, k = 3, eig = TRUE )
summary( dstMDS )

       Length Class  Mode   
points 63     -none- numeric
eig    21     -none- numeric
x       0     -none- NULL   
ac      1     -none- numeric
GOF     2     -none- numeric

str( dstMDS )

List of 5
 $ points: num [1:21, 1:3] 2290.3 -825.4 59.2 -82.8 -352.5 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:21] "Athens" "Barcelona" "Brussels" "Calais" ...
  .. ..$ : NULL
 $ eig   : num [1:21] 19538377 11856555 1528844 1118742 789347 ...
 $ x     : NULL
 $ ac    : num 0
 $ GOF   : num [1:2] 0.79 0.91

dstMDS$eig

 [1]  1.953838e+07  1.185656e+07  1.528844e+06  1.118742e+06  7.893472e+05  5.816552e+05  2.623192e+05  1.925976e+05
 [9]  1.450845e+05  1.079673e+05  5.139484e+04 -2.793968e-09 -9.496124e+03 -5.305820e+04 -1.322166e+05 -2.573360e+05
[17] -3.326719e+05 -5.162523e+05 -9.191491e+05 -1.006504e+06 -2.251844e+06

#Variabilidad acumulada
cumsum( abs( dstMDS$eig ) ) / sum( abs( dstMDS$eig ) )

 [1] 0.4690928 0.7537543 0.7904600 0.8173197 0.8362709 0.8502358 0.8565337 0.8611578 0.8646411 0.8672332 0.8684672
[12] 0.8684672 0.8686952 0.8699690 0.8731434 0.8793217 0.8873088 0.8997033 0.9217710 0.9459359 1.0000000

#Plot mapa
x <- dstMDS$points[,1]
y <- dstMDS$points[,2]
plot( x, y, type = "n" )
text( x, y, labels = rownames( dstMDS$points ), cex= .75 )

#flip-flop  vertical (para que el norte este arriba)
plot( x, -1 * y, type = "n" )
text( x, -1 * y, labels = rownames( dstMDS$points ), cex= .75 )

#girando el mapa 45º más espejado horizontal y vertical
x1 <- ( x * cos( -pi / 4 ) - y * sin( -pi / 4 ) )
y1 <-  -1 *( x * sin( -pi / 4 ) + y * cos( - pi / 4 ) )

plot( x1, y1, type= "n" )
symbols( as.vector( x1 ), as.vector( y1 ), circles = rep( 50, 21 ), 
         inches = FALSE, add = TRUE, bg = "pink" )
text( x1, y1, labels = rownames( dstMDS$points ), cex= .65 )

#3d de eurodist
dstMDS <- cmdscale( dst, k = 3, eig = TRUE )
x <- dstMDS$points[,1]
y <- dstMDS$points[,2]
z <- dstMDS$points[,3]
plot3d( x, y, z, col = c(1:50), type = "n", size = 2, box = FALSE, axes = FALSE )
text3d( x, y, z, texts = rownames( dstMDS$points ), cex= .75 )

You must enable Javascript to view this page properly.

Ejemplo MDS para clase

En este ejemplo, la tabla de datos está formada por distintos modelos de coches (observaciones) de los que tenemos algunas características (variables). Queremos ver la relación existente entre los modelos, cómo de cerca o lejos están en base a sus características.

En este caso, no disponemos de un objeto distancia, sino una matriz de datos en base a la cual se calcularán las distancias.

data( "mtcars" )
head( mtcars )

                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

dist_df <- dist( mtcars )
carMDS <- cmdscale( dist_df, k = 3, eig = TRUE )
as.matrix( dist( mtcars[ 1:5, ] ) ) # distancias de los 5 primeros elementos

                    Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive Hornet Sportabout
Mazda RX4           0.0000000     0.6153251   54.90861       98.11252          210.3374
Mazda RX4 Wag       0.6153251     0.0000000   54.89152       98.09589          210.3359
Datsun 710         54.9086059    54.8915169    0.00000      150.99352          265.0832
Hornet 4 Drive     98.1125212    98.0958939  150.99352        0.00000          121.0298
Hornet Sportabout 210.3374396   210.3358546  265.08316      121.02976            0.0000

str( carMDS )

List of 5
 $ points: num [1:32, 1:3] -79.6 -79.6 -133.89 8.52 128.69 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
  .. ..$ : NULL
 $ eig   : num [1:32] 577879.5 45113.6 292.4 52.9 25.5 ...
 $ x     : NULL
 $ ac    : num 0
 $ GOF   : num [1:2] 1 1

#variabilidad acumulada
cumsum( abs( carMDS$eig ) ) / sum( abs( carMDS$eig ) )

 [1] 0.9269989 0.9993673 0.9998362 0.9999212 0.9999620 0.9999839 0.9999886 0.9999927 0.9999958 0.9999980 1.0000000
[12] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[23] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000

# representacion 2D
plot( carMDS$points[ , 1 ], carMDS$points[ , 2 ], type = "n", axes = FALSE, xlab = "", ylab = "" )
text( carMDS$points[ , 1 ], carMDS$points[ , 2 ], rownames( carMDS$points ), cex = 0.75 )

# En 3d
x <- carMDS$points[ , 1 ]
y <- carMDS$points[ , 2 ]
z <- carMDS$points[ , 3 ]
plot3d( x, y, z , type = "n", size = 1, box = FALSE, axes = FALSE, col = c( 1:21 ) )
text3d( x, y, z, texts  = rownames( carMDS$points ), cex=0.7 )

You must enable Javascript to view this page properly.

# Animación final. Ejecutar directamente en Rstudio
plot3d( x, y, z, type = "s", size = 0.5, box = FALSE, axes = FALSE, col = c( 1:21 ) )
text3d( x, y, z, texts  = rownames( carMDS$points ), cex = .5 )
play3d( spin3d( axis = c( 0, 1, 0 ), rpm = 4 ), duration = 10 )

Manos a la obra

Es hora de hacer un ejemplito sencillo para repasar y hacernos pensar un poco, asi que… manos a la obra.

Referencias

Este tutorial, tiene algo más de matemática que nuestro curso, aún así es bastante sencillo y con buenos ejemplos de aplicaciones (visitar)
Un artículo sobre la familia de componentes principales con R aquí
Por supuesto la página del archiconocido paquete FactoMineR
Un librito muy apañado sobre análisis multivariate general con R, en su capítulo 1 trata el análisis de componentes principales (visitar).

Análisis multivariante. Ordenación I

Módulo 4

Para ponernos en situación

Análisis de componentes principales

Definiciones no formales

Instalación de las librerías

Cargar y preparar los datos

Estadísticos descriptivos

Aplicar componentes principales

Autovalores

Variables

Individuos

Representación gráfica de las variables y los individuos

Variables con factoextra

Individuos con factoextra

Todos juntos

Resultados interactivos

Variables suplementarias

Gráfico a mano

Interpretar los resultados

Visualización de 3 componentes

Visualización 3D de los individuos

Visualización 3D de las variables

Resumiendo

Escalado métrico multidimensional (MDS)

Ejemplo MDS métrico

Ejemplo MDS para clase

Manos a la obra

Referencias