####################################################
######   Exploratory Techniques R Commands   #######
######           UAM Cluster Analysis        #######
####################################################

## Cluster analysis sorts the data relative to a given set of variables. 

## The packages you will need for these techniques are: pvclust, cluster, ape, FactoMineR

	

##############################################################
### Cluster Analysis with bootstrapped p-values in pvclust ###
##############################################################

## Hierarchal Cluster Analysis with bootstrapped p-values {pvclust}

library(pvclust)
# loads package. Note that it first must be downloaded from the internet using the package manager in R. 
# This needs to be done only once.

# load the data
df <- read.table(file.choose(), header=T, sep="\t")

# convert the data to cross-table (make sure you have run the script - see file data frame to stab script)
xtab = stack.fnc(df)

# check data
xtab

# invert the data
xtab.t <- t(xtab)

PVClust <- pvclust(xtab.t, method.hclust= "ward", method.dist= "euclidean") 
# performs the cluster analysis using Euclidean distance matrix and the agglomeration method "average". 
# bootstraps the results of estimated p-values, This can take some time.

# Plots results 
plot(PVClust, frame= T) 
pvrect(PVClust, alpha=0.95)

# for a list of clustering methods and distance matrices type:
?dist 

##################################################################
## Hierarchal Cluster Analysis with a phylogenetic plot in ape ###
##################################################################

library (ape)
# loads package. Note that it first must be downloaded from the internet using the package manager in R. This needs to be done only once.
# load the data
df <- read.table(file.choose(), header=T, sep="\t")

# convert the data to cross-table (make sure you have run the script - see file data frame to stab script)
xtab = stack.fnc(df)

# check data
xtab

# converts frequency table to distance matrix
data.dist <- dist(xtab)

# perform the cluster analysis
PhyloClust <- nj(data.dist)
plot(PhyloClust, type = "u", frame= T, cex= .9)

###############################################
###  Cluster of MCA output in FactoMineR   ####
###############################################
# Run a cluster analysis on the results of the MCA

library(FactoMineR)

# load the data
df <- read.table(file.choose(), header=T, sep="\t")

# check the data
summary(df)

# Perform Multiple Correspondence Analysis
MCAFac = MCA(df, graph= F)
plot = HCPC(MCAFac)

###########################################################
###  Cluster with Agglomerative Coefficient in agnes   ####
###########################################################

In agnes we can create a silhouette plot (also called banner plot) and computer the agglomerative coefficient (also called silhouette coefficient)
library(cluster)

# load the data
df <- read.table(file.choose(), header=T, sep="\t")

# convert the data to cross-table (make sure you have run the script - see file data frame to stab script)
xtab = stack.fnc(df)

# check data
xtab

# converts frequency table to distance matrix
data.dist <- dist(xtab)

# Do the cluster analysis
AgClust = agnes(data.dist, method="single")

# Plot the cluster analysis
plot(AgClust)

###################################################
###  K-Means Cluster w/ cluster plot in clara  ####
###################################################

library(cluster)

# load the data
df <- read.table(file.choose(), header=T, sep="\t")

# perform the analysis, specify the number of clusters, here it's 3 
ClaraClust <- clara(df, 3, medoids.x = TRUE, metric = "manhattan", stand = T, samples = 500)

# plot the results
plot(ClaraClust)

# Interpreting coefficient scores
# < 0.25 no substantial structure found;
# 0.26-0.50 structure is found but it is weak; 
# 0.51-0.70 a reasonable structure identified; 
# > 0.71 a strong structure identified 
# Reference UNESCO 2013

plot(ClaraClust, color = T, cex = .8, col.clus = if(1) c(1, 3, 4, 6) else 5, shade = T)

 clusplot(ClaraClust, clus,  stand = T, lines = 1, shade = F, color = T, labels= 0, plotchar = TRUE, col.p = "dark green", col.txt = col.p,  col.clus = if(1) c(1, 3, 4, 6) else 5, cex = .8, cex.txt = cex, span = TRUE, add = FALSE, xlim = NULL, ylim = NULL,  main = paste("CLUSPLOT(", deparse(substitute(x)),")"),  sub = paste("These two components explain",  round( digits = 2), "% of the point variability."),  xlab = "Component 1", ylab = "Component 2",  verbose = getOption("verbose"))
################################################
###  Creating a “stacked” cross-tabulation  ####
################################################

# Step 1: Load the data frame
# Important! Make sure the first column in the data frame 
# is the column you want as row names 

df <- read.table(file.choose(), header=T, sep="\t")

# Step 2: Run this script 
# Copy / paste into R and enter. Only do once per session.

stack.fnc=function(df.flat){
 for(i in 1:ncol(df.flat))
 df.flat[,i]=factor(df.flat[,i])
 df.stacked.nrow=length(levels(df.flat[,1]))
 df.stacked.ncol=0
 for(i in 2:ncol(df.flat))
 df.stacked.ncol=df.stacked.ncol+length(levels(df.flat[,i]))
 k=1
 df.stacked.colnames=rep("X",df.stacked.ncol)
  for(i in 2:ncol(df.flat))
  for(j in 1:length(levels(df.flat[,i]))) {df.stacked.colnames[k]=paste(colnames(df.flat[i]),levels(df.flat[,i])[j],sep=".")
  k=k+1}
df.stacked=data.frame(matrix(nrow=df.stacked.nrow,ncol=df.stacked.ncol))
df.stacked=matrix(data=rep(0,df.stacked.nrow*df.stacked.ncol),nrow=df.stacked.nrow,ncol=df.stacked.ncol,dimnames=list(levels(df.flat[,1]),df.stacked.colnames))
  for(i in 2:ncol(df.flat))
  for(j in 1:nrow(df.flat))
  df.stacked[as.character(df.flat[j,1]),paste(colnames(df.flat)[i],df.flat[j,i], sep=".")]=
  df.stacked[as.character(df.flat[j,1]),paste(colnames(df.flat)[i],df.flat[j,i], sep=".")]+1
  rm(list=c("df.stacked.ncol","df.stacked.nrow","i","j","k"))
  return(data.frame(df.stacked))
  }

# Step 3: This creates your contingency table called "xtab"
xtab = stack.fnc(df)

# Step 4: Check it has worked
xtab