#######################################
### R scripts for Classical MDS
### BEMO 2015
### R. Castilho || CCMAR || 2015-02-17

#######################################
### PREPARATION STEPS
# install packages
install.packages("clusterSim")
install.packages("MASS)
install.packages("calibrate")

# load packages 
library(clusterSim)
library(MASS)
library(calibrate)

#######################################
### LOAD DATA
# set working directory... change this to the directory in which the species' data is stored and where you want to save output
#setwd("D:/UTILIZADORES/user-en/Desktop/BEMO")
setwd("~/Desktop/0_R_MBE")

## Classical MDS
## Reading the data (from a txt file, tab delimited, 1st row: labels, no column labels)
data.df = read.table("data.txt", header=TRUE, sep = "\t")
data.mat = as.matrix(data.df[,1:14]) #transform data.frame into matrix so that can be transposed
dim(data.mat) #X Y
tdata.mat<-t(data.mat) #transposing matrix
dim(tdata.mat) #Y X

## Obtaining distance matrix BC=Bray-Curtis
data.bcdit <- dist.BC(tdata.mat)
fit <- cmdscale(data.bcdist,eig=TRUE, k=2) # k is the number of dim
fit # view results
# plot solution 
x <- fit$points[,1]
y <- fit$points[,2]
plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2", main="Bray-Curtis MDS",   type="p", pch=16)
textxy(x, y, labs = row.names(tdata.mat), cx = 0.5, dcol = "black", m = c(0, 0))

## Delete the Lesseptian species
data.df =data.df[,-14]
data.mat = is.matrix(data.df[,1:13]) #transform data.frame into matrix so that can be transposed
dim(data.mat) #X Y
tdata.mat<-t(data.mat) #transposing matrix
dim(tdata.mat) #Y X

## Obtaining distance matrix BC=Bray-Curtis without LES
data.bcdist <- dist.BC(tdata.mat)
fit <- cmdscale(data.bcdist,eig=TRUE, k=2) # k is the number of dim
fit # view results
# plot solution 
x <- fit$points[,1]
y <- fit$points[,2]
plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2", main="Bray-Curtis MDS",   type="p", pch=16)
textxy(x, y, labs = row.names(tdata.mat), cx = 0.5, dcol = "black", m = c(0, 0))


#Obtaining distance matrix Euclidean distances
d <- dist(tdata.mat) # euclidean distances between the rows
fit <- cmdscale(d,eig=TRUE, k=2) # k is the number of dim

fit # view results
# plot solution 
x <- fit$points[,1]
y <- fit$points[,2]
plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2", main="Euclidean MDS",   type="p", pch=16)
textxy(x, y, labs = row.names(tdata.mat), cx = 0.5, dcol = "black", m = c(0, 0))