---
title: "Metrics to Assess How Well Dimension Reduction Methods Preserve Data Set Properties"
author: "Kevin R. Coombes and Polina Bombina"
data: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Preservation Metrics}
  %\VignetteKeywords{DimensionReduction,Presevation,Pairwise Distances,Neighborhoods,Paths}
  %\VignetteDepends{Preservation}
  %\VignettePackage{Preservation}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r opts, echo=FALSE}
knitr::opts_chunk$set(fig.width=7, fig.height=5)
```

# Introduction
The introduction of high dimensional technologies in molecular biology, which perform
genome-scale measurements on DNA, RNA, or protein, has been accompanied by more extensive
use of "dimension reduction" (DR) algorithms that promote visualization of important aspects
of the data in just two dimensions. However, it is critical to understand the extent to
which different DR algorithms distort the properties of the data set that are present in
the original high dimensional space. The `Preservation` R package collects together some
metrics to aid our understanding.

# Preparation
## Packages
As a first step, we must load the `Preservation` library package.
```{r packs}
library(Preservation)
```

## Sample Data
We need some data in order to illustrate the preservation metrics. Here we use a subset of the 
CAVA CD4 T-cell single-cell CITE-Seq data set [1]. This data set was previously identified as
a single coherent cluster in the much larger complete data set. It contains measurements of the
expression of 51 proteins in 3096 individual cells. We also sub-clustered this subset into four 
parts to help with visualization.
```{r sampleData}
data("prot5")
dim(prot5)
table(colvec)
```

Next, we applied four popular dimension reduction methods:

1. Multi-Dimensional Scaling (MDS), which is equivalent to principal components analysis (OCA)
   when using Euclidean distance.
2. t-Stochastic Neighbor Embedding (t-SNE or TSNE).
3. Uniform Manifold Approximation and Projection (UMAP).
4. Diffusion Maps (DM).

The first two DR algorithms (MDS and TSNE) were computed using the `Rdimtools` package; 
UMAP, using the `umap` package; and DM, using the `diffusionMap` package.

```{r drs}
data("drs")
sapply(drs, dim)
```

Here is an initial plot showing the dimension reduction results.
```{r fig01, fig.width=8, fig.height=8, fig.cap="Four views of the CAVA subset after dimension reduction."}
opar <- par(mfrow = c(2,2))
for (tag in names(drs)) {
  plot(drs[[tag]], col = colvec, pch = 16, main = tag, xlab = "X1", ylab = "X2")
}
par(opar)
```

We first note that MDS and DM produce essentially the same plot, after recognizing that DM has
switched the (arbitrary) signs of both components. By contrast, both UMAP and TSNE have tried 
to place the purple group and half of the blue group in a separate cluster. The TSNE algorithm
is slightly more extreme, since it also seems to have "flipped" that new cluster so that the two
halves of the blue group are further apart.

# Preservation of Pairwise Distances
Some DR algorithms are designed to try to preserve the pairwise distances between points in the
data set. MDS is the most explicit example of such an algorithm, since its defined goal is to 
preserve the average of these distances as well as possible.

The `Preservation` package includes six metrics to assess how well pairwise distances are
preserved.  Here we apply five of these metrics to the four different DR views of the CAVA
dataset. (We omit Spearman's Rho, since it is somewhat slower than the others.)

```{r PWmetrics}
hiD <- dist(prot5)
loD <- lapply(drs, dist)
suppressWarnings(
  scores <- sapply(c("MilnorDistortion", "SigmaDistortion", "Stress", 
                     "EarthMover"), function(metric) {
    cat(metric, "\n", file = stderr())
    FUN <- get(metric)
    sapply(names(drs), function(tag) {
      FUN(hiD, loD[[tag]])
    })
  })
)
m1d <- sapply(names(drs), function(tag) {
  M1Distortion(prot5, drs[[tag]])
})
scores <- cbind(scores, M1Distortion = m1d)
round(scores, 3)
Preservation:::PWBest(colnames(scores))
```

DM is the best DR method based on the Milnor distortion; UMAP is best by Stress 
or Spearman's Rho (not shown); MDS is best based on the Sigma distortion, M1 
distortion, or Earth Mover's Distance. These findings probably tells us more
about how well the metrics assess their desired goal than about how well the DR
algorithms perform.

# Neighborhood Preservation
Method like TSNE and UMAP are inherently non-linear, and their heuristic is to try
to preserve local neighborhoods of points rather than all pairwise distances. The 
`Preservation` package includes several metrics to assess how well neighborhoods 
are preserved. These metrics are naturally functions of the number, $k$, of nearest
neighbors defining a neighborhood.

```{r neigh}
ajd <-sapply(names(drs), function(tag) {
  AverageJaccardDistance(prot5, drs[[tag]], k = 50)
})
scores <- cbind(scores, ajd)
# too slow: qnx <- QNX(hiD, loD[["MDS"]])
quad <- c(MDS = "#B00068", TSNE = "#1C8356", UMAP = "#325A9B", DM = "#FEAF16" )
ct <- lapply(names(drs), function(tag) {
  Preservation:::compute_ct_list(prot5, drs[[tag]])
})
avtrust <- sapply(ct, function(X) X[,3])
colnames(avtrust) <- names(drs)

avcont <- sapply(ct, function(X) X[,6])
colnames(avcont) <- names(drs)
```

```{r fig.width = 9, fig.height=5}
opar <- par(mfrow = c(1,2))
plot(avtrust[,1], type = "n", ylim = c(0.1, 0.3), xlab = "Neighbors", ylab = "Trustworthiness")
for (tag in names(drs)) lines(avtrust[, tag], lwd = 2, col = quad[tag])
legend("topright", names(drs), col = quad, lwd = 2)

plot(avcont[,1], type = "n", ylim = c(0, 0.25), xlab = "Neighbors", ylab = "Continuity")
for (tag in names(drs)) lines(avcont[, tag], lwd = 2, col = quad[tag])
legend("topleft", names(drs), col = quad, lwd = 2)
par(opar)
```

# Path Preservation
We now want to see how well DR algorithms preserve paths. We define a path to be an ordered set
of points in a data set. As an example, we computed a path in the original high (51-)dimensional
space that closely tracked the first principal component. So, in the original space, the path is 
a close approximation to a line segment that crosses the full data set. The path is represented
by 301 individual cells (indicated by their indices as rows in the data set).
```{r idex}
data("path-idex")
length(idex)
highD <- prot5[idex,]
dim(highD)
```

Here is a plot showing the dimension reductions with the overlaid image of the PC1 path.
```{r fig02, fig.width=8, fig.height=8, fig.cap="Four views of the path of the first principal components of the CAVA subset after dimension reduction."}
opar <- par(mfrow = c(2,2))
for (tag in names(drs)) {
  lowD <- drs[[tag]][idex,]
  plot(drs[[tag]], col = colvec, pch = 16, main = tag, xlab = "X1", ylab = "X2")
  lines(lowD, lwd = 2)
}
par(opar)
```

Again, the paths in the MDS and DM reductions jitter about the first component axis. In the
UMAP and TSNE reductions, the path seems to do a lot of jumping between the two pseudo-clusters
in the plot. In the TSNE plot, there are also numerous jumps between the two halves of the
blue group.

```{r LD}
pscores <- sapply(c("LengthDistortion", "SegmentVariance", "Curvature", "SpatialSimilarity"),
                 function(metric) {
  FUN <- get(metric)
  sapply(names(drs), function(tag) {
    lowD <- drs[[tag]][idex,]
    FUN(highD, lowD)
  })
})
round(pscores, 3)
Preservation:::PathBest(colnames(scores))
scores <- cbind(scores, pscores)
```

Each of the four metrics we have looked at compares the measurements on original and reduced
data sets by computing their log ratio (which is why the best score is the one nearest zero).
Here are variants that use Spearman's rank correlation instead, where te best score would be
the score closest to one. here length and curvature are topped by the UMAP and TSNE methods,
while spatial imilarity is better with the DM and MDS algorithms.

```{r rk}
rscores <- sapply(c("RankLengthDistortion", "RankCurvature", "RankSpatialSimilarity"),
                 function(metric) {
  FUN <- get(metric)
  sapply(names(drs), function(tag) {
    lowD <- drs[[tag]][idex,]
    FUN(highD, lowD)
  })
})
round(rscores, 3)
scores <- cbind(scores, rscores)
```

We find it surprising that UMAP and TSNE rank so highly under many of these metrics. Perhaps we
need to look more broadly to find metrics that give results that match our expectations?

One idea is to look at a smooth version of the path in the reduced spaces.
Here is a plot showing the dimension reductions with the overlaid image of smoothed PC1 path.
```{r fig03, fig.width=8, fig.height=8, fig.cap="Four views of the smoothed PC1 path of the CAVA subset after dimension reduction."}
smoosh <- function( coords) {
  ticks <- 1:nrow(coords) # arbitrary parametrization
  ## fit smooth curve to each coordinate as a function of the parameter
  fitted <- apply(coords, 2, function(X) predict(loess(X ~ ticks, span = 0.1)))
  delta <- fitted - coords # difference between actual path and smooth fitted curve
  val <- list(coords = coords, fitted = fitted, delta = delta)
  class(val) <- "smoosh"
  val
}
fits <- lapply(names(drs), function(tag) smoosh(drs[[tag]][idex,]))
names(fits) <- names(drs)
plot.smoosh <- function(obj, main = "", col = "purple", ...) {
  plot(obj$coords, main = main, type = "b", col = "grey80", ...)
  lines(obj$fitted, col = col, lwd = 2)
  invisible(main)
}
lines.smoosh <- function(obj, col = "purple") {
  lines(obj$fitted, col = col, lwd = 2)
}
opar <- par(mfrow = c(2, 2))
for (tag in names(fits)) {
  plot(fits[[tag]], main=tag, col = "black", xlab = "X1", ylab = "X2")
  points(drs[[tag]], col = colvec, pch = ".", cex=3)
  lines(fits[[tag]], col = "black")
}
par(opar)
```

```{r morePathhMetrics}
mscores <- sapply(c("Smoothness", "Coil", "EndpointDisplacement", "Contact"),
                 function(metric) {
  cat(metric, "\n", file = stderr())
  FUN <- get(metric)
  sapply(names(drs), function(tag) {
    lowD <- drs[[tag]][idex,]
    FUN(highD, lowD)
  })
})
tyscores <- sapply(c("max", "count"),
                 function(type) {
  cat(type, "\n", file = stderr())
  FUN <- Contact
  sapply(names(drs), function(tag) {
    lowD <- drs[[tag]][idex,]
    FUN(highD, lowD, type = type)
  })
})
cnscores <-  sapply(names(drs), function(tag) {
  lowD <- drs[[tag]][idex,]
  CrossingNumber(lowD)
})
mscores <- cbind(mscores, tyscores, cnscores)
round(mscores, 3)
scores <- cbind(scores, mscores)
```

```{r twist}
round(t(scores), 3)
```

# Appendix
This computation was performed using the following software tools.
```{r si}
sessionInfo()
```