Final Project - Exploratory Multidimensional Scaling

library("MPsychoR")
library("smacof")
chps = read.csv("chps.csv")

#Random force used to add missing values. missRanger - Fast Imputation of Missing Values by Chained Random Forests

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

data set 㤼㸱chps㤼㸲 not found

mydata2<-missRanger::missRanger(data = chps, verbose=0)

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

ggplot(M) + aes(x=Value, fill=Item, label = ..count..)+geom_bar()+facet_wrap(~Item)

#Standardized the variables since not all of the variables were the same.

standardize = function(x){
  return ((x-mean(x,na.rm=TRUE))/sd(x,na.rm=TRUE))
}

std_data = as.data.frame(apply(mydata2,2,standardize))

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

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

Model: Symmetric SMACOF 
Number of objects: 25 
Stress-1 value: 0.12 
Number of iterations: 49 

plot(fit.wenchuan1, main = "CHPS Survey MDS")

#The stress is higher using the randomstress test.

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

[1] 0.32

hist(x=rsvec, main = "Histogram of randomly generated/nmultidimensional scaling stress")

#The model above is less than 0.31, which means that it is statistically significant.

mean(rsvec) - 2*sd(rsvec)

[1] 0.31

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

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

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

Observed stress value: 0.12 
p-value: <0.001 

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

plot(x=1:(n-1), y=svec, type="b", main = "MDS Scree Plot", pch=20, xlab="Number of Dimensions", ylab = "Stress")

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

[1] 0.1

fit.wenchuan1$stress ## stress (classical scaling start)

[1] 0.12

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

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

Model: Symmetric SMACOF 
Number of objects: 25 
Stress-1 value: 0.18 
Number of iterations: 21 

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

plot(fit.wenchuan1, plot.type = "Shepard",
main = "Shepard Diagram (Orginal MDS)")

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

#Configurations

plot(fit.wenchuan2, plot.type = "confplot", main= "confplot (Random Ordinal MDS)")

plot(fit.wenchuan2, plot.type = "resplot", main= "resplot (Random Ordinal MDS)")

plot(fit.wenchuan2, plot.type = "stressplot", main= "stressplot (Random Ordinal MDS)")

plot(fit.wenchuan2, plot.type = "bubbleplot", main= "bubbleplot(Random Ordinal MDS)")

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

colpal <- c(rainbow_hcl(3, c=100))
pal <-palette(colpal)
memb <- c(1,2,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,1,1,1)
plot(fit.wenchuan2, main = "Configuration", label.conf = list(label=FALSE), col = memb)
legend("topright", legend = c("School Support", "Feeling", "Efficacy"), col=colpal, pch=20)

palette(pal)
colpal <- c(rainbow_hcl(3, c=100, l=30))
pal <- palette(colpal)

text(fit.wenchuan2$conf[-7,], labels=rownames(fit.wenchuan2$conf)[-7], col=memb[-7], pos=3, cex=0.8)
text(fit.wenchuan2$conf[7.1], fit.wenchuan2$conf[7,2], labels= rownames(fit.wenchuan2$conf)[7], col=memb[7], pos=1, cex=0.8)