Escalado multidimensional

Introducción

El Análisis de Componentes Principales (ACP) se utiliza para obtener un “mapa” de los datos en dimensión reducida que preserve tanto como sea posible las distancias euclídeas entre las observaciones en su espacio original.

El objetivo del Escalado Multidimensional (Multidimensional Scaling, MDS) es similar ya que pretende obtener una representación gráfica en dimensión reducida pero a partir de una matriz de distancias entre los objetos o individuos. Esta matriz de distancias se puede calcular desde una matriz de datos \(\mathbf{X}\) sobre algunas variables medidas o puede ser una matriz de distancias o similaridades obtenida directamente como, por ejemplo, de un juicio de expertos.

Existen muchos métodos de escalado multidimensional, pero siguiendo el apartado 5.2 del libro de Brian S. Everitt (1) vamos a explicar únicamente el método clásico.

El método clásico

Como los otros métodos, el método clásico tiene el objetivo de representar una matriz de distancias (o de similaridades) en un mapa o modelo geométrico. Un modelo geométrico no es más que un conjunto de puntos \(\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_n\) en un espacio de \(q\) dimensiones de forma que cada punto represente uno de los objetos o individuos y que las distancias entre dichos puntos se ajusten lo máximo posible a las distancias entre los objetos. Así pues, el objetivo es doble, el primero es determinar la dimensión \(q\) del modelo y además obtener las coordenadas (de dimensión \(q\)) de los puntos \(\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_n\) que proporcionan un “buen” ajuste a las distancias observadas o propuestas. Ese “buen” ajuste se deberá medir con algún tipo de medida o índice. La idea es que si dos objetos tienen una distancia grande entre ellos, los puntos que los representan estén suficientemente lejos y si, por el contrario, los objetos son cercanos sus representantes también lo sean.

La solución se basa en un trabajo de Young y Householder (1938) y hay que decir que no es única. El conjunto de coordenadas cuyas distancias se aproximan a las distancias originales se puede trasladar, rotar o reflectar, sin que las distancias se distorsionen, de manera que no se pueden determinar completamente. El problema de localización se resuelve tomando el punto con las medias como el origen de coordenadas. En cuanto a las rotaciones, el algoritmo proporciona una, pero cualquier otra puede servir para facilitar la interpretación de las soluciones.

Algunos detalles matemáticos

Supongamos que la matriz de datos \(\mathbf{X}\) sea tal que las distancias euclídeas entre sus filas sea la solución que buscamos. Dada una matriz como esta, calcular las distancias euclídeas es fácil. El problema ahora es justamente el contrario: Dada una matriz de distancias \(\mathbf{D}\) \(n\times n\), ¿como hallamos \(\mathbf{X}\)?

Para empezar, definimos la matriz \(\mathbf{B}\) \(n\times n\) como \[\mathbf{B} = \mathbf{X} \mathbf{X}'\] Los elementos de \(\mathbf{B}\) son \[b_{ij} = \sum_{k=1}^q x_{ik}x_{jk}\] Se puede comprobar que el cuadrado de las distancias euclídeas entre las filas de \(\mathbf{X}\) puede escribirse en términos de los elementos de \(\mathbf{B}\) como \[d_{ij}^2 = b_{ii} + b_{jj} - 2b_{ij}\] Si en lugar de las distancias euclídeas ponemos las distancias de la matriz \(\mathbf{D}\), entonces tenemos un sistema de ecuaciones para resolver. Dicho sistema de ecuaciones no es trivial, pero con algunas restricciones los elementos de \(\mathbf{B}\) se pueden hallar en función de las distancias: \[b_{ij} = -\frac 12\left[ d_{ij}^2 - d_{i\cdot}^2 - d_{\cdot j}^2 + d_{\cdot\cdot}^2\right]\] donde \[d_{i\cdot}^2 = \frac 1n \sum_{j=1}^n d_{ij}^2\]
\[d_{\cdot j}^2 = \frac 1n \sum_{i=1}^n d_{ij}^2\] \[d_{\cdot\cdot}^2 = \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n d_{ij}^2\] Una vez hallados los elementos de \(\mathbf{B}\), sólo nos queda factorizar para obtener las coordenadas \(\mathbf{X}\). La descomposición en valores singulares de \(\mathbf{B}\) se puede escribir \[\mathbf{B} = \mathbf{U} \mathbf{\Lambda} \mathbf{U}'\] donde \(\mathbf{\Lambda} = \text{diag}(\lambda_1,\dots,\lambda_n)\) es la matriz diagonal con los valores propios de \(\mathbf{B}\) y \(\mathbf{U}\) sus correspondientes vectores propios normalizados de forma que la suma de los cuadrados de los elementos de cada columna es 1. Los valores propios se ordenan de mayor a menor de forma que \(\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_n\).

Cuando \(\mathbf{D}\) proviene de una matriz \(n\times q\) de rango máximo, entonces la matriz \(\mathbf{B}\) es de rango \(q\) y sus últimos \(n-q\) valores propios son cero. Así que \(\mathbf{B}\) se puede escribir como \[\mathbf{B} = \mathbf{U}_q \mathbf{\Lambda}_q \mathbf{U}_q'\] donde \(\mathbf{\Lambda}_q\) contiene los primeros \(q\) valores propios y \(\mathbf{U}_q\) sus correspondientes vectores propios.

Las coordenadas de los puntos que buscamos son \[\mathbf{X} = \mathbf{U}_q \mathbf{\Lambda}^{1/2}_q\] donde \(\mathbf{\Lambda}^{1/2}_q= \text{diag}(\lambda^{1/2}_1,\dots,\lambda^{1/2}_q)\).

La mejor representación de dimensión \(k\) viene dada por los primeros \(k\) valores propios de \(\mathbf{B}\) y sus correspondientes vectores propios.

El ajuste de la representación en dimensión \(k\) a los datos originales se puede medir con el criterio \[P_k = \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^q \lambda_i}\] Valores de \(P_k\) superiores a \(0.8\) sugieren un buen ajuste.

Cuando la matriz de disimilaridades no es euclídea, la matriz \(\mathbf{B}\) no es definida positiva. En estos casos algunos de los valores propios de \(\mathbf{B}\) serán negativos y, en consecuencia, algunas coordenadas serán números complejos. Sin embargo, si \(\mathbf{B}\) tiene un número pequeño de valores propios negativos y éstos son pequeños, la representación asociada a los \(k\) mayores valores propios positivos todavía es válida. En estos casos el criterio de ajuste puede ser \[P_k^{(1)} = \frac{ \sum_{i=1}^k \text{abs}(\lambda_i)}{ \sum_{i=1}^n \text{abs}(\lambda_i)}\] \[P_k^{(2)} = \dfrac{\sum_{i=1}^k \lambda^2_i}{\sum_{i=1}^n \lambda^2_i}\]

Algunas recomendaciones alternativas son:

1. Criterio de la traza: Elegir el número de coordenadas de forma que la suma de las que corresponden a valores propios positivos sumen aproximadamente los mismo que la suma de todos los valores propios (positivos y negativos).
2. Criterio de la magnitud: Aceptar como genuinamente positivos los valores propios cuya magnitud exceda subtancialmente la del mayor valor propio negativo.

Un ejemplo

Consideremos las distancias entre cinco objetos

d <- matrix(c(0, 1, 2, 4, 5, 
              1, 0, 2.5, 4, 5, 
              2, 2.5, 0, 4, 5, 
              4, 4, 4, 0, 1.2, 
              5, 5, 5, 1.2, 0), nrow = 5)
rownames(d) <- colnames(d) <- c("A", "B", "C", "D", "E")
d

##   A   B   C   D   E
## A 0 1.0 2.0 4.0 5.0
## B 1 0.0 2.5 4.0 5.0
## C 2 2.5 0.0 4.0 5.0
## D 4 4.0 4.0 0.0 1.2
## E 5 5.0 5.0 1.2 0.0

El escalado multidimensional se obtiene con la función cmdscale() de R

mds <- cmdscale(d)
mds

##        [,1]        [,2]
## A -1.837444  0.42826525
## B -1.757255  1.00509080
## C -1.620166 -1.47854596
## D  2.086176 -0.01557718
## E  3.128689  0.06076709

Las distancias obtenidas del MDS son

dist(mds)

##           A         B         C         D
## B 0.5823727                              
## C 1.9191505 2.4874173                    
## D 3.9486436 3.9766471 3.9846260          
## E 4.9797124 4.9763642 4.9921049 1.0453055

Finalmente, podemos representar los objetos como puntos de un mapa de dimensión dos:

plot(mds, type = "n", xlab = "Coord. 1", ylab = "Coord. 2")
text(mds[, 1], mds[, 2], labels = rownames(mds), cex = 0.8)

Los dos primeros valores propios son

mds <- cmdscale(d, eig = T)
mds$eig[1:2]

## [1] 23.229909  3.383652

y el ajuste

mds$GOF

## [1] 0.9806751 0.9806751

Bibliografía

[1] Everitt, Brian (2005), An Rand S-PLUS@ companion to multivariate analysis. Springer texts in statistics.