• Basic Question: How to represent any data in two dimensions?

1. Projection
2. Squeeze data onto a table: close points stay close.

• Both ideas yield sensible visualization - but show a different aspect of the data.

• Represent high-dimensional point cloud in few (usually 2) dimensions keeping distances between points similar.
• Useful tool in visualizing any big data-set, specially for clustering purposes.
• Is a popular exploratory tool. Used before any inferential procedure.

• Given pairwise dissimilarities, reconstruct a map that preserves distances.
• From any dissimilarity (no need to be a metric)
• Reconstructed map has coordinates \(\mathbf{x}_i = (x_{i1}, x_{i2})\) and the natural distance \(\Vert x_i - x_j \Vert^2\).
• MDS is a family of different algorithms, each designed to arrive at optimal low-dimensional configuration (p = 2 or 3)
• Includes: Classical MDS, Metric MDS and Non-metric MDS.

• Problem: Given Euclidean Distance between points, recover the position of the points.
• Example: Road distance between 21 european cities

library(datasets); class(eurodist)

## [1] "dist"

                Athens Barcelona Brussels Calais Cherbourg Cologne
Barcelona         3313                                            
Brussels          2963      1318                                  
Calais            3175      1326      204                         
Cherbourg         3339      1294      583    460                  
Cologne           2762      1498      206    409       785  

eurocmd <- cmdscale(eurodist)
plot(eurocmd, type = "n")
text(eurocmd, rownames(eurocmd), cex = 0.7)

• Can identify points up to shift, reflection and rotation.

plot(eurocmd[,1], -eurocmd[,2], type = "n", asp = 1)
text(eurocmd[,1], -eurocmd[,2], rownames(eurocmd), cex = 0.7)

• Can identify points up to shift, reflection and rotation.

plot(x[,1], x[,2], type = "n")
text(x[,1], x[,2], labels = rownames(x), cex = 0.7)

• Classical MDS: Euclidean distance between \(n\) objects in \(p\) dimensions.
• Output: Position of points upto rotation, reflection, shift.
• Two steps:

1. Compute inner product matrix \(B = X X^T\) from distances.
2. Compute positions from \(B\).

• Distance, dissimilarity and similarity (or proximity) are defined for any pair of objects in any space.
• In mathematics, a distance function (that gives a distance between two objects) is also called metric, satisfying:

\[ d(x,y) \ge 0 \\ d(x,y) = 0 \text{ iff } x = y \\ d(x,y) = d(y,x) \\ d(x,z) \le d(x,y) + d(y,z) \]

• Given a dissimilarity matrix \(D = (d_{ij})\), MDS seeks to find \(x_1, \ldots, x_n \in \mathbb{R}^p\), such that: \[ d_{ij} \approx \Vert x_i - x_j \Vert^2 \text{ as close as possible} \]

• Oftentimes, for some large \(p\), there exists a configuration \(x_1, \ldots, x_n\) with exact distance match \(d_{ij} = \Vert x_i - x_j \Vert^2\).
• In such a case \(d\) is called a Euclidean distance.
• There are, however, cases where the dissimilarity is distance, but there exists no configuration in any \(p\) with perfect match.

• Such a distance is called non-Euclidean distance.

• Suppose for now we have Euclidean distance matrix \(D = (d_{ij})\)
• We want to find \((x_1, \ldots, x_n)\) such that \(d_{ij} = dist(x_i,x_j)\).
• Not unique solution: \(x^* = x + c\).
• Assume the observations are centered.

• Find inner product matrix \(B = X X^T\) instead of \(X\).
• Conenct to distance: \(d_{ij}^2 = b_{ii}+b_{jj}-2b_{ij}\).
• Center points to avoid shift invariance.
• Invert relationship: \[ b_{ij} = -1/2(d_{ij}^2 - d_{i.}^2 - d_{.j}^2 + d_{..}^2) \]
• "Doubly centered" distance.

• Since \(B = X X^T\), we need the square-root of \(B\).
• \(B\) is symmetric and positive definite.
• Can be diagonalised: \(B = V \Lambda V^T\).
• \(\Lambda\) is diagonal matrix with eigenvalues \(\lambda_1 > \lambda_2 > \cdots > \lambda_n\)
• Some eignevalues are zero, drop them. \(B = V_1 \Lambda_1 V_1^T\).
• Take "square-root" \(X = V_1 \Lambda_1^T\).

• Want a 2-D plot.
• Keep only largest 2 eignevectors and eignevalues.
• The resulting \(X\) will be the low-dimensional representation we were looking for.
• Same for 3-D.

library(randomForest)
ind <- sample(2,nrow(iris),replace=TRUE,prob=c(0.7,0.3))
trainData <- iris[ind==1,]
testData <- iris[ind==2,]
iris_rf <- randomForest(Species~.,data=trainData,ntree=100,proximity=TRUE)
table(predict(iris_rf),trainData$Species)

##             
##              setosa versicolor virginica
##   setosa         35          0         0
##   versicolor      0         37         2
##   virginica       0          1        34

plot(iris_rf)

varImpPlot(iris_rf)

irisPred<-predict(iris_rf,newdata=testData)
table(irisPred, testData$Species)

##             
## irisPred     setosa versicolor virginica
##   setosa         15          0         0
##   versicolor      0         11         2
##   virginica       0          1        12

• Proximity measure: The \((i,j)\) element of the proximity matrix produced by randomForest is the fraction of trees in which elements \(i\) and \(j\) fall in the same terminal node.
• The intuition is that similar observations should be in the same terminal nodes more often than dissimilar ones.
• We can use the proximity measure to identify structure.

set.seed(71)
iris.rf <- randomForest(Species ~ ., data=iris, importance=TRUE,
                        proximity=TRUE)
## Do MDS on 1 - proximity:
iris.mds <- cmdscale(1 - iris.rf$proximity, eig=TRUE)
## Look at variable importance:
round(importance(iris.rf), 2)

##              setosa versicolor virginica MeanDecreaseAccuracy
## Sepal.Length   5.88       5.87      9.21                10.62
## Sepal.Width    5.23       0.31      4.71                 4.94
## Petal.Length  21.60      31.41     27.71                32.39
## Petal.Width   22.96      33.74     32.07                33.85
##              MeanDecreaseGini
## Sepal.Length             9.37
## Sepal.Width              2.45
## Petal.Length            42.13
## Petal.Width             45.28

voting.mds <- cmdscale(voting, eig=TRUE)
names(voting.mds)

## [1] "points" "eig"    "x"      "ac"     "GOF"

plot(voting.mds$points[,1], voting.mds$points[,2],
     type = "n", xlab = "Coordinate 1", ylab = "Coordinate 2",
     xlim = range(voting.mds$points[,1])*1.2)
text(voting.mds$points[,1], voting.mds$points[,2], 
     labels = colnames(voting))

voting_mds <- isoMDS(voting)

## initial  value 15.268246 
## iter   5 value 10.264075
## final  value 9.879047 
## converged

plot(voting_mds$points[,1], voting_mds$points[,2],
     type = "n", xlab = "Coordinate 1", ylab = "Coordinate 2",
     xlim = range(voting_mds$points[,1])*1.2)
text(voting_mds$points[,1], voting_mds$points[,2], 
     labels = colnames(voting))