################################################################################
# Multidimensional scaling
# is used to represent an observed proximity matrix geometrically
# classical multidimensional scaling is equivalent to principal components analysis
################################################################################
# Correspondence analysis
# is used to display the associations among a set of categorical variables
################################################################################
#
# Airline distances between 10 US cities
#
airline.dist<-read.table("/Users/yasu/Desktop/maFall2009/Data/airline.txt", header=T)
airline.mds<-cmdscale(as.matrix(airline.dist),k=9,eig=T)
airline.mds$eig
#
# Look if two dimensions are appropriate
#
sum(abs(airline.mds$eig[1:2]))/sum(abs(airline.mds$eig))
sum(airline.mds$eig[1:2]^2)/sum(airline.mds$eig^2)

sum(abs(airline.mds$eig[1]))/sum(abs(airline.mds$eig))
sum(airline.mds$eig[1]^2)/sum(airline.mds$eig^2)
#
# You may need to flip to obtain desired image
#
par(pty="s")
plot(-airline.mds$points[,1],-airline.mds$points[,2],type="n",xlab="Coordinate 1",ylab="Coordinate 2",
xlim=c(-2000,1500),ylim=c(-2000,1500))
text(-airline.mds$points[,1],-airline.mds$points[,2],labels=row.names(airline.dist))
#

#################################################################################
#
# Hair color and eye color of a sample of individuals
#
haireye<-matrix(c(688,116,584,188,4,326,38,241,110,3,343,84,909,412,26,98,48,403,681,81),ncol=5,byrow=TRUE)
#
colnames(haireye) <- c("Hair:Fair","Hair:Red","Hair:Medium","Hair:Dark","Hair:Black")
rownames(haireye) <- c("Eye:Light","Eye:Blue","Eye:Medium","Eye:Dark")
#
ncol<-length(haireye[1,])
nrow<-length(haireye[,1])
n<-sum(haireye)
#
rtot<-apply(haireye,1,sum) # row total
ctot<-apply(haireye,2,sum) # column total
#
# compare column - hair color
#
haireye
xctot<-rbind(ctot,ctot,ctot,ctot)
xctot
xctot<-haireye/xctot
xctot
rdot<-rtot/n
rdot
dcols<-matrix(0,ncol,ncol)
for(i in 1:ncol){
 for(j in 1:ncol){
  d<-0
  for(k in 1:nrow) d<-d+((xctot[k,i]-xctot[k,j])^2/rdot[k])
  dcols[i,j]<-sqrt(d)
 }
}
dcols
#
# compare row - eye color
#
haireye
xrtot<-cbind(rtot,rtot,rtot,rtot,rtot)
xrtot<-haireye/xrtot
xrtot
cdot<-ctot/n
cdot
drows<-matrix(0,nrow,nrow)
for(i in 1:nrow){
 for(j in 1:nrow){
  d<-0
  for(k in 1:ncol) d<-d+(xrtot[i,k]-xrtot[j,k])^2/cdot[k]
  drows[i,j]<-sqrt(d)
 }
}
drows
#
# scaling
#
r1<-cmdscale(dcols,eig=TRUE)
r1$points
r1$eig
c1<-cmdscale(drows,eig=TRUE)
c1$points
c1$eig
#
# plotting
#
par(pty="s")
plot(r1$points,xlim=range(r1$points[,1],c1$points[,1])*1.2,ylim=range(r1$points[,2],c1$points[,2])*1.2,type="n",
xlab="Coordinate 1",ylab="Coordinate 2")
text(r1$points,labels=colnames(haireye))
text(c1$points,labels=rownames(haireye),col='red')
abline(h=0,lty=2, lwd=.5)
abline(v=0,lty=2, lwd=.5)
#
haireye
#