## ----------------------------------------------------------------------------- library("knitr") library("plot3D") library("smacof") ## ----------------------------------------------------------------------------- load("cola.switch.RData") dis <- sim2diss(as.matrix(cola.switch), method = "counts") dis[is.infinite(dis)] <- 1 res2 <- mds(dis) par(mfrow = c(1, 2)) plot(res2, las = 1, col = "red") plot(res2, plot.type = "Shepard", las = 1, sub = paste("Stress-1 =", formatC(res2$stress, digits = 6))) ## ----------------------------------------------------------------------------- d <- dist(res2$conf) eta.dhat <- sum(res2$dhat^2) eta.X <- sum(d^2) rho.X <- sum(res2$dhat * d) lambda <- rho.X/(eta.dhat * eta.X)^.5 lambda lambda^2 1 - lambda^2 sum((res2$dhat - d)^2)/eta.dhat ## ----------------------------------------------------------------------------- n <- 31 x <- seq(-3, 3, length.out = n) y <- seq(-3, 3, length.out = n) Y <- matrix(rep(x, each = n), nrow = n, ncol = n) X <- matrix(rep(y, n), nrow = n, ncol = n) f <- function(X,Y){ -sqrt((X - 0)^2 + (Y - 0)^2) } Z <- f(X, Y) zlims <- range(Z) Z[Z < -3] <- NA zlims[1] <- zlims[1] - .3 * (zlims[2] - zlims[1]) phi <- 20 theta <- 35 res <- ribbon3D(x = x, y = y, z = Z, clab = "f(x)", shade = 0.2, scale = TRUE, ticktype = "detailed", theta = theta, phi = phi, xlab = expression("x[1]"), ylab = expression("x[2]"), zlab = "f(x)", contour = TRUE, zlim = zlims, cex = 1.5, colkey = FALSE) f.sup <- f(-1, -1) Z <- -(X + Y)/f.sup res <- ribbon3D(x = x, y = y, z = Z, shade = 0.2, scale = TRUE, ticktype = "detailed", theta = theta, phi = phi, contour = TRUE, zlim = zlims, cex = 1.5, colkey = FALSE, add = TRUE) lines3D(x = c(-1, -1), y = c(-1, -1), z = c(zlims[1], f.sup), type = "b", pch = 20, col = "red", add = TRUE) ## ----------------------------------------------------------------------------- res1 <- mds(dis, ndim = 1, verbose = TRUE) ## ----------------------------------------------------------------------------- multistart <- function(dis, n.starts = 1000, seed = NULL, ...){ stress <- rep(0, n.starts + 1) res <- mds(dis, init = "torgerson", ...) stress[n.starts + 1] <- res$stress stress.best <- res$stress x <- res$conf if (!is.null(seed)) set.seed(seed) for (i in 1:n.starts) { res <- mds(dis, init = "random", ...) stress[i] <- res$stress if (stress[i] < stress.best){ stress.best <- stress[i] x <- res$conf } } res <- mds(dis, init = x, ...) return(list(res = res, stress = stress)) } ## ----------------------------------------------------------------------------- res.multi <- multistart(dis, ndim = 1, seed = 1234) x <- res.multi$res$conf par(mfrow = c(2, 2)) plot(rep(0, length(x)), x, type = "p", pch = 20, xaxt = "n", xlab = "", las = 1, ylab = "Dimension 1", xlim = c(-1, 1), col = "red") text(rep(0, length(x)), x, labels = rownames(x), pos = 4) plot(res.multi$res, plot.type = "Shepard", las = 1, sub = paste("Stress-1 =", formatC(res.multi$res$stress, digits = 6))) hist(res.multi$stress, las = 1, xlab = "Stress-1") abline(v = res.multi$res$stress, col = "red") abline(v = res.multi$stress[1], col = "blue") ## ----------------------------------------------------------------------------- res.multi.2 <- multistart(dis, n.starts = 100, ndim = 2, seed = 5678) res.multi.14 <- multistart(dis, n.starts = 100, ndim = 14, seed = 9012) par(mfrow = c(1, 2)) hist(res.multi.2$stress, las = 1, xlab = "Stress-1", main = "", sub = "Dimensionality p = 2") abline(v = res.multi.2$res$stress, col = "red") abline(v = res.multi.2$stress[1], col = "blue") tt <- hist(res.multi.14$stress, las = 1, xlab = "Stress-1", main = "", sub = "Dimensionality p = 14", xlim = c(.15, .16), breaks = c(.1535, .1545)) abline(v = res.multi.14$res$stress, col = "red") abline(v = res.multi.14$stress[1], col = "blue") ## ----------------------------------------------------------------------------- res.procr <- Procrustes(res2$conf, res.multi.2$res$conf) plot(res.procr, las = 1, main = "") ## ----------------------------------------------------------------------------- w <- dissWeights(dis, type = "unif") ind <- order(as.vector(dis)) res.unw <- mds(dis) res.wgt <- mds(dis, weightmat = w) par(mfrow = c(4, 2)) plot(dis[ind], 1:length(dis)/length(dis), type = "l", las = 1, col = "blue", main = "EDF", xlab = "Dissimilarities", ylab = "Proportion") plot(dis[ind], cumsum(w[ind])/sum(w), type = "l", las = 1, col = "blue", main = "Weighted EDF", xlab = "Dissimilarities", ylab = "Proportion") plot(res.unw, plot.type = "histogram", las = 1, main = "Unweighted histogram", border = "blue") plot(res.wgt, plot.type = "histogram", las = 1, main = "Weighted histogram", border = "blue") plot(res.unw, las = 1, col = "red") plot(res.wgt, las = 1, col = "red") plot(res.unw, plot.type = "Shepard", ylim = c(0, 2)) plot(res.wgt, plot.type = "Shepard", ylim = c(0, 2)) ## ----------------------------------------------------------------------------- par(mfrow = c(2, 2)) w <- dissWeights(dis, type = "unifpower", power = 10) w <- w * (n * (n - 1) * 0.5/sum(w^2))^.5 res.5 <- mds(dis, weightmat = w, ndim = 2) plot(res.5, las = 1, main = "Configuration q = 10", col = "red") plot(res.5, plot.type = "Shepard", las = 1, main = "Shepard q = 10") w <- dissWeights(dis, type = "unifpower", power = -5) w <- w * (n * (n - 1) * 0.5/sum(w^2))^.5 res.minus.5 <- mds(dis, weightmat = w, ndim = 2) plot(res.minus.5, las = 1, main = "Configuration q = -5", col = "red") plot(res.minus.5, plot.type = "Shepard", las = 1, main = "Shepard q = -5") ## ----------------------------------------------------------------------------- par(mfrow = c(1, 2)) n <- 500 y <- runif(n, min = -1, max = 1) theta <- runif(n) x <- theta * cos(3 * pi * theta) z <- theta * sin(3 * pi * theta) lims <- c(-1.2, 1.2) points3D(x, y, z, col = "red", pch = 20, cex = 0.5, phi = 10, theta = -10, xlim = lims, ylim = lims, zlim = lims) X <- cbind(x, y, z) dis <- dist(X) w <- dissWeights(dis, type = "knn", k = 15) w <- w * (n * (n - 1) * 0.5/sum(w^2))^.5 res <- mds(dis, weightmat = w, init = "torgerson") plot(res, col = "red", pch = 20, cex = 0.5, las = 1, label.conf = list(label = FALSE)) ## ----------------------------------------------------------------------------- library("network") data("flo") c <- 0.1 dis <- flo w <- flo dis.star <- w + (1 - w)/c w.star <- w + (1 - w) * c res <- mds(dis.star, weightmat = w.star, init = "torgerson", ndim = 2) plot(res, col = "red") n <- nrow(res$conf) k <- 0 w.mat <- as.matrix(w) for (j in 1:(n - 1)){ for (i in (j + 1):n){ if (w.mat[i, j] == 1) { lines(res$conf[c(i, j), 1], res$conf[c(i, j), 2], col = "lightpink") } } }