multidimensional scaling of variables

dimension reduction of geometric data

Because ordination is grounded in geometric data analysis, it is easiest and most straightforward on naturally geometric data. For example, the mineral compositional measurements analyzed by Freestone, Gorin–Rosen, and Hughes (2000) consists of five primary (and several secondary) concentrations for each of 68 artifacts. While all the variables lie on the same (percentage) scale, they vary at different absolute rates; in order to represent the cases (artifacts) geometrically we might separately center and scale them. The most basic ordination method, principal components analysis (PCA), computes a singular value decomposition \(X = U D V^\top\) of the centered and scaled data matrix \(X \in \mathbb{R}^{n\times p}\):

data(glass, package = "ordr")
x <- scale(glass[, c("SiO2", "Al2O3", "FeO", "MgO", "CaO")],
           center = TRUE, scale = TRUE)
s <- svd(x)
r <- length(s$d)

The matrix factors \(U \in \mathbb{R}^{n\times r}\) and \(V \in \mathbb{R}^{p\times r}\) arise from eigendecompositions of \(X X^\top\) and of \(X^\top X\), respectively, which have the same set of eigenvalues \(\lambda_1,\ldots,\lambda_r\). The square roots of these eigenvalues make up the diagonal matrix \(D \in \mathbb{R}^{r\times r}\). The conventional PCA biplot uses principal coordinates for the cases, given by the rows of \(U D\), and standardized coordinates for the variables, given by the rows of \(V\).

The columns of \(U = \left[\,u_1\,\cdots\,u_r\,\right]\) and \(V = \left[\,v_1\,\cdots\,v_r\,\right]\) comprise the eigenvectors of \(X X^\top\) and of \(X^\top X\). They are orthonormal, which means that \(U^\top U = I_r = V^\top V\) and that the total variance (called inertia) in each matrix is \(\sum_{j=1}^{r}{ {v_j}^2 } = r = \sum_{j=1}^{r}{ {v_j}^2 }\). Meanwhile, the diagonal \(D\) contains all of the inertia of the centered and scaled data matrix \(X\):

# inertia of the (scaled) data
sum(x^2)

## [1] 335

# inertia of the case and variable factors
sum(s$u^2)

## [1] 5

sum(s$v^2)

## [1] 5

# inertia of the diagonal factor
sum(s$d^2)

## [1] 335

For the (2-dimensional, i.e. planar) biplot, this inertia is conferred on the cases (converting standardized coordinates to principal coordinates) so that the planar distances between the points in the plot representing the cases are approximately equal to their Euclidean distances in \(X\). Because the biplot only accommodates two principal components (it is a projection onto a subset of orthogonal coordinates), the planar distances underestimate the true distances. Here is a scatterplot of the true distances and their biplot approximations:

# distances between cases
x.dist <- dist(x)
# distances between cases (principal coordinates)
s.dist <- dist(s$u[, 1:2] %*% diag(s$d[1:2]))
# scatterplot
plot(
  x = as.vector(x.dist),
  y = as.vector(s.dist),
  xlim = c(0, 10), ylim = c(0, 10),
  asp = 1, pch = 19, cex = .5,
  xlab = "Case distances in centered and scaled data",
  ylab = "Case point distances in planar biplot"
)
lines(x = c(0, 10), y = c(0, 10))

In contrast, the variables are better understood through their correlations, which are approximately preserved by their standardized coordinates. Writing \(X = [\,y_1\,\cdots\,y_m\,]\) as an array of column variables, the covariance between \(y_i\) and \(y_j\) is proportional to their inner product \[\textstyle \operatorname{cov}(y_i,y_j) = \frac{1}{n} y_i \cdot y_j = \frac{1}{n} \lVert y_i\rVert\lVert y_j\rVert\cos\theta_{ij}\text,\] so that the cosine of the angle \(\theta_{ij}\) between them equals their correlation: \[\cos\theta_{ij} = \frac{\operatorname{cov}(y_i,y_j)}{\sqrt{\operatorname{cov}(y_i,y_i)\operatorname{cov}(y_j,y_j)}/n} = \frac{\operatorname{cov}(y_i,y_j)}{\sigma_i\sigma_j} = r_{ij}\]

Here the cosines between the vectors in the biplot are plotted against the variable correlations \(r_{ij}\) in the centered and scaled data. While the vectors are shorter in the biplot for the same reason that distances are shorter, their consines may be larger or smaller depending on their respective coordinates in the remaining coordinates.

# correlations between variables
x.cor <- cor(x)
# magnitudes of variable vectors
s.len <- apply(s$v[, 1:2] %*% diag(s$d[1:2]), 1, norm, "2")
# cosines between variables (principal coordinates)
s.cor <- (s$v[, 1:2] / s.len) %*% diag(s$d[1:2]^2) %*% t(s$v[, 1:2] / s.len)
# scatterplot
plot(
  x = as.vector(x.cor[lower.tri(x.cor)]),
  y = as.vector(s.cor[lower.tri(s.cor)]),
  xlim = c(-1, 1), ylim = c(-1, 1),
  asp = 1, pch = 19, cex = .5,
  xlab = "Variable correlations in centered and scaled data",
  ylab = "Variable vector cosines in planar biplot"
)
lines(x = c(-1, 1), y = c(-1, 1))

multidimensional scaling of distance data

The faithful approximation of inter-case distances by principal coordinates is the idea behind classical multidimensional scaling (CMDS), which can be applied to a data set of distances \(\delta_{ij},\ 1\leq i\leq j\leq n\) in the absence of coordinates (variables). CMDS produces a set of artificial coordinates for the cases that yield nested best approximations of the inter-case distances in terms of the sum of squared errors. The technique uses the eigendecomposition of a doubly-centered matrix of squared distances, which produces a matrix \(U \Lambda^{1/2}\) whose first \(r\) coordinates—for any \(r\leq n\)—minimize the variance of \((U \Lambda^{1/2}) (U \Lambda^{1/2})^\top - \Delta = U \Lambda U^\top - \Delta\), where \(\Delta = (\delta_{ij}) \in \mathbb{R}^{n\times n}\), to obtain the approximation. The columns of \(U \Lambda^{1/2}\) are then called principal coordinates, and the first \(r\) most faithfully recover the \(\delta_{ij}\) using points in \(\mathbb{R}^r\) of any possible \(r\). For example, because the road distances between several U.S. cities arise from a roughly 2-dimensional process, the point distances in a CMDS are very close approximations:

d <- as.matrix(UScitiesD)
cent <- diag(1, nrow(d)) - matrix(1/nrow(d), nrow(d), nrow(d))
d.cent <- -.5 * cent %*% (d^2) %*% cent
d.cmds <- svd(d.cent)
d.coord <- d.cmds$u[, 1:2] %*% diag(sqrt(d.cmds$d[1:2]))
# scatterplot
plot(
  x = as.vector(UScitiesD),
  y = as.vector(dist(d.coord)),
  xlim = c(0, max(UScitiesD)), ylim = c(0, max(UScitiesD)),
  asp = 1, pch = 19, cex = .5,
  xlab = "City road distances",
  ylab = "Point distances in planar CMDS"
)
lines(x = c(0, max(UScitiesD)), y = c(0, max(UScitiesD)))

In practice, the goal of CMDS is usually to produce a scatterplot in which the distances \(\sqrt{(x_j-x_i)^2+(y_j-y_i)^2}\) between the points that represent the \(n\) cases approximate their original distances \(\delta_{ij}\). In this case, the artificial coordinates approximately recover the geographic arrangement. By chance, this CMDS rotated the conventional cardinal directions by about \(\pi\) radians (though at least the map is recognizable from above rather than below the surface of the Earth):

plot(
  d.coord, pch = NA, asp = 1,
  xlab = "First principal coordinate", ylab = "Second principal coordinate"
)
text(d.coord, labels = rownames(d), cex = .9)

multidimensional scaling of covariance data

The faithful approximation of inter-variable correlations by the angles between their coordinate vectors provides a dual CMDS procedure. Suppose we have data that consist not of distances between cases but of covariances \(\operatorname{cov}(y_i,y_j),\ 1\leq i\leq j\leq p\) between variables. Again the data are coordinate-free, so PCA is inapplicable. But were the data derived from a (not necessarily centered or scaled) case–variable matrix \(X\), then the covariance matrix \(C=(\operatorname{cov}(y_i,y_j))\) would have been obtained as \(C=\frac{1}{n}X^\top X\). This, up to scalar, is the matrix whose eigenvectors would be given by \(V\) in the SVD \(X = U D V^\top\). Therefore, we can obtain artificial coordinates for these variables that approximate what we know of their geometry—thinking of the variables as unknown vectors whose magnitudes and angles are encoded in \(C\)—via an eigendecomposition \(C = V \Lambda V^\top\): Take \(Y = V \Lambda^{1/2} \in \mathbb{R}^{p\times r}\), so that \(Y^\top Y \approx C\).

While covariances of the unscaled glass measurements are not meaningful, they can be used to validate the technique. Again because the eigendecomposition is intrinsically ordered by variance, the first \(r\) eigenvectors provide the most faithful \(r\)-dimensional approximation; i’ll take \(r=2\) in anticipation of a biplot:

# covariances and standard deviations
c <- cov(x)
s <- diag(sqrt(diag(c)))
# eigendecomposition of covariance matrix
c.eigen <- eigen(c)
# artificial coordinates
c.coord <- c.eigen$vectors[, 1:2] %*% diag(sqrt(c.eigen$values[1:2]))
# scatterplot
c.inner <- c.coord %*% t(c.coord)
plot(
  x = as.vector(c[lower.tri(c)]),
  y = as.vector(c.inner[lower.tri(c.inner)]),
  xlim = range(c[lower.tri(c)]), ylim = range(c[lower.tri(c)]),
  asp = 1, pch = 19, cex = .5,
  xlab = "Measurement covariances in unscaled data",
  ylab = "Vector inner products in planar CMDS"
)
lines(x = range(c[lower.tri(c)]), y = range(c[lower.tri(c)]))

Thus, whereas CMDS of cases approximates distances, CMDS of variables approximates covariances. This equips us with an analogous visualization technique for correlation data: Represent each variable \(y_i\) as a unit vector \(\hat{y}_i\) in such a way that their angle cosines \(\cos\theta_{ij}\) equal their correlations \(r_{ij}\), then project these vectors onto the plane in which their variance is maximized. This is exactly how variables are commonly represented in row-principal PCA biplots, by the equivalence of \(V\) in the SVD \(X = U D V^\top\) with \(E\) in the eigendecomposition \(X^\top X = E \Lambda E^\top\). Here is the result for the five compositional measurements on the glass artifacts:

c <- cor(glass[, c("SiO2", "Al2O3", "FeO", "MgO", "CaO")])
c.eigen <- eigen(c)
c.coord <- c.eigen$vectors[, 1:2] %*% diag(sqrt(c.eigen$values[1:2]))
plot(
  c.coord, pch = NA, asp = 1,
  xlab = "First principal coordinate", ylab = "Second principal coordinate"
)
arrows(0, 0, c.coord[, 1L], c.coord[, 2L])
text(c.coord, labels = rownames(c), cex = .9)

We can immediately see that the concentrations of silica and ferrous oxide are highly correlated, that they are anticorrelated with that of quicklime, and that all of these vary independently (in the statistical sense) of those of aluminium oxide and magnesia, which are anticorrelated with each other.

In this setting, we could have obtained this visualization directly from the case–variable data \(X\), as in PCA. The technique becomes uniquely useful when we have correlations without underlying coordinates.

use case: rankings of universities

library(ordr)
library(dplyr)

A natural use case for CMDS of variables is the analysis of multiple rankings of the same set of objects in terms of their concordance. Rankings’ concordance is often measured using rank correlations such as Kendall’s \(\tau\), which may be general correlation coefficients in the sense proposed by Kendall but are not associated with an underlying geometry. In this setting, there is no original \(X\) nor Euclidean coordinates. Nevertheless, we can use CMDS to represent these rankings as unit vectors in Euclidean space whose pairwise cosines approximate their rank correlations!

A real-world example is provided by the Quacquarelli Symonds World University Rankings, which include rankings of hundreds of world universities along six dimensions: academic reputation, employer reputation, faculty–student ratio, citations per faculty, international faculty ratio, and international student ratio. QS weight these rankings differently in their overall assessment, but our analysis will compare the rankings to each other across universities, ignoring these weights. The subset qswur_usa installed with ordr include U.S.-based universities ranked in the years 2017–2020, their classifications by QS, and their six integer-valued rankings. (Scores used to generate the rankings are included with the QS data files but omitted from qswur_usa.)

data(qswur_usa, package = "ordr")
print(qswur_usa)

## # A tibble: 612 × 13
##     year institution      size  focus res     age status rk_academic rk_employer
##    <int> <chr>            <fct> <fct> <fct> <int> <chr>        <int>       <int>
##  1  2017 MASSACHUSETTS I… M     CO    VH        5 B                6           4
##  2  2017 STANFORD UNIVER… L     FC    VH        5 A                5           5
##  3  2017 HARVARD UNIVERS… L     FC    VH        5 B                1           1
##  4  2017 CALIFORNIA INST… S     CO    VH        5 B               23          90
##  5  2017 UNIVERSITY OF C… L     FC    VH        5 B               13          47
##  6  2017 PRINCETON UNIVE… M     CO    VH        5 B               10          32
##  7  2017 YALE UNIVERSITY  L     FC    VH        5 B                9          12
##  8  2017 CORNELL UNIVERS… L     FC    VH        5 B               21          49
##  9  2017 JOHNS HOPKINS U… L     FC    VH        5 B               41         139
## 10  2017 UNIVERSITY OF P… L     FC    VH        5 B               33          35
## # ℹ 602 more rows
## # ℹ 4 more variables: rk_ratio <int>, rk_citations <int>,
## #   rk_intl_faculty <int>, rk_intl_students <int>

For this example, i’ll focus only on rankings for the year 2020 for which all rankings were available. This leaves me with only 38 universities, so my conclusions must be taken with caution! Since the rankings were subsetted from the full international data set, they are not contiguous (some integers between rankings never appear). To resolve this, i’ll also recalibrate the rankings by matching each vector of ranks to the vector of its sorted unique values:

qswur_usa %>%
  filter(year == 2020L) %>%
  select(institution, starts_with("rk_")) %>%
  mutate_at(
    vars(starts_with("rk_")),
    ~ match(., sort(unique(as.numeric(.))))
  ) %>%
  filter_at(vars(starts_with("rk_")), ~ ! is.na(.)) ->
  qswur_usa2020
print(qswur_usa2020)

## # A tibble: 38 × 7
##    institution     rk_academic rk_employer rk_ratio rk_citations rk_intl_faculty
##    <chr>                 <int>       <int>    <int>        <int>           <int>
##  1 MASSACHUSETTS …           3           2        8            4               1
##  2 STANFORD UNIVE…           2           3        5            7               2
##  3 HARVARD UNIVER…           1           1       17            5               8
##  4 CALIFORNIA INS…          12          17        1            2               3
##  5 UNIVERSITY OF …           9          11       19           22              20
##  6 PRINCETON UNIV…           7           7       35            1              26
##  7 CORNELL UNIVER…          11          13       40            9               5
##  8 UNIVERSITY OF …          14          12       11           31              12
##  9 YALE UNIVERSITY           6           4        2           46              17
## 10 COLUMBIA UNIVE…           8           8       10           47              51
## # ℹ 28 more rows
## # ℹ 1 more variable: rk_intl_students <int>

This subset of universities is now contiguously ranked along the six dimensions described above. The Kendall correlation \(\tau_{ij}\) between two rankings measures their concordance. To calculate it, every pair of universities contributes either \(+1\) or \(-1\) according as the rankings \(i\) and \(j\) place that pair in the same order, and the sum is scaled down by the number of pairs \({n\choose 2}\) so that the result lies between \(-1\) and \(1\). We interpret \(\tau_{ij}=1\) as perfect concordance (the rankings are equivalent), \(\tau_{ij}=-1\) as perfect discordance (the rankings are reversed), and \(\tau_{ij}=0\) as independence (the rankings are independent).

The six QS rankings are

• academic reputation (rk_academic),
• employer reputation (rk_employer)
• faculty–student ratio (rk_ratio)
• citations per faculty (rk_citations)
• international faculty ratio (rk_intl_faculty)
• international student ratio (rk_intl_students)

These not variations on a single theme (measures of the same construct), like different measures of guideline adherence or positive affect. Though they do all seem potentially sensitive to a university’s resources, including finance and prestige. I intuit that the two reputation indexes should be positively correlated, and that the two international personnel ratios should be as well. I also wonder if the faculty–student ratio might be anti-correlated with the number of citations per faculty, separating more research-focused institutions from more teaching-focused ones.

corr <- cor(select(qswur_usa2020, starts_with("rk_")), method = "kendall")
heatmap(corr, scale = "none")

eigen_ord(corr) %>%
  as_tbl_ord() %>%
  augment_ord() %>%
  mutate_rows(metric = stringr::str_remove(name, "rk_")) %>%
  confer_inertia(1) ->
  c_eigen
c_eigen %>%
  ggbiplot() +
  theme_minimal() +
  geom_unit_circle() +
  geom_rows_vector(aes(label = metric)) +
  scale_x_reverse(expand = expansion(add = .4)) +
  scale_y_continuous(expand = expansion(add = .3)) +
  ggtitle("Kendall correlations between university rankings",
          "CMDS correlation monoplot")

c_eigen %>%
  fortify(.matrix = "rows") %>%
  select(-name, -.matrix)

## # A tibble: 6 × 8
##      EV1     EV2     EV3     EV4     EV5      EV6 metric        .element
##    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>    <dbl> <chr>         <chr>   
## 1 -0.834 -0.0907  0.412   0.0430 -0.0206  0.351   academic      active  
## 2 -0.795 -0.0964  0.477  -0.0416 -0.181  -0.311   employer      active  
## 3 -0.517  0.771  -0.0480  0.331   0.158  -0.0372  ratio         active  
## 4 -0.731 -0.352  -0.239  -0.0278  0.528  -0.0685  citations     active  
## 5 -0.631 -0.233  -0.521   0.392  -0.352   0.00783 intl_faculty  active  
## 6 -0.603  0.262  -0.324  -0.665  -0.140   0.0312  intl_students active

c_eigen %>%
  ggbiplot(aes(x = 2, y = 3)) +
  theme_minimal() +
  geom_unit_circle() +
  geom_rows_vector(aes(label = metric)) +
  scale_x_continuous(expand = expansion(add = .5)) +
  scale_y_continuous(expand = expansion(add = .5)) +
  ggtitle("Kendall correlations between university rankings",
          "CMDS correlation monoplot, second & third principal coordinates")