#############################################################################
######################                                                      #           
# R example code for multidimensional scaling and correspondence analysis:  #
######################                                                      #
#############################################################################


# The built-in R data set called eurodist contains a "dist" object 
# with road distances between 21 European cities.
# For information, type:

help(eurodist)

# To see the (partial) distance matrix, type:

eurodist

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

######################################
# Classical multidimensional scaling #
######################################

# Classical MDS for the European distances using the cmdscale function:
# The default number of dimensions is k = 2:

euro.mds.2 <- cmdscale(eurodist, eig=T)
euro.mds.2

# The first number in the GOF section is what I called P_k, which is 
# (almost) the P^(1) criterion given near the top of page 96.

# Changing the number of dimensions k:

euro.mds.4 <- cmdscale(eurodist, k=4, eig=T)
euro.mds.4

# Using the P_k criterion for a variety of values of k
# to choose the appropriate amount of dimension reduction:

# If we are considering values of k from 1 to, say, 12:

max.k <- 12
P.k <- rep(0,max.k)
SStress <- rep(0,max.k)
for (kk in 1:max.k){
my.mds.kk <- cmdscale(eurodist,k=kk,eig=T)
P.k[kk] <- my.mds.kk$GOF[1]
#SStress[kk] <- ( sum( (eurodist^2 - (dist(my.mds.kk$points))^2)^2 )/sum(eurodist^4) )^0.5
}
cbind(1:max.k,P.k)


# It looks like 2 or 3 dimensions would be reasonable, and 4 gives quite a good fit.

# A 2-D representation of the solution for k=2:

par(pty="s")  # Creates a square plotting region

## The reason for the (-2500, 2500) limits below is because the coordinates in the first two dimensions,
## which can be seen by typing:  euro.mds.2$points
## range from -2048.449113 to 2290.274680 
## Making the axes go from -2500 to 2500 allows all the labels to fit nicely on the plot.
## For other data sets, you'll have to adjust these numbers.
## You can look at the '$points' component of the output and see the smallest and largest values
## and use limits that add a little space beyond this range.

plot(euro.mds.2$points[,1], euro.mds.2$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(-2500,2500), ylim=c(-2500,2500) )
text(euro.mds.2$points[,1], euro.mds.2$points[,2], labels=labels(eurodist) )

# Same thing, but using abbreviations:

euro.abb <- abbreviate(labels(eurodist))

plot(euro.mds.2$points[,1], euro.mds.2$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(-2500,2500), ylim=c(-2500,2500) )
text(euro.mds.2$points[,1], euro.mds.2$points[,2], labels=euro.abb )

# Maybe rotating across the x-axis would produce a better reflection of reality!

plot(euro.mds.2$points[,1], -euro.mds.2$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(-2500,2500), ylim=c(-2500,2500) )
text(euro.mds.2$points[,1], -euro.mds.2$points[,2], labels=euro.abb )

### A 3-D plot of the 3-D solution, including labels of the city names:

library(scatterplot3d)

   s3d <- scatterplot3d(euro.mds.4$points[,1],euro.mds.4$points[,2],euro.mds.4$points[,3],        # x y and z axis
                 color="blue", pch=19,        # filled blue circles
                 type="h",                    # vertical lines to the x-y plane
                 main="3-D Scatterplot of 3-D MDS Solution",
                 xlab="Coordinate 1",
                 ylab="Coordinate 2",
                 zlab="Coordinate 3")
s3d.coords <- s3d$xyz.convert(euro.mds.4$points[,1],euro.mds.4$points[,2],euro.mds.4$points[,3]) # convert 3D coords to 2D projection
text(s3d.coords$x, s3d.coords$y,             # x and y coordinates
    labels=labels(eurodist),               # text to plot
    cex=.7, pos=4)           # shrink text and place to right of points)

### This plot works well.


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


# Cola dissimilarity matrix:
# [Data from Table 5.1 (Subject 1)]


cola.diss <- matrix(c(
0,16,81,56,87,60,84,50,99,16,
16,0,47,32,68,35,94,87,25,92,
81,47,0,71,44,21,98,79,53,90,
56,32,71,0,71,98,57,73,98,83,
87,68,44,71,0,34,99,19,52,79,
60,35,21,98,34,0,99,92,17,44,
84,94,98,57,99,99,0,45,99,24,
50,87,79,73,19,92,45,0,84,18,
99,25,53,98,52,17,99,84,0,98,
16,92,90,83,79,44,24,18,98,0
),ncol=10,nrow=10,byrow=T)


# Using the P_k criterion for a variety of values of k
# to choose the appropriate amount of dimension reduction:

# If we are considering values of k from 1 to, say, 8:

max.k <- 8
P.k <- rep(0,max.k)
SStress <- rep(0,max.k)
for (kk in 1:max.k){
my.mds.kk <- cmdscale(cola.diss,k=kk,eig=T)
P.k[kk] <- my.mds.kk$GOF[1]
#SStress[kk] <- ( sum( (cola.diss^2 - (dist(my.mds.kk$points))^2)^2 )/sum(cola.diss^4) )^0.5
}
cbind(1:max.k,P.k)

# Plotting SStress against k for several values of k:

#plot(1:max.k, SStress, type='b')


# It looks like we need at least 4 dimensions,
# certainly no more than 5.

cola.mds.4 <- cmdscale(cola.diss,k=4,eig=T)
cola.mds.4

# A representation of the 2-dimensional solution:


par(pty="s")  # Creates a square plotting region

plot(cola.mds.4$points[,1], cola.mds.4$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(-55,55), ylim=c(-55,55) )
text(cola.mds.4$points[,1], cola.mds.4$points[,2])

# Interesting:  Note the positions of colas 1, 8, and 10 on the map
# and examine the pairwise dissimilarities for these three colas...

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


# Reading in the College data set from the 'ISLR' packages:

# install.packages("ISLR")  # if needed

library(ISLR)
data(College)

LargeColleges <- subset(College, Apps >= 7600)

coll.abbs<- abbreviate(row.names(LargeColleges))

LargeColleges.numeric <- LargeColleges[,-1]
attach(LargeColleges.numeric)



std <- apply(LargeColleges.numeric,2,sd) # finding standard deviations of variables

# dividing each variable by its standard deviation:
LargeColleges.std <- sweep(LargeColleges.numeric,2,std,FUN="/")  

# Calculating Euclidean distances between these 80 scaled observations:

coll.dist <- dist(LargeColleges.std)


# Using the P_k criterion for a variety of values of k
# to choose the appropriate amount of dimension reduction:

# If we are considering values of k from 1 to, say, 6:

max.k <- 6
P.k <- rep(0,max.k)
SStress <- rep(0,max.k)
for (kk in 1:max.k){
my.mds.kk <- cmdscale(coll.dist,k=kk,eig=T)
P.k[kk] <- my.mds.kk$GOF[1]
SStress[kk] <- ( sum( (coll.dist^2 - (dist(my.mds.kk$points))^2)^2 )/sum(coll.dist^4) )^0.5
}
cbind(1:max.k,P.k)

# Plotting SStress against k for several values of k:
# This makes more sense when the MDS is based on Euclidean distances:

plot(1:max.k, SStress, type='b')

# Using 4 or 5 dimensions seems appropriate here.

coll.mds.5 <- cmdscale(coll.dist,k=5,eig=T)
coll.mds.5

# A representation of the 2-dimensional solution:

par(pty="s")  # Creates a square plotting region

## The value of 'cex' controls how large the labels are printed on the plot.

plot(coll.mds.5$points[,1], coll.mds.5$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(-7,7), ylim=c(-7,7) )
text(coll.mds.5$points[,1], coll.mds.5$points[,2], labels=coll.abbs, cex=0.7)

# Zooming in on the cluster of points in the middle:

plot(coll.mds.5$points[,1], coll.mds.5$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(-7,7), ylim=c(-4,4) )
text(coll.mds.5$points[,1], coll.mds.5$points[,2], labels=coll.abbs, cex=0.7)

# Is it better to use the full row names as labels?

# The Zoomed-in plot (doesn't show Rutgers):

plot(coll.mds.5$points[,1], coll.mds.5$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(-7,7), ylim=c(-4,4) )
text(coll.mds.5$points[,1], coll.mds.5$points[,2], labels=row.names(LargeColleges), cex=0.5)

## Zooming in even more to see our favorite university:

plot(coll.mds.5$points[,1], coll.mds.5$points[,2], type='n',
     xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(0,4), ylim=c(-2,0) )
text(coll.mds.5$points[,1], coll.mds.5$points[,2], labels=row.names(LargeColleges), cex=0.7)

### A 3-D plot of the 3-D solution, including labels of the college names:

library(scatterplot3d)

   s3d <- scatterplot3d(coll.mds.5$points[,1],coll.mds.5$points[,2],coll.mds.5$points[,3],        # x y and z axis
                 color="blue", pch=19,        # filled blue circles
                 type="h",                    # vertical lines to the x-y plane
                 main="3-D Scatterplot of 3-D MDS Solution",
                 xlab="Coordinate 1",
                 ylab="Coordinate 2",
                 zlab="Coordinate 3")
s3d.coords <- s3d$xyz.convert(coll.mds.5$points[,1],coll.mds.5$points[,2],coll.mds.5$points[,3]) # convert 3D coords to 2D projection
text(s3d.coords$x, s3d.coords$y,             # x and y coordinates
    labels=row.names(LargeColleges),               # text to plot
    cex=.5, pos=4)           # shrink text 50% and place to right of points)

### Honestly there are too many observations for this plot to work well.

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

library(MASS)  # loading the MASS package

# Nonmetric scaling on the cola distances:

cola.iso<-isoMDS(cola.diss, k = 2)

plot(cola.iso$points, type = "n",xlab="Coordinate 1", ylab="Coordinate 2")
text(cola.iso$points, labels = as.character(1:nrow(cola.diss)))
abline(h=0); abline(v=0)


###############################
# Correspondence analysis     #
###############################

# Several contingency tables from the book:

girls <- matrix(c(
21,21,14,13,8,
8,9,6,8,2,
2,3,4,10,10),byrow=T, ncol=5, nrow=3,
dimnames = list(c('nbf','bfns','bfs'),
                               c('AG1', 'AG2', 'AG3', 'AG4', 'AG5')))

smokemoms <- matrix(c(
50,315,24,4012,
9,40,6,459,
41,147,14,1594,
4,11,1,124), byrow=T, ncol=4, nrow=4,
dimnames = list(c('YNS','YS','ONS','OS'),
                               c('pd','pa','ftd', 'fta')))

hodgkin <- matrix(c(
74,18,12,
68,16,12,
154,54,58,
18,10,44), byrow=T, ncol=3, nrow=4,
dimnames = list(c('LP','NS','MC','LD'),
                               c('positive','partial','none')))

# One that's not in the book --- a snoring / heart disease data set:

snore.heart.data <- matrix(c(24,35,51,1355,603,416), nrow=2, ncol=3, byrow=TRUE, 
dimnames = list(c("Yes", "No"), c("Never", "Occasionally", "~ Every Night")))

as.table(snore.heart.data)

#########################################################
## Code to get the matrices of chi-squared distances 
## for both rows and columns 
#########################################################
#################  BEGIN CODE ###########################

# Enter name of data matrix / contingency table:
mymat <- hodgkin
rr <- nrow(mymat); cc <- ncol(mymat)
row.sums <- apply(mymat,1,sum)
col.sums <- apply(mymat,2,sum)
N<-sum(row.sums)

pijrow <- matrix(0,nrow=rr,ncol=cc); pijcol <- matrix(0,nrow=rr,ncol=cc)
distmat.row <- matrix(0,nrow=rr,ncol=rr); distmat.col <- matrix(0,nrow=cc,ncol=cc)

for (i in 1:rr){
pijcol[i,] <- mymat[i,]/row.sums[i]
}
for (j in 1:cc){
pijrow[,j] <- mymat[,j]/col.sums[j]
}

for (i in 1:rr){
for (ii in 1:(i-1)) {
my.hold <- sqrt( sum( (N/col.sums)*(pijcol[i,]-pijcol[ii,])^2 ) )
distmat.row[i,ii] <- my.hold; distmat.row[ii,i] <- my.hold
}
}

for (j in 1:cc){
for (jj in 1:(j-1)) {
my.hold <- sqrt( sum( (N/row.sums)*(pijrow[,j]-pijrow[,jj])^2 )  )
distmat.col[j,jj] <- my.hold; distmat.col[jj,j] <- my.hold
}
}

round(distmat.col,digits=2)
round(distmat.row,digits=2)

###################  END CODE ###########################
#########################################################

# A straightforward chi-square test for independence of rows and columns:

chisq.test(hodgkin)

# We see that histological classification and response to treatment are clearly associated.

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

# A quick way to do correspondence analysis using the "corresp" function
# in the MASS library:

library(MASS) # loading the MASS package

corresp(hodgkin,nf=2)  # The two-dimensional solution.

biplot(corresp(hodgkin,nf=2)); abline(h=0); abline(v=0)

# It apparently uses a slightly different algorithm than the book's method:
# While the results are close, they are not identical.

# Some conclusions:

# LP and NS have very similar row profiles.
# "No response" occurs especially often with LD.
# "Positive response" is associated with LP and with NS.
# "Partial response" is linked with MC.

# Examining the 1-dimensional solution can be useful in interpretation as well:

corresp(hodgkin,nf=1) 

# Note "LD" and "none" both have large positive coordinates.
# Note "LP" (and "NS") and "positive" both have large negative coordinates.
# But "LD" and "positive" have opposite signs.

# To (nearly) replicate the book's examples in section 5.3 and in 5.3.1,
# replace "hodgkin" with "girls" or "smokemoms" in the above code.

# What are the substantive conclusions in the "smoking / motherhood" analysis?