## Remark: Depending on the architecture/BLAS some solutions can be rotated differently.

## ----intelligence---------------------------------------------------------------------------------
library("smacof")
idiss <- sim2diss(intelligence[,paste0("T", 1:8)])  
fitrat <- mds(idiss)                                
fitint <- mds(idiss, type = "interval")             
fitord <- mds(idiss, type = "ordinal")              
fitspl <- mds(idiss, type = "mspline")              
round(c(fitrat$stress, fitint$stress, fitord$stress, fitspl$stress), 3)


## ----shepard-plot---------------------------------------------------------------------------------
op <- par(mfrow = c(2,2))
plot(fitrat, plot.type = "Shepard", main = "Ratio MDS")
plot(fitint, plot.type = "Shepard", main = "Interval MDS")
plot(fitord, plot.type = "Shepard", main = "Ordinal MDS")
plot(fitspl, plot.type = "Shepard", main = "Spline MDS")
par(op)


## ----lawlerfit------------------------------------------------------------------------------------
WishD <- sim2diss(wish, method = 7)  
fitWish <- mds(WishD)                


## ----rstarts--------------------------------------------------------------------------------------
set.seed(123)
stressvec <- rep(NA, 100)
fitbest <- mds(WishD, init = "random")
stressvec[1] <- fitbest$stress
for(i in 2:100) {
  fitran <- mds(WishD, init = "random")
  stressvec[i] <- fitran$stress
  if (fitran$stress < fitbest$stress) fitbest <- fitran
}
round(fitbest$stress, 4)



## ----ic-plot--------------------------------------------------------------------------------------
set.seed(123)
icWish <- icExplore(WishD, nrep = 100, returnfit = TRUE)
plot(icWish, main = "IC Plot Wish")


## ----lawler---------------------------------------------------------------------------------------
LawlerD <- sim2diss(Lawler, to.dist = TRUE)
fitLaw <- mds(LawlerD)


## ----rstress--------------------------------------------------------------------------------------
set.seed(123)
rstress <- randomstress(n = 9, ndim = 2, nrep = 500, type = "ratio")


## ----mrstress-------------------------------------------------------------------------------------
bound <- mean(rstress) - 2*sd(rstress)
round(bound, 3)


## ----permlaw--------------------------------------------------------------------------------------
set.seed(123)
permLaw <- permtest(fitLaw, nrep = 500, verbose = FALSE)
permLaw


## ----wenmds---------------------------------------------------------------------------------------
library("MPsychoR")
data("Wenchuan", package = "MPsychoR")
Wdelta <- dist(t(Wenchuan))                                        
fitWen <- mds(Wdelta, type = "interval")        
round(fitWen$stress, 3)


## ----wenperm--------------------------------------------------------------------------------------
set.seed(123)
permWen <- permtest(fitWen, data = Wenchuan, method.dat = "euclidean", 
  nrep = 1000, verbose = FALSE)
permWen


## ----wenperm-plot---------------------------------------------------------------------------------
op <- par(mfrow = c(1,2))
plot(permWen, xlim = c(0.19, 0.30))
perm5 <- quantile(permWen$stressvec, probs = 0.05)    
hist(permWen$stressvec, xlim = c(0.10, 0.40), xlab = "Stress Values", 
     main = "Wenchuan Permutation Histogram")
abline(v = perm5, lty = 2, col = "gray", lwd = 2)             
abline(v = fitWen$stress, col = "gray")       
par(op)


## ----jackmds--------------------------------------------------------------------------------------
data("Rogers", package = "MPsychoR")
RogersSub <- Rogers[,1:16]
RogersD <- dist(t(RogersSub))
fitRogers <- mds(RogersD, type = "ordinal")
jackRogers <- jackmds(fitRogers)
jackRogers



## ----jack-plot------------------------------------------------------------------------------------
plot(jackRogers, legend = TRUE, cex.legend = 0.8, inset = c(-0.3, 0))


## ----bootmds1-------------------------------------------------------------------------------------
set.seed(123)
bootRogers <- bootmds(fitRogers, RogersSub, method.dat = "euclidean",
  nrep = 500)
bootRogers

## ----boot-plot------------------------------------------------------------------------------------
plot(bootRogers)


## ----confmds--------------------------------------------------------------------------------------
fitRogers2 <- mds(RogersD)               
confRogers <- confEllipse(fitRogers2)    

## ----cell-plot------------------------------------------------------------------------------------
plot(confRogers, eps = 0.01, ylim = c(-0.11, 0.11), 
  ell = list(lty = 1, col = "gray"))


## ----simplebi-------------------------------------------------------------------------------------
fitFace <- mds(FaceExp, type = "ordinal")    
ext <- data.frame(PU = FaceScale[,1])
biFace <- biplotmds(fitFace, extvar = ext)   
coef(biFace)


## ----biaxprep-------------------------------------------------------------------------------------
library("calibrate")
biFace2 <- biplotmds(fitFace, extvar = ext, scale = FALSE)           
coef(biFace2)  
PUc <- scale(ext, scale = FALSE)
tm <- seq(floor(min(ext)), ceiling(max(ext)), by = 1)           
tmc_PU <- tm - mean(tm)  
X <- fitFace$conf


## ----simplebi-plot--------------------------------------------------------------------------------
plot(biFace, main = "Biplot Vector Representation", vecscale = 0.8, 
  xlim = c(-1.5, 1.5), vec.conf = list(col = "brown"), pch = 20, cex = 0.5)
plot(fitFace, main = "Biplot Axis Representation", xlim = c(-1.5, 1.5))
abline(h = 0, v = 0, lty = 2, col = "gray")
calPU <- calibrate(coef(biFace2), PUc, tm = tmc_PU, tmlab = tm, Fr = X,
  dp = TRUE, axiscol = "brown", axislab = "PU", labpos = 3, verb = FALSE)


## ----neuralbi-------------------------------------------------------------------------------------
data("NeuralActivity")       
data("NeuralScales")         
NeuralD <- Reduce("+", NeuralActivity)/length(NeuralActivity)
fitNeural <- mds(NeuralD, type = "mspline")
biNeural <- biplotmds(fitNeural, NeuralScales[,1:8])
round(biNeural$R2vec, 3)


## ----neuralbi-plot--------------------------------------------------------------------------------
plot(biNeural, col = "gray", label.conf = list(col = "gray", cex = 0.7),
  vec.conf = list(col = "brown"), main = "Neural Biplot",
  xlim = c(-1.55, 1.55), ylim = c(-1.55, 1.55))


## ----facecons_int---------------------------------------------------------------------------------
fitFace <- mds(FaceExp, type = "ordinal")  
Z <- FaceScale[, c(1,3)]                 
fitFaceC1 <- smacofConstraint(FaceExp, type = "ordinal", 
  constraint = "diagonal", external = Z, constraint.type = "interval", 
  init = fitFace$conf)
round(fitFaceC1$C, 3)



## ----faceint-plot---------------------------------------------------------------------------------
Zc <- scale(Z, scale = FALSE)
X1 <- fitFaceC1$conf
tm <- seq(1, 9, by = 1)
tmc_PU <- tm - mean(Z[,1])
tmc_TS <- tm - mean(Z[,2])
op <- par("xpd" = TRUE, mar = c(5, 4, 7, 5) + 0.1)
plot(fitFaceC1, main = "Constrained MDS (Interval)", xlim = c(-1.4, 1.4), ylim = c(-1, 1))
calPU <- calibrate(c(1,0), Zc[,"PU"], tm = tmc_PU, tmlab = tm, Fr = X1, dp = FALSE, 
  axiscol = "brown", axislab = "PU", labpos = 3, verb = FALSE, shiftvec = c(0, par("usr")[4]), 
  shiftfactor = 1, where = 2, laboffset = c(-0.2, 0.15), cex.axislab = 0.9, reverse = TRUE)
calTS <- calibrate(c(0,1), Zc[,"TS"], tm = tmc_TS, tmlab = tm, Fr = X1, dp = FALSE, 
  axiscol = "brown", axislab = "TS", labpos = 3, verb = FALSE, shiftvec = c(par("usr")[2], 0), 
  shiftfactor = 1, where = 2, laboffset = c(0.2, -0.25), cex.axislab = 0.9)
par(op)


## ----fitface_ord----------------------------------------------------------------------------------
fitFaceC2 <- smacofConstraint(FaceExp, type = "ordinal", 
  constraint = "diagonal", external = Z, constraint.type = "ordinal", 
  init = fitFace$conf)
round(fitFaceC2$C, 3)


## ----trafo-plot-----------------------------------------------------------------------------------
op <- par(mfrow = c(1,2))
ind1 <- order(Z[,1])
plot(Z[ind1, 1], fitFaceC2$external[ind1, 1], type = "s", lwd = 2,
  xlab = "Original Scale", ylab = "Transformed Scale", 
  main = "Pleasant-Unpleasant Transformation")
ind2 <- order(Z[,2])
plot(Z[ind2, 2], fitFaceC2$external[ind2, 2], type = "s", lwd = 2,
  xlab = "Original Scale", ylab = "Transformed Scale", 
  main = "Tension-Sleep Transformation")
par(op)


## ----faceord-plot---------------------------------------------------------------------------------
X2 <- fitFaceC2$conf
Zhat <- fitFaceC2$external
tm <- seq(-0.8, 0.8, by = 0.2)
op <- par("xpd" = TRUE, mar = c(5, 4, 7, 5) + 0.1)
plot(fitFaceC2, main = "Constrained MDS (Ordinal)", xlim = c(-1.4, 1.4), ylim = c(-1, 1))
calPU <- calibrate(c(1,0), Zhat[,"PU"], tm = tm, tmlab = tm, Fr = X2, dp = FALSE,
  axiscol = "brown", axislab = "PU transformed", labpos = 3, verb = FALSE, shiftvec = c(0, par("usr")[4]), 
  shiftfactor = 1, where = 2, laboffset = c(-0.2, 0.15), cex.axislab = 0.9, reverse = TRUE)
calPU <- calibrate(c(0,1), Zhat[,"TS"], tm = tm, tmlab = tm, Fr = X2, dp = FALSE,
  axiscol = "brown", axislab = "TS transformed", labpos = 3, verb = FALSE, shiftvec = c(par("usr")[2], 0),
  shiftfactor = 1, where = 2, laboffset = c(0.25, -0.2), cex.axislab = 0.9)
par(op)


## ----face2b---------------------------------------------------------------------------------------
Z <- FaceScale
fitFaceC3 <- smacofConstraint(FaceExp, type = "ordinal",
  constraint = "unrestricted", external = Z, constraint.type = "interval",
  init = fitFace$conf)
round(fitFaceC3$C, 3)



## ----face2-plot-----------------------------------------------------------------------------------
Zc <- scale(Z, scale = FALSE)
X <- fitFaceC3$conf
biFaceZ <- biplotmds(fitFaceC3, extvar = Z, scale = FALSE)
betas <- coef(biFaceZ)
tm <- seq(1, 9, by = 1)
tmc_PU <- tm - mean(Z[,1])
tmc_AR <- tm - mean(Z[,2])
tmc_TS <- tm - mean(Z[,3])
plot(fitFaceC3, main = "Constraint MDS Biplot", xlim = c(-1.3, 1.3), ylim = c(-1, 1))
abline(h = 0, v = 0, lty = 2, col = "gray")
calPU <- calibrate(betas[,"PU"], Zc[,"PU"], tm = tmc_PU, tmlab = tm, Fr = X,
  dp = FALSE, axiscol = "brown", axislab = "PU", labpos = 3, verb = FALSE)
calAR <- calibrate(betas[,"AR"], Zc[,"AR"], tm = tmc_AR, tmlab = tm, Fr = X,
  dp = FALSE, axiscol = "brown", axislab = "AR", labpos = 3, verb = FALSE)
calTS <- calibrate(betas[,"TS"], Zc[,"TS"], tm = tmc_TS, tmlab = tm, Fr = X,
  dp = FALSE, axiscol = "brown", axislab = "TS", labpos = 3, verb = FALSE)


## ----carconf--------------------------------------------------------------------------------------
library("prefmod")
carconf1 <- carconf[1:100, 1:6]
head(carconf1)


## ----ordinal-unfolding----------------------------------------------------------------------------
unf_ord <- unfolding(carconf1, type = "ordinal") 
unf_ord


## ----carunf-plot----------------------------------------------------------------------------------
op <- par(mfrow = c(1,2))
plot(unf_ord, main = "Unfolding Configuration Car Preferences", 
  ylim = c(-1, 1))
plot(unf_ord, plot.type = "Shepard", 
  main = "Shepard Diagram Car Preferences")
par(op)


## ----row-unfolding--------------------------------------------------------------------------------
startconf <- list(unf_ord$conf.row, unf_ord$conf.col) 
unf_cond <- unfolding(carconf1, type = "ordinal", conditionality = "row", 
  eps = 6e-5, init = startconf)
unf_cond


## ----condunf-plot---------------------------------------------------------------------------------
op <- par(mfrow = c(1,2))
plot(unf_cond, main = "Conditional Unfolding Configuration Car Preferences")
plot(unf_cond, plot.type = "Shepard", 
  main = "Shepard Diagram Car Preferences", col.dhat = "gray", 
  xlim = c(0.9, 6.1))
par(op)


## ----valcirc--------------------------------------------------------------------------------------
tuv <- matrix(NA, nrow = ncol(PVQ40agg), ncol = 2) 
alpha <- -360/10
for (i in 1:10){ 
  alpha <- alpha+360/10
  tuv[i,1] <- cos(alpha*pi/180)
  tuv[i,2] <- sin(alpha*pi/180) 
}


## ----extunf---------------------------------------------------------------------------------------
delta <- (max(PVQ40agg) + 1) - PVQ40agg 
unf_pvq <- unfolding(delta, type = "ordinal")        
unf_pvqext <- unfolding(delta, type = "ordinal", fixed = "column", 
  fixed.coord = tuv)             




## ----unfext-plot----------------------------------------------------------------------------------
op <- par(mfrow = c(1,2))
plot(unf_pvq, label.conf.columns = list(cex = 1), ylim = c(-1.1, 1.1), main = "PVQ Unrestricted")
plot(unf_pvqext, label.conf.columns = list(cex = 1), ylim = c(-1.1, 1.1), main = "PVQ Fixed Value Circle")
par(op)


## ----unfcircle------------------------------------------------------------------------------------
unf_vals <- unfolding(indvalues)                      
unf_valsc <- unfolding(indvalues, circle = "column")  


## ----schwartzcir-plot-----------------------------------------------------------------------------
op <- par(mfrow = c(1,2))
plot(unf_vals, main = "Values Unfolding Configuration (Unrestricted)", ylim = c(-1.5, 1.1))
abline(h = 0, v = 0, col = "gray", lty = 2)
plot(unf_valsc, main = "Values Unfolding Configuration (Circular Restriction)", ylim = c(-1.5, 1.1))
abline(h = 0, v = 0, col = "gray", lty = 2)
rad <- sqrt(sum(unf_valsc$conf.col[1,]^2))
draw.circle(0, 0, radius = rad, border = "cadetblue", lty = 2)
par(op)


## ----vmu------------------------------------------------------------------------------------------
fit_vmu <- vmu(PVQ40agg)        
fit_vmu


## ----vmubi-plot-----------------------------------------------------------------------------------
plot(fit_vmu, col = c("black", "coral3"), cex = c(1, 0.7)) 


## ----plato----------------------------------------------------------------------------------------
PlatoD <- dist(t(Plato7))
fitPlato <- uniscale(PlatoD, verbose = FALSE)  
round(sort(fitPlato$conf), 3)


## -------------------------------------------------------------------------------------------------
gopD <- gravity(GOPdtm)$gravdiss


## ----gravity-plot---------------------------------------------------------------------------------
fitGOP <- mds(gopD, type = "ordinal")    
plot(fitGOP, plot.type = "bubbleplot", bubscale = 20)


## ----drift-plot-----------------------------------------------------------------------------------
morseDrift <- driftVectors(morse2, type = "ordinal") 
plot(morseDrift, main = "Drift Vectors Morse Codes", 
  col.drift = "darkgray")


## ----maryam---------------------------------------------------------------------------------------
artD <- sim2diss(VaziriXu$artificialR)
fitart <- mds(artD, type = "ordinal")    
realD <- sim2diss(VaziriXu$realR)
fitnat <- mds(realD, type = "ordinal")   

## ----maryam-plot----------------------------------------------------------------------------------
op <- par(mfrow = c(1,2))
plot(fitart, main = "Configuration Artificial")
plot(fitnat, main = "Configuration Natural")
par(op)


## ----procnat--------------------------------------------------------------------------------------
fitproc <- Procrustes(fitart$conf, fitnat$conf)
fitproc


## ----proc-plot------------------------------------------------------------------------------------
plot(fitproc, legend = list(labels = c("artificial", "natural")))


## ----stress0--------------------------------------------------------------------------------------
stress0(rectangles, init = rect_constr)  


## ----proc2-plot-----------------------------------------------------------------------------------
fitrect <- mds(rectangles, type = "ordinal", init = rect_constr)
procRect <- Procrustes(rect_constr, fitrect$conf)
plot(procRect, legend = list(labels = c("theoretical", "observed")),
  xlim = c(2, 7))