```{r rectangle_shepard_plot, fig.align = "center", echo = FALSE} par(mfrow = c(1,2), pty = "s") plot(h1$delta, h1$dist, xlab = "dissimilarities", ylab = "distances", col = "RED") plot(h2$delta, h2$dist, xlab = "dissimilarities", ylab = "distances", col = "RED") ```

```{r rectangle_analysis_norm, echo = FALSE} h3 <- smacofShepard62(delta, norm = 2) ``` We also compute the solution using loss function $\eqref{E:KC}$ with $\alpha=0$ instead of $\eqref{E:final}$. After `r h3$itel` iterations we find loss `r h3$stress`. The gradient at the solution is ```{r rectangle_grad_3, echo = FALSE} mprint (h3$grad, f = "+") ``` The solution in figure `r figure_nums("rectangle_shep_norm", display = "n")` is virtually indistinguishable from the one in figures `r figure_nums("rectangle_conf", display = "n")` and `r figure_nums("rectangle_shep", display = "n")`.

```{r rectangle_shepard_plot_norm, fig.align = "center", echo = FALSE} par(mfrow = c(1,2), pty = "s") plot(h3$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h3$x, as.character (1:16), col = "RED") plot(h3$delta, h3$dist, xlab = "dissimilarities", ylab = "distances", col = "RED") ```

##Dutch Political Parties 1967 As the next illustration we use data from @degruijter_67, with average dissimilarity judgments between Dutch political parties in 1967. The data are ```{r poldist-data, 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) poldata <- as.vector (poldist) polnames <- attr (poldist, "Labels") ``` ```{r pol_analysis, echo = FALSE} h <- smacofShepard62 (poldata, verbose = FALSE, eps = 1e-15, itmax=1000) ``` The process converges in `r h$itel` iterations to loss function value `r h$stress`. The gradient at the solution is ```{r pol_grad, echo = FALSE} mprint (h$grad, f = "+") ```

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

```{r pol_plot_shep, fig.align = "center", echo = FALSE} par (pty = "s") plot (h$delta, h$dist, xlab = "dissimilarities", ylab = "distances", col = "RED") ```

##Ekman Color Data The final example are the color data from @ekman_54. ```{r ekman_data, echo = FALSE} data (ekman, package = "smacof") delta <- as.vector ((1 - ekman) ^ 3) colnames <- as.character (attributes(ekman)$Labels) ``` ```{r ekman_analysis, echo = FALSE} h <- smacofShepard62 (delta, verbose = FALSE, eps = 1e-15, itmax=1000) ``` The process converges in `r h$itel` iterations to loss function value `r as.character(h$stress)`. The gradient at the solution is ```{r ekman_grad, echo = FALSE} mprint (h$grad, f = "+") ``` ```{r ekman_plot_conf, fig.align = "center", echo = FALSE} par (pty = "s") plot (h$x, type = "n", xlab = "dim 1", ylab = "dim 2", xlim = c(-.5,.5), ylim = c(-.5,.5)) text (h$x, colnames, col = "RED") ```

```{r ekman_plot_shep, fig.align = "center", echo = FALSE} par (pty = "s") plot (h$delta, h$dist, xlab = "dissimilarities", ylab = "distances", col = "RED") ```

#Discussion As an aside the Shepard loss function is similar to the loss function used in the Guttman-Lingoes programs (@guttman_68a, @lingoes_roskam_73). $$ \sigma_{GL}^\star(x):=\frac12\sum_{k=1}^K(d_k(x)-d_k^\star(x))^2=\sum_{k=1}^K(d_k(x)-d_k^\star(x))d_k(x), $$ where the $d_k^\star(x)$ are the *rank-images*, i.e. the $d_k(x)$ permuted in such a way that they are monotonic with the dissimilarities. Thus instead of `sort(delta)[rank(dist)]` we use `sort(dist)[rank(delta)]`. There is no arithmetic done with dissimilarities, only their rank order is used. But the rank images do not have any clear optimality properties. As a consequence they are considerably less smooth than loss based on monotone regression or on Shepard's rearrangements. The report @deleeuw_R_73g is of interest in this context, in part because the pdf includes a copy of the report agressively annotated by Joseph Kruskal, as well as some obfuscating correspondence with Jim Lingoes. Also see Appendix 2 of @kruskal_77. If one wants to use rank images, it would make more sense to use $$ \sigma^\star_{LG}:=\sum_{k=1}^K(d_k^\star-d_k)\delta_k, $$ because, again, by the rearrangement inequality this is a directionally differentiable difference of two convex functions similar to the one used by Shepard. Shepard himself *evaluates* the departure from monotonicity, after convergence, using $$ \frac12\sum_{k=1}^K(\delta_k-\hat\delta_k)^2=\sum_{k=1}^K(\delta_k-\hat\delta_k)\delta_k=\sum_{k=1}^K(\hat\delta_k-\delta_k)\hat\delta_k, $$ which is again hopelessly non-smooth. #Appendix: Code ##shepard.R ```{r file_auxilary, code = readLines("shepard.R")} ``` ##auxilary.R ```{r file_auxilary, code = readLines("auxilary.R")} ``` ##mdsUtils.R ```{r file_auxilary, code = readLines("mdsUtils.R")} ``` ##smacofShepard62.R ```{r file_auxilary, code = readLines("smacofShepard62.R")} ``` #References