Module 5: Multidimensional Scaling

options(digits=2, scipen=999)
data(list="Wenchuan", package="MPsychoR")

Wenchuan <-missRanger::missRanger(data = Wenchuan, verbose=0)

for(j in 1:ncol(Wenchuan))Wenchuan[,j]<-round(Wenchuan[,j])
M <- tidyr::gather(data=Wenchuan,key="Item", value = "Value") 
M$Item <-factor(x=M$Item, levels = names(sort(tapply(X=M$Value, INDEX=M$Item, FUN=mean))))
require(ggplot2)

#Items are reordered and cleaned.

The ggplot below is something I found using YouTube. I love it!

ggplot(M) + aes(x=Value, fill=Item, label = ..count..)+geom_bar()+geom_text(stat="count", position = position_stack(0.9))+facet_wrap(~Item)

library("MPsychoR")
library("smacof")
data(Wenchuan)
Wdelta <- dist(t(Wenchuan)) ## Euclidean distances
fit.wenchuan1 <- mds(Wdelta, type = "ordinal") ## MDS fit
fit.wenchuan1

Call:
mds(delta = Wdelta, type = "ordinal")

Model: Symmetric SMACOF 
Number of objects: 17 
Stress-1 value: 0.133 
Number of iterations: 35 

plot(fit.wenchuan1, main = "Wenchuan MDS")

Here we can see variables that are very similar and others that are very different, looking at the different dimensions.

set.seed(123)
rsvec <- randomstress(n = attr(Wdelta, "Size"), ndim = 2,
nrep = 500, type = "ordinal")
mean(rsvec)

[1] 0.28

#Here is a random stress calculation that gives you a distribution of observations to determine statistical significance. It takes random values and puts it into the function, checks distances, then returns stress. The random stress average is 0.28.

#To detemine a statistical signifcant results, use the following and see values less than it.

mean(rsvec) - 2*sd(rsvec)

[1] 0.25

set.seed(123)
permmds <- permtest(fit.wenchuan1, data = Wenchuan,
method.dat = "euclidean", nrep = 500,
verbose = FALSE)
permmds

Call: permtest.smacof(object = fit.wenchuan1, data = Wenchuan, method.dat = "euclidean", 
    nrep = 500, verbose = FALSE)

SMACOF Permutation Test
Number of objects: 17 
Number of replications (permutations): 500 

Observed stress value: 0.13 
p-value: <0.001 

#Here the p-value is less than 0.001, therefore the model is more informative than a randomly generated reduction. *****

n <- attr(Wdelta, "Size")

v_stress <-rep(NA, n-1)
for(i in 1:(n-1))v_stress[i]<-smacof::mds(delta=Wdelta,ndim=i,type="ordinal")$stress
plot(x=1:(n-1), y=v_stress, type="b", main="MDS Scree Plot", pch=20, xlab="Number of Dimensions", ylab="Stress")

svec <- NULL
for (i in 1:(n-1)) {
svec[i] <- mds(Wdelta, ndim = i, type = "ordinal")$stress
}

#You may choose 2 or 3 dimensions.

set.seed(123)
fit.wenchuan <- NULL ## 100 random starts
for(i in 1:100) {
fit.wenchuan[[i]] <- mds(Wdelta, type = "ordinal",
init = "random")
}
## extract the best solution
ind <- which.min(sapply(fit.wenchuan,
function(x) x$stress))
fit.wenchuan2 <- fit.wenchuan[[ind]]
fit.wenchuan2$stress ## lowest stress (random start)

[1] 0.13

## [1] 0.1327943
fit.wenchuan1$stress ## stress (classical scaling start)

[1] 0.13

## [1] 0.1328058

lowest stress (random start) 0.1327943

stress (classical scaling start) 0.1328058

fit.wenchuan3 <- mds(Wdelta, type = "interval")
fit.wenchuan3

Call:
mds(delta = Wdelta, type = "interval")

Model: Symmetric SMACOF 
Number of objects: 17 
Stress-1 value: 0.18 
Number of iterations: 9 

#Interval will be more restrictive. However, the stress is worst.

plot(fit.wenchuan2, plot.type = "Shepard",
main = "Shepard Diagram (Ordinal MDS)")

plot(fit.wenchuan3, plot.type = "Shepard",
main = "Shepard Diagram (Interval MDS)")

plot(fit.wenchuan2, plot.type = "stressplot",
main = "Wenchuan Stress-per-Point")

#Here the stress plot (which will add up to 100) is being compared to the other variables.

library("colorspace")
set.seed(123)
bootWen <- bootmds(fit.wenchuan2, data = Wenchuan,
method.dat = "euclidean", nrep = 100)
cols <- rainbow_hcl(17, l = 60)
plot(bootWen, col = cols)

#Above you can see variations of each object across multiple samples. It shows how stable the MDS configuration is across multiple samples. The data ellipse shows the distribution of a confidence. The size of the ellipse matters - the larger the ellipse the less confidence we are where it falls within the dimensions.