Técnicas de ordenación: PCA, MDS
Curso 013: Introducción al análisis de datos multivariantes
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/00Rteam/pub/clas/data/biom2003.dat")
head(biom)
Grupo Peso Altura Pie Hombros Brazos Caderas Sexo Ojos
1 1 60 163 37 41 68 95 1 1
2 1 52 166 37 37 70 87 1 2
3 1 61 172 39 39 69 91 1 2
4 1 73 181 43 50 78 101 2 2
5 1 53 172 39 39 72 89 1 1
6 1 63 169 40 37 66 96 1 2
Tipo
1 2
2 1
3 1
4 1
5 1
6 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
Min. :40.00 Min. :150.0 Min. :35.00
1st Qu.:56.00 1st Qu.:163.0 1st Qu.:38.00
Median :61.50 Median :170.0 Median :39.00
Mean :63.16 Mean :169.7 Mean :39.61
3rd Qu.:70.00 3rd Qu.:175.8 3rd Qu.:42.00
Max. :95.00 Max. :188.0 Max. :47.00
Hombros Brazos Caderas
Min. :32.00 Min. :63.00 Min. : 73.00
1st Qu.:38.00 1st Qu.:69.00 1st Qu.: 90.00
Median :40.00 Median :71.00 Median : 95.00
Mean :40.86 Mean :71.55 Mean : 95.21
3rd Qu.:44.00 3rd Qu.:74.75 3rd Qu.:100.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
Peso 1.0000000 0.7071186 0.7647123 0.6461115 0.5669994
Altura 0.7071186 1.0000000 0.8539068 0.5457743 0.7868916
Pie 0.7647123 0.8539068 1.0000000 0.6260485 0.6852809
Hombros 0.6461115 0.5457743 0.6260485 1.0000000 0.5446283
Brazos 0.5669994 0.7868916 0.6852809 0.5446283 1.0000000
Caderas 0.6891848 0.4202303 0.4289266 0.3785138 0.2786997
Caderas
Peso 0.6891848
Altura 0.4202303
Pie 0.4289266
Hombros 0.3785138
Brazos 0.2786997
Caderas 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
comp 1 4.0267249 67.112082
comp 2 0.8656409 14.427349
comp 3 0.5177578 8.629296
comp 4 0.3117111 5.195184
comp 5 0.1664386 2.773977
comp 6 0.1117267 1.862112
cumulative percentage of variance
comp 1 67.11208
comp 2 81.53943
comp 3 90.16873
comp 4 95.36391
comp 5 98.13789
comp 6 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
1 -1.2376961 0.6087399 0.5472820 0.14791608 -0.04442171
2 -1.9038719 -0.7226803 -0.2127355 0.15627710 0.04957810
3 -0.6381043 -0.2495886 -0.1913736 -0.49914831 0.06857999
4 3.1519143 -0.3780732 0.7518615 0.41096376 0.30339598
5 -0.7921862 -1.0059793 -0.2826909 0.07894064 0.44568833
6 -0.7792484 0.6756044 -0.5011026 -1.02488689 0.15276054
Dim.6
1 -0.232230101
2 -0.243800640
3 -0.470101232
4 -0.004718868
5 -0.167220217
6 -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
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
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
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
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
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 )
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.
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 )
- Hacer componentes principales
dfPCA <- PCA( df, scale.unit = TRUE )
- Visualizar resultados
fviz_pca_ind( dfPCA ); fviz_pca_var( dfPCA )
- Interpretar resultados
Escalado métrico multidimensional (MDS)
El escalado multidimensional (MDS) o coordenadas principales es una aplicación del PCA cuando la información que disponemos de 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) cuando la matriz es propiamente una matriz de distancias, y el MDS no métrico cuando no son distancias, sino similitudes o disimilitudes.
Ejemplo MDS métrico
Eurodist contiene las distancias en Km entre diferentes ciudades europeas,
#Cargar datos
data( eurodist )
dst <- eurodist #Directamente objeto dist no matrix ni data.frame
str( dst )
Class 'dist' atomic [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
( 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
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
[5] 7.893472e+05 5.816552e+05 2.623192e+05 1.925976e+05
[9] 1.450845e+05 1.079673e+05 5.139484e+04 -3.259629e-09
[13] -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
[21] -2.251844e+06
#Variabilidad acumulada
cumsum( abs( dstMDS$eig ) ) / sum( abs( dstMDS$eig ) )
[1] 0.4690928 0.7537543 0.7904600 0.8173197 0.8362709
[6] 0.8502358 0.8565337 0.8611578 0.8646411 0.8672332
[11] 0.8684672 0.8684672 0.8686952 0.8699690 0.8731434
[16] 0.8793217 0.8873088 0.8997033 0.9217710 0.9459359
[21] 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
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
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0
Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0
gear carb
Mazda RX4 4 4
Mazda RX4 Wag 4 4
Datsun 710 4 1
Hornet 4 Drive 3 1
Hornet Sportabout 3 2
Valiant 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
Mazda RX4 0.0000000 0.6153251 54.90861
Mazda RX4 Wag 0.6153251 0.0000000 54.89152
Datsun 710 54.9086059 54.8915169 0.00000
Hornet 4 Drive 98.1125212 98.0958939 150.99352
Hornet Sportabout 210.3374396 210.3358546 265.08316
Hornet 4 Drive Hornet Sportabout
Mazda RX4 98.11252 210.3374
Mazda RX4 Wag 98.09589 210.3359
Datsun 710 150.99352 265.0832
Hornet 4 Drive 0.00000 121.0298
Hornet Sportabout 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
[6] 0.9999839 0.9999886 0.9999927 0.9999958 0.9999980
[11] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[16] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[21] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[26] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[31] 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 )