Proximity Scaling
For proximity scaling (PS) or multidimensional scaling (MDS) the input is typically an \(N\times N\) matrix \(\Delta^*=f(\Delta)\), a matrix of proximities with elements \(\delta^*_{ij}\), that is a function of a matrix of observed non-negative dissimilarities \(\Delta\) with elements \(\delta_{ij}\). \(\Delta^*\) usually is symmetric (but does not need to be). The main diagonal of \(\Delta\) is 0. We call a \(f: \delta_{ij} \mapsto \delta^*_{ij}\) a proximity transformation function. In the MDS literature these \(\delta_{ij}^*\) are often called dhats or disparities. The problem that proximity scaling solves is to locate an \(N \times M\) matrix \(X\) (the configuration) with row vectors \(x_i, i=1,\ldots,N\) in low-dimensional space \((\mathbb{R}^M, M \leq N)\) in such a way that transformations \(g(d_{ij}(X))\) of the fitted distances \(d_{ij}(X)=d(x_i,x_j)\)—i.e., the distance between different \(x_i, x_j\)—approximate the \(\delta^*_{ij}\) as closely as possible. We call a \(g: d_{ij}(X) \mapsto d_{ij}^*(X)\) a distance transformation function. In other words, proximity scaling means finding \(X\) so that \(d^*_{ij}(X)=g(d_{ij}(X))\approx\delta^*_{ij}=f(\delta_{ij})\). There may also be weights transformed weight \(w^*_{ij}=h(w_{ij})\) with \(h: w_{ij} \mapsto w_{ij}^*\) being a weight transformation function.
This approximation \(D^*(X)\) to the matrix \(\Delta^*\) is found by defining a fit criterion (the loss function), \(\sigma_{PS}(X)=L(\Delta^*,D^*(X))\), that is used to measure how closely \(D^*(X)\) approximates \(\Delta^*\).
The loss function used is then minimized to find the vectors \(x_1,\dots,x_N\), i.e., \[\begin{equation} \label{eq:optim} \arg \min_{X}\ \sigma_{PS}(X). \end{equation}\]
Stress Models
The first popular type of PS supported in stops is based
on the loss function type stress (Kruskall 1964), employing a
quadratic loss function.
A general formulation of a loss function based on a quadratic loss is \[\begin{equation} \label{eq:stress} \sigma_{PS}(X)=\sigma_{stress}(X)=\sum^N_{i=1}\sum^N_{j=1} w^*_{ij}\left[d^*_{ij}(X)-\delta^*_{ij}\right]^2=\sum^N_{i=1}\sum^N_{j=1} h(w_{ij})\left[g\left(d_{ij}(X)\right)-f(\delta_{ij})\right]^2 \end{equation}\] Here, the \(w_{ij}\) are finite weights, with \(w_{ij}=0\) if the entry is missing and \(w_{ij}=1\) otherwise. There are a number of optimization techniques one can use to solve this optimization problem.
The distances used is usually some type of Minkowski distance (\(p > 0\)) as the distance fitted to the points in the configuration, \[\begin{equation} \label{eq:dist} d_{ij}(X) = ||x_{i}-x_{j}||_p=\left( \sum_{m=1}^M |x_{im}-x_{jm}|^p \right)^{1/p} \ i,j = 1, \dots, N. \end{equation}\] Typically, the norm used is the Euclidean norm, so \(p=2\). In standard MDS \(g(\cdot)=f(\cdot)=I(\cdot)\), the identity function.
This formulation enables one to express a large number of PS methods
many of which are implemented in stops. In
stops we allow to use specific choices for \(f(\cdot)\), \(g(\cdot)\) and \(h(\cdot)\) from the family of power
transformations so one can fit the following stress models:
- Explicitly normalized stress: \(w^*_{ij}=(\sum_{ij}\delta^{*2}_{ij})^{-1}\), \(\delta_{ij}^*=\delta_{ij}\), \(d^*_{ij}(X)=d_{ij}(X)\)
- Stress-1: \(w^*_{ij}=(\sum_{ij} d^{*2}_{ij}(X))^{-1}\), \(\delta_{ij}^*=\delta_{ij}\), \(d_{ij}(X)^*=d_{ij}(X)\)
- Sammon stress (Sammon 1969): \(w^*_{ij}=\delta^{*-1}_{ij}\) , \(\delta_{ij}^*=\delta_{ij}\), \(d_{ij}(X)^*=d_{ij}(X)\)
- Elastic scaling stress (McGee 1966): \(w_{ij}=\delta^{*-2}_{ij}\), \(\delta_{ij}^*=\delta_{ij}\), \(d_{ij}(X)^*=d_{ij}(X)\)
- S-stress (Takane et al. 1977): \(\delta^*_{ij}=\delta_{ij}^2\) and \(d^*_{ij}(X)=d^2_{ij}(X)\), \(w^*_{ij}=1\)
- R-stress (de Leeuw, 2014): \(\delta^*_{ij}=\delta_{ij}\) and \(d^*_{ij}=d^{2r}_{ij}\), \(w^*_{ij}=1\)
- Power MDS (Buja et al. 2008, Rusch, Mair & Hornik 2021): \(\delta^*_{ij}=\delta_{ij}^\lambda\) and \(d^*_{ij}=d^\kappa_{ij}\), \(w^*_{ij}=1\)
- Power elastic scaling (Rusch, Mair & Hornik 2023): \(w^*_{ij}=\delta^{*-2}_{ij}\), \(\delta^*_{ij}=\delta_{ij}^\lambda\) and \(d^*_{ij}=d^\kappa_{ij}\)
- Power Sammon mapping (Rusch, Mair & Hornik 2023): \(w^*_{ij}=\delta^{*-1}_{ij}\), \(\delta^*_{ij}=\delta_{ij}^\lambda\) and \(d^*_{ij}=d^\kappa_{ij}\)
- Restricted powerstress (Rusch, Mair & Hornik, 2021): \(\delta^*_{ij}=\delta_{ij}^\lambda\), \(d^*_{ij}=d^\kappa_{ij}\) but with the restriction \(\kappa=\lambda\). The \(w^*_{ij}=w_{ij}^\rho\) for arbitrary \(w_{ij}\) (e.g., a function of the \(\delta_{ij}\))
- Power stress (encompassing all previous models; Buja et al. 2008, Rusch, Mair & Hornik 2021): \(\delta^*_{ij}=\delta_{ij}^\lambda\), \(d^*_{ij}=d^\kappa_{ij}\) and \(w^*_{ij}=w_{ij}^\rho\) for arbitrary \(w_{ij}\) (e.g., a function of the \(\delta_{ij}\))
- Approximate power stress (Rusch, Mair & Hornik, 2021): This approximates power stress by \[\begin{equation} \label{eq:apstress} \sigma_{PS}(X)=\sigma_{apstress}(X)=\sum_{i<j} \delta_{ij}^\upsilon \left(\delta_{ij}^\tau - d_{ij}(X)\right)^2 \end{equation}\] where the relation to p-stress is so that \(\upsilon=\rho + 2\lambda(1 -1/\kappa)\) and \(\tau=\lambda/\kappa\). The approximation of p-stress by ap-stress works well in cases when for the \(x_i,x_j\) for which \(w_{ij}\) is large, the error \(\delta_{ij}^\lambda - d_{ij}(X)^\kappa\) is small, so that \(d_{ij}(X)^\kappa\) is approximated reasonably well by \(\delta_{ij}^\lambda\) and, equivalently, \(d_{ij}(X)^\kappa\) can be approximated well by \(d_{ij}(X) \delta_{ij}^{(\lambda (\kappa-1)/\kappa)}\).
- Local MDS (LMDS; Chen & Buja, 2009): Let \(N_k\) define the symmetric set of nearby pairs of points \((i,j)\) so that \((i,j)\in N_k\) iff \(i\) is among the \(k-\)nearest neighbours of \(j\) or the other way round. Then we have a stress model with \[\begin{equation} \delta^*_{ij}=\begin{cases} \delta_{ij} &\mbox{ if } (i,j) \in N_k\\ \sqrt{w}\delta_{\infty} &\mbox{ if } (i,j) \in N_k\\ \end{cases} \end{equation}\] and \[\begin{equation} d^*_{ij}(X)=\begin{cases} d_{ij}(X) &\mbox{ if } (i,j) \in N_k\\ \sqrt{w}d_{ij}(X) &\mbox{ if } (i,j) \in N_k\\ \end{cases} \end{equation}\] where \(\delta_{\infty} \rightarrow \infty\) is a large ``imputed’’ dissimilarity that is constant and \(w\) a small weight. For simplification, one takes \(w \sim 1/\delta_{\infty}\) in the standard lMDS objective, so one uses the tuning parameter \(\tau=2w\delta_{\infty}\) for given \(k\).
For all of the above models but local MDS one can use the function
powerStressMin which uses majorization to find the solution
(De Leeuw, 2014) . The function allows to specify a kappa,
lambda and nu argument as well as a
weightmat (the \(w_{ij}\)). LMDS can be fitted with the
lmds function. For approximate power stress there also is a
convenience function apStressMin.
Note that the implementation in powerStressMin is more a
proof-of-concept than optimal. Majorizing this type of stress is a
pretty hard problem and the code we use relies on a while loop in pure
R. We plan to speed the loop up with a re-implementation in C in the
future.
Strain Models
The second popular type of PS supported in stops is
based on the loss function type . Here the \(\Delta^*\) are a transformation of the
\(\Delta\), \(\Delta^*= f (\Delta)\) so that \(f(\cdot)=-(h\circ l)(\cdot)\) where \(l\) is any function and \(h(\cdot)\) is a double centering operation,
\(h(\Delta)=\Delta-\Delta_{i.}-\Delta_{.j}+\Delta_{..}\)
where \(\Delta_{i.}, \Delta_{.j},
\Delta_{..}\) are matrices consisting of the row, column and
grand marginal means respectively. These then get approximated by
(functions of) the inner product matrices of \(X\) \[\begin{equation}
\label{eq:dist2}
d_{ij}(X) = \langle x_{i},x_{j} \rangle
\end{equation}\] We can thus express classical scaling as a
special case of the general PS loss with \(d_{ij}(X)\) as an inner product, \(g(\cdot) = I(\cdot)\) and \(f(\cdot)=-(h \circ I)(\cdot)\).
If we again allow power transformations for \(g(\cdot)\) and \(f(\cdot)\) one can fit the following strain
models with stops
- Classical scaling (Torgerson, 1958): \(\delta^*_{ij}=-h(\delta_{ij})\) and \(d^*_{ij}=d_{ij}\)
- Powerstrain (Buja et al. 2008, Rusch, Mair & Hornik 2021): \(\delta^*_{ij}=-h(\delta_{ij}^\lambda)\), \(d^*_{ij}=d_{ij}\) and \(w_{ij}=w_{ij}^\rho\) for arbitrary \(w_{ij}\)
In stops we have a wrapper to cmdscale
(overloading the base function) which extend functionality
by offering an object that matches smacofP objects with
corresponding methods.
Energy Models
The third type of PS supported in stops is based on the
energy type losses of pairwise attraction and repulsion between objects
(Chen & Buja, 2014, Rusch, Mair, Hornik 2023a). It is related to
graph drawing: if vertices are equated with objects, the edge weights
are the \(\delta^*_{ij}\), the node
positions in the layout is the \(X\)
and the distance between nodes in the layout are \(d_{ij}(X)\).
One can express badness-of-fit as comprising a repulsion part \(\propto -\delta^*_{ij}d^*_{ij}(X)\) and an
attraction part \(\propto
d^{*\mu}_{ij}(X)\) (Chen & Buja, 2014). This is
\[\begin{equation}
\label{eq:energy}
\sigma_{PS}(X)=\sigma_{EM}(X) = \sum_{i<j} w^*_{ij} \left(a
d^{*\mu}_{ij}(X) - b \delta^*_{ij}d^*_{ij}(X) \right)
\end{equation}\] with \(a,b\)
being constants and \(d^*_{ij}(X)=g(d_{ij}(X))\) and \(\delta^*_{ij}=f(d_{ij}(X))\). For \(\mu=2, a=1\) and \(b=2\) this is a stress model where terms
depending solely on \(\delta^*_{ij}\)
are disregarded for finding \(X\).
Specific instances of energy models that can be fitted in
stops are:
- Box-Cox MDS (Chen & Buja, 2014). This is a three-parameter family \(\theta=(\mu^{(BC)},\lambda^{(BC)},\rho^{(BC)})^\top\): \[\begin{equation} \label{eq:bcstress} \sigma_{PS}(X)=\sigma_{BCstress}(X) = \sum_{i<j} \delta_{ij}^{\rho} \left( BC_{\mu+\lambda}(d_{ij}(X)) - \delta_{ij}^{\lambda} BC_{\mu}(d_{ij}(X)) \right) \end{equation}\] with \(\mu, \rho\in \mathbb{R}\) and \(\lambda \in \mathbb{R}_+\). Here \(BC_\alpha\) is the one-parameter Box-Cox transformation with parameter \(\alpha\), \[\begin{equation} \label{eq:bc} BC_{\alpha}(d)=\begin{cases} \frac{d^{\alpha}-1}{\alpha} \mbox{ if } \alpha \neq 0\\ \log(d) \mbox{ if } \alpha = 0\\ \end{cases} \end{equation}\] For \(w_{ij}=\delta_{ij}, \mu=\lambda=1\) this is very similar to the stress-equivalent energy MDS model. Note that here the distance transformation functions in the attraction and repulsion parts need not be equal. \(\sigma_{\text{BC}}(X)\) shares close similarities with stress with power transformations (power stress): They will often lead to a similar configurations for the same parameters (even equivalent if the exponents for repulsion and attraction parts are the same).
Usage of flexible MDS models
The objects returned from powerStressMin,
lmds, bcStressMin and apStressMin
inherit from class smacofP which extends the
smacof classes (De Leeuw & Mair, 2009) to allow for the
transformations. Apart from that all the objects are made so that they
are compatible to methods from smacof and also inherit from
smacof. Accordingly, the following S3 methods are
available:
| Method | Description |
|---|---|
| Prints the object | |
| summary | A summary of the object |
| plot | 2D Plots of the object |
| plot3d | Dynamic 3D configuration plot |
| plot3dstatic | Static 3D configuration plot |
| residuals | Residuals |
| coef | Model Coefficients |
Let us illustrate the usage.
First we setup our dissimilarity matrix \([\delta_{ij}]\) and call it
dis.
dis<-as.matrix(smacof::kinshipdelta)We now fit a series of MDS models available in
stops.
- A standard MDS (
stress)
res1<-powerStressMin(dis,kappa=1,lambda=1)
res1##
## Call:
## powerStressMin(delta = dis, kappa = 1, lambda = 1)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.264
## Number of iterations: 5352- A
sammonmapping
res2<-powerStressMin(dis,kappa=1,lambda=1,nu=-1,weightmat=dis)
res2##
## Call:
## powerStressMin(delta = dis, kappa = 1, lambda = 1, nu = -1, weightmat = dis)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.289
## Number of iterations: 81666Alternatively, one can use the faster sammon function
from MASS (Venables & Ripley, 2002) for which we
provide a wrapper that adds class attributes and methods (and overloads
the function).
res2a<-sammon(dis)## Initial stress : 0.17053
## stress after 3 iters: 0.10649res2a##
## Call: sammon(d = dis)
##
## Model: Sammon Scaling
## Number of objects: 15
## Stress: 0.1064925- An
elasticscaling
res3<-powerStressMin(dis,kappa=1,lambda=1,nu=-2,weightmat=dis)
res3##
## Call:
## powerStressMin(delta = dis, kappa = 1, lambda = 1, nu = -2, weightmat = dis)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.305
## Number of iterations: 1e+05- A
sstressmodel
res4<-powerStressMin(dis,kappa=2,lambda=2)
res4##
## Call:
## powerStressMin(delta = dis, kappa = 2, lambda = 2)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.346
## Number of iterations: 47130- An
rstressmodel (with \(r=1\) as \(r=\kappa/2\))
res5<-powerStressMin(dis,kappa=2,lambda=1)
res5##
## Call:
## powerStressMin(delta = dis, kappa = 2, lambda = 1)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.404
## Number of iterations: 21201- A
powermdsmodel (i.e., \(\rho=1\))
res6<-powerStressMin(dis,kappa=2,lambda=1.5)
res6##
## Call:
## powerStressMin(delta = dis, kappa = 2, lambda = 1.5)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.367
## Number of iterations: 50564- A
powersammonmodel
res7<-powerStressMin(dis,kappa=2,lambda=1.5,nu=-1,weightmat=dis)
res7##
## Call:
## powerStressMin(delta = dis, kappa = 2, lambda = 1.5, nu = -1,
## weightmat = dis)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.436
## Number of iterations: 1e+05- A
powerelasticscaling
res8<-powerStressMin(dis,kappa=2,lambda=1.5,nu=-2,weightmat=dis)
res8##
## Call:
## powerStressMin(delta = dis, kappa = 2, lambda = 1.5, nu = -2,
## weightmat = dis)
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.519
## Number of iterations: 1e+05- A
powerstressmodel
res9<-powerStressMin(dis,kappa=2,lambda=1.5,nu=-1.5,weightmat=2*1-diag(nrow(dis)))
res9##
## Call:
## powerStressMin(delta = dis, kappa = 2, lambda = 1.5, nu = -1.5,
## weightmat = 2 * 1 - diag(nrow(dis)))
##
## Model: Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.367
## Number of iterations: 57237summary(res9)##
## Configurations:
## D1 D2
## Aunt -0.1225 0.2498
## Brother 0.1964 -0.1400
## Cousin 0.0525 0.3099
## Daughter -0.2050 -0.1256
## Father 0.1639 -0.1822
## Granddaughter -0.2358 -0.0531
## Grandfather 0.2146 -0.1336
## Grandmother -0.2360 -0.0868
## Grandson 0.2147 -0.1079
## Mother -0.2098 -0.1363
## Nephew 0.1707 0.2110
## Niece -0.1231 0.2442
## Sister -0.2219 -0.0967
## Son 0.1702 -0.1707
## Uncle 0.1710 0.2181
##
##
## Stress per point:
## SPP SPP(%)
## Niece 0.0008 4.6861
## Nephew 0.0008 4.7097
## Aunt 0.0008 4.9367
## Uncle 0.0008 5.0269
## Daughter 0.0010 6.0875
## Son 0.0010 6.2785
## Father 0.0011 6.4569
## Mother 0.0011 6.6921
## Cousin 0.0011 6.8096
## Sister 0.0012 7.2765
## Brother 0.0012 7.2882
## Grandson 0.0013 7.9701
## Granddaughter 0.0013 8.0950
## Grandmother 0.0015 8.8103
## Grandfather 0.0015 8.8761- An
lmdsmodel
res10<-lmds(dis,tau=0.5,k=3)
res10##
## Call:
## lmds(delta = dis, k = 3, tau = 0.5)
##
## Model: Local MDS
## Number of objects: 15
## Stress-1 value: 0.302
## Number of iterations: 80- An
apstressmodel (approximate power stress)
res11<-apStressMin(dis,tau=0.5,ups=2)
res11##
## Call:
## apStressMin(delta = dis, tau = 0.5, ups = 2)
##
## Model: Approximate Power Stress SMACOF
## Number of objects: 15
## Stress-1 value: 0.217
## Number of iterations: 16We can visualize the obtained configurations with plot. For
smacofP objects we have the following plots.
plot(res9)
plot(res9,"transplot")
plot(res9,"Shepard")
plot(res9,"resplot")
plot(res9,"bubbleplot")
powerstrainmodel is rather easy to fit with simply subjecting the dissimilarity matrix to some power. Here we use \(\lambda=3\).
resc<-cmdscale(kinshipdelta^3)
resc##
## Call: cmdscale(d = kinshipdelta^3)
##
## Model: Torgerson-Gower Scaling
## Number of objects: 15
## GOF: 0.4257747 0.6281985summary(resc)##
## Configurations:
## D1 D2
## Aunt 178193.10 204986.70
## Brother -174369.39 -94357.29
## Cousin -48355.46 265057.17
## Daughter 169149.17 -109936.27
## Father -145389.13 -168604.23
## Granddaughter 187039.12 -44850.80
## Grandfather -180116.08 -103668.07
## Grandmother 199145.19 -83038.71
## Grandson -169767.75 -72358.82
## Mother 185797.71 -138677.31
## Nephew -195963.98 173123.81
## Niece 168278.13 208870.31
## Sister 182224.16 -70764.01
## Son -149876.11 -135883.18
## Uncle -205988.70 170100.68summary(resc)
plot(resc)- A BC MDS with
bcStressMin
resbc<-bcStressMin(dis,mu=3,lambda=2,nu=2)
resbc##
## Call:
## bcStressMin(delta = dis, mu = 3, lambda = 2, nu = 2)
##
## Model: Box-Cox Stress MDS
## Number of objects: 15
## Stress-1 value: 0.251
## Number of iterations: 67plot(resbc)
