Definición general

Los datos, a menudo pueden representarse, además de la matriz/tabla de registros y campos,también como matriz de proximidades (englobando matriz de distancias y matriz de similaridades).

El Escalado Multidimensional pretende descubrir estructuras que pudiera esconder la matriz de proximidades, deteminando:


Escalado Multidimensional Métrico Euclídeo y no Euclídeo

Objetivo: Reconstruir la matriz de datos a partir la matriz de proximidades dada.

Reconstrucción de la matriz de datos a partir de la matriz de distancias

Dada una matriz de distancias \(\textbf{D}=(d_{ij})_{n\times n}\), pretendemos obtener la matriz de datos \(\textbf{X}_{n\times p}\).

La matriz de datos \(\text{X}\)

\[ \textbf{X}= \begin{array} {|cc|} x_{11} & x_{1j} & x_{1p} \\ x_{i1} & x_{ij} & x_{ip} \\ x_{n1} & x_{nj} & x_{np} \end{array} \]

La matriz \(\textbf{B}\) generada de \(\textbf X\)

\[ \textbf B=\textbf X\textbf X^t \]

La relación de la matriz \(\textbf B\) con la matriz de distancias \(\textbf D\)

Se puede demostrar que \(\textbf B\) está relacionada con \(\textbf D\) por la función:

\[ b_{ij}=-\frac{1}{2}(d_{ij}^2-\frac{1}{n}\sum_{1=j}^nd_{ij}^2- \frac{1}{n}\sum_{1=i}^nd_{ij}^2 +\frac{1}{n^2}\sum_{1=j}\sum_{1=i}^nd_{ij}^2 ) \]

La descomposición espectral de \(\textbf B\)

\[ \textbf B = \textbf{VΔV}^t \]

\(\textbf{Δ}\) es la matriz diagonal de autovalores (eigenvalues) de \(\textbf{B}\), \(\lambda_1\geq\lambda_2\geq ...\geq \lambda_n\).

\(\textbf V\) es la matriz de autovectores de \(\textbf B\).

Con una matriz inicial \(\textbf X\) de rango completo, el rango de \(\textbf B\) sería \(p\) y sus últimos \(n-p\) autovalores serían \(0\) y \(\textbf B\) se puede escribir como

\[ \textbf B =\textbf U_1Δ_1\textbf{U}_1^t \]

\(\textbf U_1\) es la matriz que contiene los primeros autovectores de \(\textbf B\).

\(\textbf{Δ}_1\) son autovalores de \(\textbf B\) distintos de \(0\).

Lo anterior quedará finalmente así:

\[ \textbf X =\textbf U_1\textbf{Δ}_1^{1/2} \]

Donde \(\textbf{Δ}_1^{1/2}\) es una matriz diagonal de valores \(\sqrt\lambda_1\geq\sqrt\lambda_2\geq ...\geq \sqrt\lambda_n\).

Estas coordenadas que se determinan para \(\textbf X\) son los scores que se hubieran obtenido en un análisis de componentes principales de \(\textbf X\). Por ello el Escalado Multidimensional también se puede denominar Análisis de Coordenadas Principales.


Matriz de Proximidades Euclídea y no Euclídea

Si la matriz \(\textbf B\) obtenida a partir de \(\textbf D\) es euclídea sus \(n-1\) primeros autovalores \(\lambda_i\) serán positivos y el último será \(0\). Decidiremos representar los datos en un espacio de \(k\) dimensiones si el cociente

\[ \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^{n-1}\lambda_i} \]

será suficientemente grande (0.8 o más).

En muchos casos, la matriz de proximidades será no Euclídea y así los autovalores no serán todos positivos. Los últimos autovalores de \(\textbf B\) serán negativos. Sin embargo, si son pocos y pequeños, se pueden utilizar los autovectores de los primeros autovalores como coordenadas.

Criterios para elección de la dimensión \(k\):

  1. \[ \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^{n}|\lambda_i|} \] es suficientemente grande, o

  2. \[ \frac{\sum_{i=1}^k \lambda_i^2}{\sum_{i=1}^{n}\lambda_i^2} \] es suficientemente grande.


Teoría según Método Nottingham

Ver R. Wilkinson, Applied multivariate statistics 6.1, University of Nottingham, 2026.

En el MDS clásico, no trabajamos directamente con la matriz de distancia, sino con la matriz de producto interno centrado.

La matriz de producto interno centrado

Definición

Dada una matriz de distancias \(\textbf D =\{d_{ij}\}_{i,j=1}^n\), la matriz de producto interno centrado (también llamada matriz de Gram centrada) es

\[ \begin{equation} {\mathbf B}={\mathbf H} \mathbf A{\mathbf H} \end{equation} \]

Matriz de distancias cuadradas negativas

donde \(\mathbf A\) es la matriz de distancias cuadradas negativas dividida por dos:

\[ \begin{equation} \mathbf A=\{a_{ij}\}_{i,j=1}^n, \quad \text{donde} \qquad a_{ij}=-\frac{1}{2}d_{ij}^2 \end{equation} \]

Matriz de centrado

y \(\mathbf H\) es la matriz de centrado \(n × n\) (ver 2.4 del mismo manual).

\[ \begin{equation} \mathbf H=\mathbf I_n - \frac{1}{n} {\mathbf 1}_n {\mathbf 1}_n^\top. \tag{2.2} \end{equation} \]

Matriz de distancias cuadradas negativas - otra definición

Otra manera de escribir la ecuación es

\[ \mathbf A= -\frac{1}{2}\mathbf D\odot \mathbf D \]

donde \(\odot\) denota el producto Hadamard (o producto elemento a elemento de dos matrices (que sería un producto no matricial).

La ecuación se puede mejor presentar de esta forma \[ \mathbf B=-\frac{1}{2}\mathbf H(\mathbf D\odot\mathbf D)\mathbf H \]

Teorema

Siendo \(\mathbf D\) una matriz de distancia \(n×n\) con una matriz de producto interno centrada correspondiente \(\mathbf B=-\frac{1}{2}\mathbf H(\mathbf D\odot\mathbf D)\mathbf H\),

  1. Si \(\textbf D\) es una matriz de distancia euclidiana para la muestra de \(n\) vectores \(x_1,x_2,...x_n\), entonces

\[ \begin{equation} \mathbf B= ({\mathbf H} {\mathbf X})({\mathbf H} {\mathbf X})^\top \end{equation} \] Así, \(\textbf B\) es positiva semi-definida.

  1. Suponiendo que que \(\textbf B\) es semidefinida positiva con valores propios \(\lambda_1 \geq \lambda_2 \cdots \geq \lambda_k > 0\) y que tiene descomposición espectral \(\mathbf B={\mathbf U} {\pmb \Lambda}{\mathbf U}^\top\), donde \({\pmb \Lambda}=\text{diag}\{\lambda_1 \ldots \lambda_k\}\) y \(\textbf U\), es \(n × k\) y satisface \({\mathbf U}^\top {\mathbf U}={\mathbf I}_k\). Entonces \[ {\mathbf X}=[\mathbf x_1, \ldots , \mathbf x_n]^\top={\mathbf U}{\pmb \Lambda}^{1/2} \] es una matriz de datos \(n \times k\) con una matriz de distancias euclídeas \(\textbf D\).


La matriz de producto interno centrado con R

Datos

X <- matrix(c(0,0,
              1,0,
              0,1,
              -1,0,
              0,-1), nc=2, byrow=TRUE)
X
##      [,1] [,2]
## [1,]    0    0
## [2,]    1    0
## [3,]    0    1
## [4,]   -1    0
## [5,]    0   -1

Matriz de distancas

# Matriz de distancias
D <- as.matrix(dist(X, upper=T, diag=T))
D
##   1        2        3        4        5
## 1 0 1.000000 1.000000 1.000000 1.000000
## 2 1 0.000000 1.414214 2.000000 1.414214
## 3 1 1.414214 0.000000 1.414214 2.000000
## 4 1 2.000000 1.414214 0.000000 1.414214
## 5 1 1.414214 2.000000 1.414214 0.000000

Matriz de distancias cuadradas negativas

# A: matriz de distancias cuadradas negativas

A <- - D*D/2 # Producto negativo elemento a elemento, dividido por 2
A
##      1    2    3    4    5
## 1  0.0 -0.5 -0.5 -0.5 -0.5
## 2 -0.5  0.0 -1.0 -2.0 -1.0
## 3 -0.5 -1.0  0.0 -1.0 -2.0
## 4 -0.5 -2.0 -1.0  0.0 -1.0
## 5 -0.5 -1.0 -2.0 -1.0  0.0

Matriz de centrado

# H: Matriz de centrado
H <- diag(5) - 1/5 * matrix(rep(1,5), nc=1)%*%matrix(rep(1,5), nr=1)
H
##      [,1] [,2] [,3] [,4] [,5]
## [1,]  0.8 -0.2 -0.2 -0.2 -0.2
## [2,] -0.2  0.8 -0.2 -0.2 -0.2
## [3,] -0.2 -0.2  0.8 -0.2 -0.2
## [4,] -0.2 -0.2 -0.2  0.8 -0.2
## [5,] -0.2 -0.2 -0.2 -0.2  0.8

Matriz de producto interno centrado

B <- H %*% A %*% H
B
##               [,1]          [,2]          [,3]          [,4]          [,5]
## [1,] -2.775558e-17  2.775558e-17  2.775558e-17  4.163336e-17  5.551115e-17
## [2,] -3.122502e-17  1.000000e+00 -2.255141e-16 -1.000000e+00 -2.081668e-16
## [3,]  2.775558e-17 -1.387779e-16  1.000000e+00 -1.665335e-16 -1.000000e+00
## [4,]  5.204170e-17 -1.000000e+00 -1.422473e-16  1.000000e+00 -1.526557e-16
## [5,]  8.326673e-17 -1.387779e-16 -1.000000e+00 -1.665335e-16  1.000000e+00

Control

Redondeando a 4 decimales la matriz \(\textbf B\) obtenida (R prouce valores infinitesimales de redondeo), obtenemos la siguiente

#check this matches 
round(B,4)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    0
## [2,]    0    1    0   -1    0
## [3,]    0    0    1    0   -1
## [4,]    0   -1    0    1    0
## [5,]    0    0   -1    0    1

Según el Teorema es

\[ \begin{equation} \mathbf B= ({\mathbf H} {\mathbf X})({\mathbf H} {\mathbf X})^\top \end{equation} \]

Lo que queda, confirmado por los datos de R:

(H %*% X) %*% t(H %*% X)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    0
## [2,]    0    1    0   -1    0
## [3,]    0    0    1    0   -1
## [4,]    0   -1    0    1    0
## [5,]    0    0   -1    0    1


Escalado Multidimensional con R

La obtención de la matriz de datos X a partir de la matriz de proximidades puede puede hacers con la función \[\texttt{cmdscale(d,k,eig)} \]

d: Matriz de proximidades.

k: dimensiones del espacio en que representar los datos (2, por defecto).

eig: Es la opción de obtener o no (T/F) los autovalores mencionados.


Ejemplos

Ejemplo 1 de obtención de coordenadas a a partir de la función \(\texttt{cmdscale()}\)

Datos obtenidos en R. Wilkinson, Applied multivariate statistics 6., University of Nottingham, 2026.

Datos y manipulaciones previas

library(dplyr)

dfcities <- as.data.frame(readxl::read_xlsx("EAAR/UKcitydist.xlsx"))

rownames(dfcities) <- dfcities$...1

dmatrix <- dfcities %>% 
  select(-1) %>% 
  as.matrix()
dmatrix
##             London Birmingham Manchester Leeds Newcastle Liverpool Portsmouth
## London           0        163        262   273       403       286        103
## Birmingham     163          0        114   149       282       125        195
## Manchester     262        114          0    58       174        50        308
## Leeds          273        149         58     0       135       105        335
## Newcastle      403        282        174   135         0       199        469
## Liverpool      286        125         50   105       199         0        317
## Portsmouth     103        195        308   335       469       317          0
## Southampton    112        179        293   323       458       299         24
## Nottingham     175         73         94    98       231       132        239
## Bristol        170        124        227   271       401       219        127
## Sheffield      228        105         53    47       181       101        288
## Norwich        159        217        255   230       328       299        261
## Plymouth       308        281        369   420       542       346        221
##             Southampton Nottingham Bristol Sheffield Norwich Plymouth
## London              112        175     170       228     159      308
## Birmingham          179         73     124       105     217      281
## Manchester          293         94     227        53     255      369
## Leeds               323         98     271        47     230      420
## Newcastle           458        231     401       181     328      542
## Liverpool           299        132     219       101     299      346
## Portsmouth           24        239     127       288     261      221
## Southampton           0        229     103       276     268      202
## Nottingham          229          0     193        53     169      353
## Bristol             103        193       0       228     296      162
## Sheffield           276         53     228         0     203      382
## Norwich             268        169     296       203       0      453
## Plymouth            202        353     162       382     453        0

Análisis de autovalores con \(\texttt{cmdscale()}\)

cmdscale(dmatrix, eig = T)$eig
##  [1]  2.832068e+05  9.662516e+04  2.996738e+02  2.558621e+02  1.423145e+02
##  [6]  9.662094e+01  5.613665e+01  3.410386e+01 -8.753887e-12 -6.758398e+01
## [11] -1.780337e+02 -3.042158e+02 -3.431856e+02

Como podemos ver, en las últimas posiciones existen valores negativos, lo que nos indica que las proximidades no son perfectamente euclídeas.

Vamos pues a analizar los cocientes vistos arriba, para \(k=2\):

\[ \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^{n-1}|\lambda_i|} \]

sum(cmdscale(dmatrix, eig = T)$eig[1:2])/sum(abs(cmdscale(dmatrix, eig = T)$eig)[-13])
## [1] 0.9962374

\[ \frac{\sum_{i=1}^k \lambda_i^2}{\sum_{i=1}^{n-1}\lambda_i^2} \]

sum(cmdscale(dmatrix, eig = T)$eig[1:2]^2)/sum(abs(cmdscale(dmatrix, eig = T)$eig^2)[-13])
## [1] 0.9999964

En ambos casos, se puede apreciar como las \(k=2\) dimensiones capturan una fracción altísima de la suma de autovalores. Es por lo tanto viable la representación bidimensional de las \(k\) coordenadas.

Representación gráfica de las coordenadas

plot(cmdscale(dmatrix)[,1],cmdscale(dmatrix)[,2],
     xlab = "Primera coordenada",ylab = "Segunda coordenada")
text(cmdscale(dmatrix)[,1], cmdscale(dmatrix)[,2],rownames(cmdscale(dmatrix)))

El mapa es substancialmente correcto pero está rotado de \(90^o\) y las coordenadas tiene promedio \(0\).

Ejemplo 2

Datos obtenidos de EAAR, capítulo 41.

Se trata de datos de proximidad de 12 países generados por encuesta en España:

paises <- c("C. de Marfil", "Cuba", "Venezuela", "Egipto", "UK", "India", 
            "Israel", "Japón", "China", "Rusia", "USA", "España")
matrix(paises)
##       [,1]          
##  [1,] "C. de Marfil"
##  [2,] "Cuba"        
##  [3,] "Venezuela"   
##  [4,] "Egipto"      
##  [5,] "UK"          
##  [6,] "India"       
##  [7,] "Israel"      
##  [8,] "Japón"       
##  [9,] "China"       
## [10,] "Rusia"       
## [11,] "USA"         
## [12,] "España"

Matriz de distancias

p <- matrix(scan("datosr/paises.txt"),ncol = 12)
p
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,]  0.0  7.5  8.5  6.0  9.5  7.0  8.0  9.5  5.5   8.5   9.5   9.0
##  [2,]  7.5  0.0  2.0  5.0  7.5  5.0  2.0  9.0  5.5   6.0   9.0   6.0
##  [3,]  8.5  2.0  0.0  7.5  8.0  8.5  2.0  9.0  9.0   8.0   7.5   5.5
##  [4,]  6.0  5.0  7.5  0.0  8.0  5.5  7.5  7.0  7.5   8.0   8.0   9.0
##  [5,]  9.5  7.5  8.0  8.0  0.0  6.0  2.0  5.0  7.5   6.5   1.5   2.0
##  [6,]  7.0  5.0  8.5  5.5  6.0  0.0  7.5  7.5  6.0   7.5   7.5   9.0
##  [7,]  8.0  2.0  2.0  7.5  2.0  7.5  0.0  2.0  2.0   2.5   3.0   0.5
##  [8,]  9.5  9.0  9.0  7.0  5.0  7.5  2.0  0.0  2.5   6.0   2.0   4.5
##  [9,]  5.5  5.5  9.0  7.5  7.5  6.0  2.0  2.5  0.0   4.5   6.5   6.5
## [10,]  8.5  6.0  8.0  8.0  6.5  7.5  2.5  6.0  4.5   0.0   5.0   6.0
## [11,]  9.5  9.0  7.5  8.0  1.5  7.5  3.0  2.0  6.5   5.0   0.0   2.5
## [12,]  9.0  6.0  5.5  9.0  2.0  9.0  0.5  4.5  6.5   6.0   2.5   0.0

Análisis de los autovalores

cmdscale(p,eig = T)$eig
##  [1]  1.101409e+02  7.176177e+01  4.345015e+01  2.937512e+01  2.148003e+01
##  [6]  1.501705e+01  4.592181e+00  1.776357e-14 -3.241329e+00 -1.007454e+01
## [11] -1.998470e+01 -2.437084e+01
  1. \[ \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^{n-1}|\lambda_i|} \]
sum(cmdscale(p, eig = T)$eig[1:2])/sum(abs(cmdscale(p, eig = T)$eig))
## [1] 0.5145928
  1. \[ \frac{\sum_{i=1}^k \lambda_i^2}{\sum_{i=1}^{n-1}\lambda_i^2} \]
sum(cmdscale(p, eig = T)$eig[1:2]^2)/sum(abs(cmdscale(p, eig = T)$eig^2))
## [1] 0.7910665

Tal y como se puede observar, el primer cociente es bastante bajo. Sin embargo, el segundo es suficiente para poder representar las \(k=2\) dimensiones.

Representación gráfica de las coordenadas

plot(cmdscale(p)[,1],cmdscale(p)[,2],
     xlab = "Primera coordenada",ylab = "Segunda coordenada")
text(cmdscale(p)[,1], cmdscale(p)[,2])

Representación gráfica de las coordenadas con \(\texttt{ggplot2}\)

library(ggplot2)

rownames(p) <- paises


p1 <- data.frame(cmdscale(p)) %>% 
  tibble::rownames_to_column("País")

p1 %>% 
  ggplot(aes(x=X1, y=X2)) +
  geom_point() +
  ggrepel::geom_text_repel(aes(label= País)) +
  xlab("Dimensión 1") +
  ylab("Dimensión 2")

Escalado Multidimensional no Métrico

Hasta este punto, la matriz de proximidades tenía un valor métrico: una distancia 2 se entendía como el doble de una distancia 1. En muchos casos (ciencias del comportamiento), en los que se trata de juicios de valor, la matriz representa ya no distancias sino ordenaciones.

En este caso se trabajará con la función de la librería \(MASS\) \(,\texttt{isoMDS(d,k)}\). \(\texttt{d}\) es la matriz, \(\texttt{k}\) es el número de dimensiones en el espacio en que representar los puntos.

Una medida de la bondad de la \(k\) elegida, es \(stress\).

Ejemplo 3 (sobre datos de ejemplo 2)

Ejecución de la función

library(MASS)

isoMDS(p)
## initial  value 22.013391 
## iter   5 value 16.278431
## final  value 16.235616 
## converged
## $points
##                    [,1]       [,2]
## C. de Marfil -6.4218745  1.8066001
## Cuba         -2.0913650 -1.8517701
## Venezuela    -0.6750485 -5.1552200
## Egipto       -4.6640368 -0.5992696
## UK            3.8877607  0.1376858
## India        -2.6655991  3.0347480
## Israel        1.1433001 -1.5836915
## Japón         3.2008829  2.8951639
## China        -0.6211791  2.5329630
## Rusia         1.8665246  1.5309132
## USA           4.7873080 -0.4208432
## España        2.2533267 -2.3272795
## 
## $stress
## [1] 16.23562

Como se puede ver, la función nos facilita un dato de stress razonable (16.2) y las coordenadas que se pueden representar.

Representación gráfica

p2 <- data.frame(isoMDS(p)$points) %>% 
  tibble::rownames_to_column("País")
## initial  value 22.013391 
## iter   5 value 16.278431
## final  value 16.235616 
## converged
p2 %>% 
  ggplot(aes(x=X1, y=X2)) +
  geom_point() +
  ggrepel::geom_text_repel(aes(label= País)) +
  xlab("Dimensión 1") +
  ylab("Dimensión 2")

Con una representación bastante similar a la realizada a partir de la función cmdscale().





  1. Siempre que citamos EAAR en este documento, nos referimos a: A. García Perez, Estadística aplicada avanzada con R, Madrid, 2022↩︎