```{r vegetables_plot, fig.align = "center", echo = FALSE} plot(1:9, h0$x, type = "n", xlab = "vegetable", ylab = "scale") text(1:9, h0$x, nms, col = "BLUE") lines(1:9, h1$x) lines(1:9, h2$x) lines(1:9, h3$x) lines(1:9, h4$x) lines(1:9, h5$x) lines(1:9, h0$x, col = "RED", lwd = 2) ```

##Nested MDS In principal component analysis and classical MDS solutions for different dimensionalities are *nested*, in the sense that the $p$-dimensional solution defines the first $p$ dimensions of the $p+1$-dimensional solution. One possible way around this problem is to define the $p$-dimensional solution as the first $p$ principal components of the full-dimensional MDS solution (@deleeuw_groenen_mair_E_16e). But we can also use additive components to go from the $p$-dimensional to the $p+1$-dimensional solution. Simply set the $f_{ij}$ equal to the distances of the $p$-dimensional configuration, and use $$ B(X^{(k)})=\sum_{i=1}^n\sum_{j=1}^nw_{ij}\ \xi_{ij}(X_{ij}^{(k)})A_{ij}, $$ with $$ \xi_{ij}(X_{ij}^{(k)})=\frac{\delta_{ij}}{\sqrt{(x_i^{(k)}-x_j^{(k)})^2+f_{ij}^2}}. $$ to compute the next dimension. As an illustration we use data from @degruijter_67, with average dissimilarity judgments between Dutch political parties in 1967. ```{r poldist_data_twice, echo = FALSE} poldist <- structure(c(5.63, 5.27, 4.6, 4.8, 7.54, 6.73, 7.18, 6.17, 6.72, 5.64, 6.22, 5.12, 4.59, 7.22, 5.47, 5.46, 4.97, 8.13, 7.55, 6.9, 4.67, 3.2, 7.84, 6.73, 7.28, 6.13, 7.8, 7.08, 6.96, 6.04, 4.08, 6.34, 7.42, 6.88, 6.36, 7.36), Labels = c("KVP", "PvdA", "VVD", "ARP", "CHU", "CPN", "PSP", "BP", "D66"), Size = 9L, call = quote(as.dist.default(m = polpar)), class = "dist", Diag = FALSE, Upper = FALSE) print(poldist) pmedian <- median (poldist) ``` ```{r poldist_run, echo = FALSE} delta <- as.matrix (poldist) nms <- rownames (delta) h <- asmacof (delta, f = 0, p = 1) f <- sqrt(.5) * as.matrix(dist (h$x)) g <- asmacof (delta, f, p = 1) r <- asmacof (delta, f = 0) xsuc <- cbind(h$x, g$x) xsim <- r$x ssuc <- sum ((delta - as.matrix (dist (xsuc))) ^ 2) ssim <- sum ((delta - as.matrix (dist (xsim))) ^ 2) ``` The two-dimensional successive solution is in figure `r figure_nums("gruijter_suc", display = "n")`. It has stress `r ssuc`. The simultaneous two-dimensional solution is in figure `r figure_nums("gruijter_sim", display = "n")`, with stress `r ssim`.

```{r poldist_plot_suc, echo = FALSE, fig.align = "center"} plot(xsuc, type = "n", xlab = "dim 1", ylab = "dim 2") text (xsuc, nms, col = "RED") ```

```{r poldist_plot_sim, echo = FALSE, fig.align = "center"} plot(xsim, type = "n", xlab = "dim 1", ylab = "dim 2") text (xsim, nms, col = "RED") ```

##Additive Constant In regularized MDS $\epsilon$ is a known constant. We can also design a version in which $\epsilon$ is estimated. The iteration for updating $\epsilon$ turns out to be $$ \epsilon^{(k+1)}=\epsilon^{(k)}\frac{\sum_{i=1}^n\sum_{j=1}^nw_{ij}\xi_{ij}(X_{ij}^{(k)})}{\sum_{i=1}^n\sum_{j=1}^nw_{ij}}, $$ with $$ \xi_{ij}(X_{ij}^{(k)})=\frac{\delta_{ij}}{\sqrt{d_{ij}^2(X^{(k)})+\{\epsilon^{(k)}\}^2}}. $$ The De Gruijter data are suitable for illustrating the additive constant algorithm, because they are averaged dissimilarity judgments over a large number of students, and thus they are likely to exhibit regression to the mean. ```{r poldist_additive_constant_run, echo = FALSE, cache = TRUE} h <- bsmacof (delta) rownames (h$x) <- nms ``` The configuration ```{r poldist_additive_constant_cconf, echo = FALSE} h$x ``` is given in figure `r figure_nums("gruijter_add", display = "n")`. It requires `r h$itel` iterations and has stress `r h$stress`. Allowing estimation of an additive constant does not change the configuration very much, but stress is much lower than before. We do see the three Christian Democrat parties KVP, ARP, and CHU moving much closer together.

```{r poldist_additive_constant_plot, fig.align = "center", echo = FALSE} plot(h$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h$x, nms, col = "RED") ```

##MDS with Restrictions Our examples so far can also be analyzed using the framework of @deleeuw_heiser_C_80. If all $X_{ij}$ are the same $$ \rho(X,X^{(k)})=\mathbf{tr}\ X'VX-2\mathbf{tr}\ X'B(X^{(k)})X^{(k)}+\mathbf{tr}\ (X^{(k)})'H(X^{(k)})X^{(k)}, $$ with $$ H(X^{(k)}):=\sum_{i=1}^n\sum_{j=1}^nw_{ij}\xi_{ij}^2(X^{(k)})A_{ij}. $$ Thus a `smacof` step minimizes $$ \mathbf{tr}\ (X-V^+B(X^{(k)})X^{(k)})'V(X-V^+B(X^{(k)})X^{(k)}). $$ Thus if there are restrictions on the possible configurations, say $X\in\Omega$, then we first compute the Guttman transform and then its (weighted) least squares projection on the constraint set $\Omega$. For all problems in which the common configuration is partitioned into two parts $\begin{bmatrix}X_1&\mid&X_2\end{bmatrix}$ we have separation of the loss components \begin{align*} \rho(X,X^{(k)})&=\mathbf{tr}\ (X_1-V^+B(X^{(k)})X^{(k)}_1)'V(X_1-V^+B(X^{(k)})X^{(k)}_1)\\&+\mathbf{tr}\ (X_2-V^+B(X^{(k)})X^{(k)}_2)'V(X_2-V^+B(X^{(k)})X^{(k)}_2). \end{align*} For the regularized MDS and additive constant problems $X$ is of the form $\begin{bmatrix}X_1&\mid&\epsilon I\end{bmatrix}$, for nested MDS we have $\begin{bmatrix}X_1&\mid&X_2\end{bmatrix}$, with $X_1$ the fixed lower-dimensional solution. Thus for the regularized MDS problem and for the additive constant problem $X^{(k+1)}_1=V^+B(X^{(k)})X^{(k)}_1$, and for the additiive constant problem $$ \epsilon^{(k+1)}=\frac{\mathbf{tr}\ VV^+B(X^{(k)})X^{(k)}_2}{\mathbf{tr}\ V}=\epsilon^{(k)}\frac{\mathbf{tr}\ B(X^{(k)})}{\mathbf{tr}\ V}. $$ For nested MDS simply $X_2^{(k+1)}=V^+B(X^{(k)})X^{(k)}_2$. #Partitioned Dissimilarities ##Piecewise MDS Consider the following situation. The dissimilarities are partitioned into, say, $m$ groups, and the configurations $X_\ell$ are constrained as $X_\ell=X\Lambda_\ell$ with $\Lambda_\ell$ diagonal. The groups of dissimilarities are arbitrary. Perhaps dissimilarities come from different sources, or have different degrees of precision. We shall treat the general case in which both $X$ and the $\Lambda_\ell$ must be estimated. The case in which the $\Lambda_\ell$ are fixed and known is included as a special case in the general treatment. In iteration $k$ we must minimize \begin{equation}\label{E:ploss} \sum_{\ell=1}^m\sum_{(i,j)\in\mathcal{I}_\ell}w_{ij}\ \mathbf{tr}\ (X\Lambda_\ell-\xi_{ij}(X_\ell^{(k)})X_\ell^{k})'A_{ij}(X\Lambda_\ell-\xi_{ij}(X_\ell^{(k)})X_\ell^{k}), \end{equation} which becomes $$ \sum_{\ell=1}^m\mathbf{tr}\ \Lambda_\ell X'V_\ell X\Lambda_\ell-2\sum_{\ell=1}^m\mathbf{tr}\ \Lambda_\ell X'B_\ell(X_\ell^{(k)})X_\ell^{(k)}+\sum_{\ell=1}^m\mathbf{tr}\ (X_l^{(k)})'H_\ell X_l^{(k)}. $$ The stationary equations are \begin{equation}\label{E:stat1} \sum_{\ell=1}^m V_\ell X\Lambda_\ell^2=\sum_{\ell=1}^mB_\ell(X_\ell^{(k)})X_\ell^{(k)}\Lambda_\ell \end{equation} and, in the case we are also minimizing over $\Lambda_\ell$, \begin{equation}\label{E:stat2} \Lambda_\ell=\frac{\mathbf{diag}(X'B_\ell(X_\ell^{(k)})X_\ell^{(k)})}{\mathbf{diag}(X'V_\ell X)} \end{equation} There is little hope to solve these two equations analytically, so we resort to alternating least squares. It is easy enough to solve $\eqref{E:stat2}$ for the $\Lambda_\ell$ for given $X$. To solve $\eqref{E:stat1}$ for $X$ for given $\Lambda_\ell$ consider column $s$ from both sides of the equation. This gives $$ x_s=\left[\sum_{\ell=1}^m\Lambda_{ls}^2V_\ell^{\ }\right]^+u_s, $$ where $u_s$ is column $s$ of the matrix on the right-hand side of $\eqref{E:stat1}$. Note that if at some point a diagonal element of $\Lambda_\ell$ is zero, it will remain zero in subsequent iterations. Thus we can impose zeroes in the $\Lambda_\ell$ by fiLambdang them in the initial estimate. The De Gruijter data show the familar bending or horseshoe effect, where the left right scale is bended to bring the extreme left and right closer than they would otherwise be. We first analyze the data with the `psmacof` function, using a single common dimension and with another single dimension for the largest dissimilarities. The cutoff is the median of the $\delta_{ij}$, which is `r pmedian`. Thus there are two groups: for the smaller dissimilarities below the median both only the first dimension gets unit weight, for the larger dissimilarities both dimensions have unit weight. The $\Lambda_\ell$ are fixed, and there is no need to minimize over them. ```{r piecewise_constant, echo = FALSE, cache = TRUE} pattern <- array (0, c(9, 9, 2)) pattern[, , 1] <- ifelse (delta <= pmedian, 1, 0) pattern[, , 2] <- ifelse (delta > pmedian, 1, 0) lambda <- matrix (c(1, 1, 0, 1), 2, 2) h<-psmacof(delta, pattern=pattern, p=2, lambda=lambda, verbose=FALSE, do_lambda=FALSE, itmax=10000) rownames(h$x) <- nms ``` After `r h$itel` iterations we convergence on the configuration ```{r piecewise_constant_conf, echo = FALSE} h$x ``` which is plotted in figure `r figure_nums("pcc", display = "n")` and has stress `r h$stress`.

```{r piecewise_constant_plot, echo = FALSE, fig.align = "center"} plot(h$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h$x, nms, col = "RED") ```

In the next `psmacof` analysis, of the same data, we have two common dimensions and one extra dimension for the large dissimilarities. We also allow optimization over the $\Lambda_\ell$. ```{r piecewise_variable, echo = FALSE, cache = TRUE} lambda <- matrix (c(1, 1, 1, 1, 0, 1), 2, 3) h<-psmacof(delta,pattern=pattern,p=3,lambda=lambda,verbose=FALSE,do_lambda=TRUE,itmax=10000) rownames(h$x) <- nms ``` The configuration is ```{r piecewise_variable_conf, echo = FALSE} h$x ``` and we plot it in figure `r figure_nums("pvc", display = "n")`. It requires `r h$itel` iterations and has stress `r h$stress`.

```{r piecewise_variable_fig, echo = FALSE, fig.align = "center"} plot(h$x[,1:2], type = "n", xlab = "dim 1", ylab = "dim 2") text (h$x[,1:2], nms, col = "RED") ```

The algorithm becomes only slightly more complicated if we do not require the $\Lambda_\ell$ to be diagonal, but allow them to be full matrices. ##INDSCAL/PARAFAC An interesting special case of $\eqref{E:ploss}$, where the dissimilarity groups are defined by different individuals or occasions or time-points, is $$ \sum_{\ell=1}^m\sum_{i=1}^n\sum_{j=1}^n w_{ij\ell}\ \mathbf{tr}\ (X\Lambda_\ell-\xi_{ij\ell}(X_\ell^{(k)})X_\ell^{k})'A_{ij}(X\Lambda_\ell-\xi_{ij\ell}(X_\ell^{(k)})X_\ell^{k}), $$ with $$ \xi_{ij\ell}(X_\ell)=\frac{\delta_{ij\ell}}{d_{ij}(X_\ell)} $$ Now piecewise MDS with $m$ groups is the `smacof` algorithm for INDSCAL, in which we optimize over the $\Lambda_\ell$ and use the same common configuration for each group. Again, this is easily extended to matrices $\Lambda_\ell$ which are not restricted to be diagonal, thus giving the `smacof` algorithm for IDIOSCAL. We use the color data of @helm_59 to illustate the `ismacof` function. There are 12 individuals, the first ten with normal color vision (N), the last two with color deficiency (CD). For individual N6 and individual CD2 the experiment was replicated. ```{r helm_data, echo = FALSE} data (helm, package = "smacof") xnms <- rownames(as.matrix(helm[[1]])) lnms <- attr (helm, "names") delta <- array (0, c(10, 10, 16)) for (i in 1:16) delta[, , i] <- as.matrix (helm[[i]]) ``` ```{r helm_analysis, echo = FALSE, cache = TRUE} h <- ismacof (delta) ``` After `r h$itel` iterations we find stress `r h$stress`. In figure `r figure_nums("helm_x", display = "n")` we give the common configuration $X$, and in figure `r figure_nums("helm_xi", display = "n")`. Note in figure `r figure_nums("helm_xi", display = "n")` the color deficient subjects are in blue, the color normal subjects are in red. In a real analysis of these data we would probably conclude that two dimensions was not enough to show important variation. The stress measures for each of the 16 individuals, including two replications, are ```{r helm_svec, echo = FALSE} names(h$svec) <- lnms mprint (h$svec) ``` and thus stress is much larger for color deficient individuals.

```{r helm_xplot, fig.align = "center", echo = FALSE} par(pty = "s") plot (h$x, type = "n", xlab = "dim 1", ylab = "dim 2", xlim = c(-7, 7), ylim = c(-7, 7)) text (h$x, xnms, col = "RED") ```

```{r helm_lplot, fig.align = "center", echo = FALSE} plot (h$lambda, xlab = "dim 1", ylab = "dim 2", type = "n") text (h$lambda[1:11, ], lnms[1:11], col = "RED") text (h$lambda[12:16, ], lnms[12:16], col = "BLUE") ```

#Modeling Asymmetry ##Slide-vector The slide vector model for asymmetric dissimilarities was proposed, but never published, by Kruskal around 1973. It was discussed in detail by @zielman_heiser_93. It assumes, in our notation, that $X_{ij} = X + \frac12(e_i-e_j)z'$ and $X_{ji} = X + \frac12(e_j-e_i)z' = X - \frac12(e_i-e_j)z'$, so that $\frac12(X_{ij}+X_{ji})=X$ and $(X_{ij}-X_{ji})=(e_i-e_j)z'$. In terms of the squared distances $$ d_{ij}^2(X_{ij})=\mathbf{tr}\ X_{ij}'A_{ij}^{\ }X_{ij}^{\ }=(x_i-x_j+z)'(x_i-x_j+z). $$ Define $$T:=\begin{bmatrix}X\\z\end{bmatrix}$$ and $u_{ij}$ is a vector of length $n+1$ with zeroes, except for elements $i$ and $n+1$, which are $+1$, and element $j$, which is $-1$. Thus $T'u_{ij}=(x_i-x_j+z)$ and $d_{ij}^2(X_{ij})=u_{ij}'TT'u_{ij}$. It follows that a `smacof` step is \begin{equation}\label{E:slide} T^{(k+1)}=\left\{\sum_{i=1}^n\sum_{j=1}^nw_{ij}u_{ij}u_{ij}'\right\}^+\left\{\sum_{i=1}^n\sum_{j=1}^nw_{ij}\xi_{ij}(X_{ij}^{(k)})u_{ij}u_{ij}'\right\}T^{(k)}, \end{equation} with $$ \xi_{ij}(X_{ij}^{(k)})=\frac{\delta_{ij}}{d_{ij}(X_{ij}^{(k)})}. $$ This can be clarified somewhat, by using $$ u_{ij}u_{ij}'=\begin{bmatrix}A_{ij}&(e_i-e_j)\\(e_i-e_j)'&1\end{bmatrix}, $$ so that, for example, $$ \sum_{i=1}^n\sum_{j=1}^nw_{ij}u_{ij}u_{ij}'=\begin{bmatrix}V&\sum_{i=1}^n\sum_{j=1}^n w_{ij}(e_i-e_j)\\\sum_{i=1}^n\sum_{j=1}^nw_{ij}(e_i-e_j)'&\sum_{i=1}^n\sum_{j=1}^n w_{ij}\end{bmatrix}. $$ Note that if $W=\{w_{ij}\}$ is symmetric, then $\sum_{i=1}^n\sum_{j=1}^nw_{ij}(e_i-e_j)=0$, which simplifies the matrix generalized inverse in $\eqref{E:slide}$. We apply the function `ssmacof` to the brand switching of nine green and seven blended bottled teas given by @tanioka_yadohisa_16. The frequencies of switching $n_{ij}$ are transformed to asymmetric dissimilarities with $\delta_{ij}=\sqrt{n_{ii}+n_{jj}-2n_{ij}}$. ```{r slide_vector_data, echo = FALSE} tea <- structure(list(DG = c(32L, 2L, 2L, 0L, 0L, 2L, 0L, 3L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L), IG1 = c(0L, 283L, 41L, 0L, 0L, 0L, 0L, 0L, 4L, 15L, 0L, 0L, 0L, 0L, 0L, 0L), IG2 = c(0L, 0L, 177L, 0L, 0L, 0L, 0L, 0L, 5L, 11L, 0L, 0L, 0L, 0L, 0L, 0L), IG3 = c(0L, 18L, 9L, 18L, 0L, 0L, 0L, 0L, 0L, 5L, 0L, 0L, 0L, 0L, 0L, 0L), KaG = c(0L, 2L, 0L, 2L, 58L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), SaG = c(0L, 2L, 2L, 0L, 0L, 24L, 1L, 3L, 0L, 0L, 0L, 3L, 5L, 0L, 0L, 0L), SG1 = c(0L, 1L, 1L, 3L, 1L, 0L, 14L, 8L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L), SG2 = c(0L, 40L, 13L, 9L, 3L, 0L, 0L, 246L, 3L, 11L, 0L, 0L, 0L, 2L, 0L, 0L), `7G` = c(0, 0, 0, 0, 0, 0, 0, 0, 56, 4, 0, 0, 0, 0, 0, 0), ABl = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 137, 0, 0, 0, 0, 0, 0), CBl1 = c(0, 0, 0, 1, 0, 0, 0, 2, 0, 2, 22, 4, 7, 0, 0, 0), CBl2 = c(0, 5, 2, 2, 1, 0, 1, 2, 1, 5, 0, 114, 22, 3, 0, 0), CBl3 = c(0, 18, 10, 8, 1, 0, 2, 26, 1, 16, 0, 0, 201, 2, 0, 0), KaBl = c(0, 0, 0, 0, 6, 0, 0, 0, 0, 2, 0, 0, 0, 28, 0, 0), KBl1 = c(1, 17, 6, 9, 0, 2, 8, 21, 2, 5, 3, 13, 23, 1, 192, 0), KBl2 = c(0, 2, 0, 0, 1, 0, 0, 1, 0, 2, 1, 7, 5, 2, 13, 5)), .Names = c("DG", "IG1", "IG2", "IG3", "KaG", "SaG", "SG1", "SG2", "7G", "ABl", "CBl1", "CBl2", "CBl3", "KaBl", "KBl1", "KBl2"), row.names = c("DG", "IG1", "IG2", "IG3", "KaG", "SaG", "SG1", "SG2", "7G", "ABl", "CBl1", "CBl2", "CBl3", "KaBl", "KBl1", "KBl2"), class = "data.frame") tea <- as.matrix (tea) nms <- row.names (tea) mprint (tea, 0, 3) ``` ```{r slide_vector_run, echo = FALSE, cache = TRUE} delta <- sqrt (outer (diag (tea), diag (tea), "+") - 2 * tea) h <- ssmacof (delta, itmax = 10000) ``` The slide vector algorithm takes `r h$itel` iterations to find a solution with stress `r h$stress`. The slide vector $z$ is `r mprint(h$z)`, and the configuration ```{r silde_vector_conf, echo = FALSE} row.names (h$x) <- nms mprint (h$x, w = 6, d =4, f = "+") ``` is plotted in figure `r figure_nums("slide_vector_tea", display = "n")` (green teas in blue, blended teas in red).

```{r slide_vector_fig, echo = FALSE, fig.align = "center"} plot (h$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h$x[1:9, ], nms[1:9], col = "BLUE") text (h$x[10:16, ], nms[10:16], col = "RED") ```

##Asymmetric Scaling of Dimensions Another possible set of constraints that can be used for asymmetric dissimilarities is $X_{ij}=X\Lambda_j$. Thus $X_{ji}=X\Lambda_i$, and $X_{ij}+X_{ji}=X(\Lambda_i+\Lambda_j)$ while $X_{ij}-X_{ji}=X(\Lambda_i-\Lambda_j)$. In terms of the squared distances $$ d_{ij}^2(X_{ij})=\sum_{s=1}^p\lambda_{js}^2(x_{is}-x_{js})^2. $$ Not surprisingly, the stationary equation are very similar to equations $\eqref{E:stat1}$ and $\eqref{E:stat2}$ in the section on piecewise MDS. We can use the iterative algorithm $$ \sum_{j=1}^mV_jX^{(k+1)}(\Lambda_j^{(k)})^2=\sum_{j=1}^mB_j(X^{(k)}_{ij})X^{(k)}(\Lambda_j^{(k)})^2, $$ $$ \Lambda_j^{(k+1)}=\frac{\mathbf{diag}((X^{(k+1)})'B_j(X^{(k)}_{ij})X^{(k)}_{ij})}{\mathbf{diag}((X^{(k+1)})'V_jX^{(k+1)})}. $$ Here \begin{align*} V_j&=\sum_{i=1}^n w_{ij}A_{ij},\\ B_j(X_{ij}^{(k)})&=\sum_{i=1}^n w_{ij}\xi_{ij}(X_{ij}^{(k)})A_{ij}. \end{align*} The solution, implemented in `zsmacof` is virtually the same as the solution for $X$ and $\Lambda$ from $\eqref{E:stat1}$ and $\eqref{E:stat2}$. It uses alternating least squares in exactly the same way. ```{r dimscal_run, echo = FALSE, cache = TRUE} delta <- sqrt (outer (diag (tea), diag (tea), "+") - 2 * tea) h <- zsmacof (delta, itmax = 10000) ``` We apply the method to the tea data. The algorithm now takes `r h$itel` iterations to find a solution with stress `r h$stress`. The weights $\Lambda_j$ are ```{r dimscal_weights, echo = FALSE} mprint (h$lambda, w = 6, d =4, f = "+") ``` and they are plotted in figure `r figure_nums("dimscal_fig_weights", display = "n")`.

```{r dimscal_fig_weights, echo = FALSE, fig.align = "center"} plot (h$lambda, type = "n", xlab = "dim 1", ylab = "dim 2") text (h$lambda[1:9, ], nms[1:9], col = "BLUE") text (h$lambda[10:16, ], nms[10:16], col = "RED") ```

The configuration is ```{r dimscal_conf, echo = FALSE} row.names (h$x) <- nms mprint (h$x, w = 6, d =4, f = "+") ``` It is plotted in figure `r figure_nums("dimscal_fig_conf", display = "n")`.

```{r dimscal_fig_conf, echo = FALSE, fig.align = "center"} plot (h$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h$x[1:9, ], nms[1:9], col = "BLUE") text (h$x[10:16, ], nms[10:16], col = "RED") ```

This particular approach to asymmetry may be most interesting in unfolding, where we have $$ d_{ij}^2(X_{ij},Y_{ij})=\sum_{s=1}^p\lambda_{js}^2(x_{is}-y_{js})^2, $$ and where the asymmetric treatment of rows and columns is more natural. #Appendix: Code The code for the six different versions of anarchic MDS is both redundant and relatively inefficient. All six programs have the same structure and, for example, share the same voluminous code for debugging I/O. ##auxilary.R ```{r file_auxilary, code = readLines("auxilary.R")} ``` ##asmacof.R ```{r file_asmacof, code = readLines("asmacof.R")} ``` ##bsmacof.R ```{r file_bsmacof, code = readLines("bsmacof.R")} ``` ##psmacof.R ```{r file_psmacof, code = readLines("psmacof.R")} ``` ##ismacof.R ```{r file_ismacof, code = readLines("ismacof.R")} ``` ##ssmacof.R ```{r file_ssmacof, code = readLines("ssmacof.R")} ``` ##zsmacof.R ```{r file_zsmacof, code = readLines("zsmacof.R")} ``` #References