setwd("C:\\Users\\authorized user\\Dropbox\\STAT445\\week11\\tutorials")

euromat = as.matrix(eurodist)

###############################################################
#classic Multidimensional scaling (MDS) ising 'cmdscale'
#with 2 dimensions
mds = cmdscale(eurodist, k=2,eig=TRUE)


x=mds$points[,1]
y= mds$points[,2]

#the scaling represents geographical map of European cities.
#however Athens is on North and Stockholm is on South
#they should be reversed
#So we are inverting the y axis to obtain more realistic map
#of the Eropean cities
plot(x,y, xlab="Coordinate 1", ylab="Coordinate 2", 
     main="Metric  MDS")
text(x, y, labels = labels(eurodist), cex=.9,col=rainbow(10))


# by adding the argument ylim=rev(range(y)) in the plot
#function we can obtain the realisitc map 
#xlim and ylim arguments are helpfull when you want
#to sep up the limits for each of the axes x and y
plot(x,y, xlab="Coordinate 1", ylab="Coordinate 2", 
     main="Metric  MDS",ylim=rev(range(y)))
text(x, y, labels = labels(eurodist), cex=.9,col=rainbow(10))

#apply the Mardia et.al (1979) criterion to see if
#the 2 -dimensional solution fit well
(pm1=cumsum(abs(mds$eig))/sum(abs(mds$eig)))
(pm2=cumsum((mds$eig)^2)/sum((mds$eig)^2))
###############################################################



###############################################################
#Multidimensional scaling (MDS) using 'cmdscale'
#with 3 dimensions
library(scatterplot3d)
mds3 = cmdscale(eurodist, k=3,eig=T)

x= mds3$points[,1]
y= mds3$points[,2]
z= mds3$points[,3]

sp=scatterplot3d(x,y,z,pch=19,color="blue",type="h")
# convert 3D coords to 2D projection
sp.coords <- sp$xyz.convert(mds3$points)
#add text labels
text(sp.coords$x, sp.coords$y,labels = labels(eurodist), cex=.9,col=rainbow(10))



sp=scatterplot3d(x,y,z,pch=19,color="blue",type="h",ylim=rev(range(y)))
# convert 3D coords to 2D projection
sp.coords <- sp$xyz.convert(mds3$points)
#add text labels
text(sp.coords$x, (sp.coords$y),labels = labels(eurodist), cex=.9,col=rainbow(10))


#apply the Mardia et.al (1979) criterion to see if
#the 2 -dimensional solution fit well
(pm1=cumsum(abs(mds3$eig))/sum(abs(mds3$eig)))
(pm2=cumsum((mds3$eig)^2)/sum((mds3$eig)^2))
###############################################################

###############################################################
#nonmetric MDS
#dim=2
#function isoMDS uses Kruskal nonmetric MDS method
#base on stress function 
#so results returned from isoMDS function
#includes coordinates (in the list object points)
#and a value for calculated stress function
#requires the MASS library

library(MASS)
nMDS2=isoMDS(d=eurodist,k = 2)

x=nMDS2$points[,1]
y=nMDS2$points[,2]

plot(x,y, xlab="Coordinate 1", ylab="Coordinate 2", 
     main="Nonmetric  MDS")
text(x, y, labels = labels(eurodist), cex=.9,col=rainbow(10))


plot(x,y, xlab="Coordinate 1", ylab="Coordinate 2", ylim=rev(range(y)),
     main="Nonmetric  MDS")
text(x, y, labels = labels(eurodist), cex=.9,col=rainbow(10))

nMDS2$stress

#plot of the stress function against the MDS dimension
stress_vector=rep(NA,7)

stress_vector[2]=nMDS2$stress/100

for (i in (3:10))
{
  stress_vector[i]=isoMDS(d=eurodist,k = i)$stress/100
}
plot(1:10,stress_vector,type="l",xlab="q",ylab="Stress",main="Stress function")

###############################################################