## ----------------------------- Supplementary Materials Borg & Mair AJS -------------------------------------- require(smacof) ## MDS package require(wordcloud) ## plotting tools ## ----------------- Part I: Torgerson (+ error) starting configuration on artificial grid -------------- ## generate rectangular grid n <- 9 ## number of objects m <- 2 ## number of dimensions x <- matrix(c(0, 0, 1, 0, 2, 0, 0, 1, 1, 1, 2, 1, 0, 2, 1, 2, 2,2 ), nrow = n, ncol = m, byrow = TRUE) ## true common space rownames(x) <- letters[1:n] D2 <- diag(c(1,2)) X2 <- x %*% D2 ## weighting the dimensions plot(X2, asp = 1, pch = 16, xlab = "Dimension 1", ylab = "Dimension 2") ## generate Figure 1 for (i in 1:2) segments(X2[i,1], X2[i,2], X2[i+1,1], X2[i+1,2], lty = 2) for (i in 4:5) segments(X2[i,1], X2[i,2], X2[i+1,1], X2[i+1,2], lty = 2) for (i in 7:8) segments(X2[i,1], X2[i,2], X2[i+1,1], X2[i+1,2], lty = 2) for (i in 1:3) segments(X2[i,1], X2[i,2], X2[i+3,1], X2[i+3,2], lty = 2) for (i in 4:6) segments(X2[i,1], X2[i,2], X2[i+3,1], X2[i+3,2], lty = 2) text(X2, label = rownames(X2), asp = 1, pos = 2, cex = 0.8) ## MDS fit with Torgerson starts (default) delta <- sqrt(dist(X2)) fit1 <- mds(delta, type = "ordinal") ## ordinal 2D MDS fit1 op <- par(mfrow = c(1,2)) ## Figure 2 plot(fit1) plot(fit1, plot.type = "Shepard") par(op) ## MDS fit with Torgerson + error starts nmod <- 3 ## number of models to be fitted ic <- torgerson(delta, p = 2) ## classical MDS fit_list <- list() ## initialize list for model fit storage stressvec <- NULL set.seed(1) for (i in 1:nmod) { z <- matrix(rnorm(n*m, mean = 0, sd = i-1), n, m) ## normal errors (increasing in each loop) ic_error <- ic + z ## adding errors to Torgerson fit_list[[i]] <- mds(delta, type = "ordinal", eps = 1e-20, init = ic_error, itmax = 1e5) stressvec[i] <- fit_list[[i]]$stress } round(stressvec, 2) ## vector with stress values ## plot some of the solutions op <- par(mfrow = c(1,2)) plot(fit_list[[2]]) plot(fit_list[[2]], plot.type = "Shepard") par(op) op <- par(mfrow = c(1,2)) plot(fit_list[[3]]) plot(fit_list[[3]], plot.type = "Shepard") par(op) ## ------------------------------ Part II: local minima analysis on artificial grid ------------------- nrep <- 100 ## number of random starts fit_list <- list() ## initialize list for model fit storage configs <- list() stressvec <- NULL set.seed(3) for (i in 1:nrep) { cat(i, "\n") z <- matrix(runif(n*2, min = 0, max = 1), n, m) fit_list[[i]] <- mds(delta, type = "ordinal", init = z, itmax = 1e5) stressvec[i] <- fit_list[[i]]$stress configs[[i]] <- fit_list[[i]]$conf } round(stressvec, 2) ## compute parwise similarities among all configurations sim_mat <- matrix(0, ncol = nrep, nrow = nrep) ## initialize correlation matrix for (i in 2:nrep){ je <- i-1 for (j in 1:je) { a <- configs[[i]] ## extract configuration pairs b <- configs[[j]] Y.hat <- Procrustes(a,b)$Yhat ## Procrustes transformation sim_mat[i,j] <- sim_mat[j,i] <- cor(c(a), c(Y.hat)) ## compute configuration correlations } } ## checking the similarity of the solutions with MDS and cluster analysis delta <- sim2diss(sim_mat, method = "corr") ## convert correlations into dissimilarities fit2 <- mds(delta, type = "interval") ## fitting an interval MDS plot(fit2, plot.type = "bubbleplot", bubscale = 100*stressvec) ## bubbleplot fitclust <- hclust(delta, method = "ward.D2") ## hierarchical cluster analysis plot(fitclust, xlab = "MDS solutions", sub = "", cex = 0.8) ## dendrogram ## ------------------------------ Part III: real-life applications --------------------- ## --- Wish data data(wish) ## similarities delta <- sim2diss(wish, method = 7) ## convert into dissimilarities fitwish <- mds(delta, type = "ordinal") ## Togerson IC fitwish$stress permtest(fitwish) ## permutation test nrep <- 100 ## number of replications fit_list <- list() ## initialize list for model fit storage configs <- list() stressvec <- NULL set.seed(3) for (i in 1:nrep) { cat(i, "\n") fit_list[[i]] <- mds(delta, type = "ordinal", init = "random") stressvec[i] <- fit_list[[i]]$stress configs[[i]] <- fit_list[[i]]$conf } ## stress distribution based on random starts hist(stressvec, main = "Stress Distribution Wish", xlab = "Stress Values", breaks = 20) min(stressvec) ## slightly smaller than Torgerson start ind <- which.min(stressvec) op <- par(mfrow = c(1,2)) ## plot both configurations plot(fit_list[[ind]], main = "Best Random IC") plot(fitwish, main = "Torgerson IC") par(op) ## Procustes with Torgerson IC fit as target fitproc <- Procrustes(fitwish$conf, fit_list[[ind]]$conf) plot(fitproc) ## the configurations are very similar ## Procedure proposed in the article sim_mat <- matrix(0, ncol = nrep, nrow = nrep) ## initialize correlation matrix for (i in 2:nrep){ je <- i-1 for (j in 1:je) { a <- configs[[i]] ## extract configuration pairs b <- configs[[j]] Y.hat <- Procrustes(a,b)$Yhat ## Procrustes transformation sim_mat[i,j] <- sim_mat[j,i] <- cor(c(a), c(Y.hat)) ## compute configuration correlations } } ## checking the similarity of the solutions with MDS and cluster analysis delta <- sim2diss(sim_mat, method = "corr") ## convert correlations into dissimilarities fit2 <- mds(delta, type = "interval") ## fitting an interval MDS plot(fit2, plot.type = "bubbleplot", bubscale = 100*stressvec) ## bubbleplot fitclust <- hclust(delta, method = "ward.D2") ## hierarchical cluster analysis plot(fitclust, xlab = "MDS solutions", sub = "", cex = 0.8) ## dendrogram ## --- Ekman data data(ekman) ## similarities delta <- sim2diss(ekman, method = 1) ## convert into dissimilarities ## Ordinal MDS with Torgerson IC fitekman <- mds(delta, type = "ordinal") fitekman permtest(fitekman) plot(fitekman) ## Random IC nrep <- 100 ## number of replications fit_list <- list() ## initialize list for model fit storage configs <- list() stressvec <- NULL set.seed(3) for (i in 1:nrep) { cat(i, "\n") fit_list[[i]] <- mds(delta, type = "ordinal", init = "random") stressvec[i] <- fit_list[[i]]$stress configs[[i]] <- fit_list[[i]]$conf } ## Procedure proposed in the article sim_mat <- matrix(0, ncol = nrep, nrow = nrep) ## initialize correlation matrix for (i in 2:nrep){ je <- i-1 for (j in 1:je) { a <- configs[[i]] ## extract configuration pairs b <- configs[[j]] Y.hat <- Procrustes(a,b)$Yhat ## Procrustes transformation sim_mat[i,j] <- sim_mat[j,i] <- cor(c(a), c(Y.hat)) ## compute configuration correlations } } ## checking the similarity of the solutions with MDS and cluster analysis delta <- sim2diss(sim_mat, method = "corr") ## convert correlations into dissimilarities fit2 <- mds(delta, type = "interval") ## fitting an interval MDS plot(fit2, plot.type = "bubbleplot", bubscale = 100*stressvec) ## bubbleplot fitclust <- hclust(delta, method = "ward.D2") ## hierarchical cluster analysis plot(fitclust, xlab = "MDS solutions", sub = "", cex = 0.8) ## dendrogram ## we see that there is one stress outlier