####################################################
######   Exploratory Techniques R Commands   #######
######       Correspondence Analysis         #######
####################################################

## Correspondence analysis helps identify associations in complex data

## The packages needed for today’s tutorial are: 
MASS
ca
FactoMineR
explor

## Don’t forget you need to **install** packages just once, 
## but need to **load** them every time you restart R (if you use them)

## If we get time, we will briefly look at Loglinear Analysis as well.
For this, you will need the packages
 
library(vcd)
library(MASS)     ## this comes pre-installed, so there is no need to download it

## I don’t know whether to assume knowledge of Chi Squared Test here, but the commands for that are:

#########################################
###  Loading and Examining your data  ###
#########################################

## You do not need packages for this.

## Load data from data frame in text file format, comma delimited

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

## Check your data is in order
summary(df)

## If you have a lot of levels in a single factor
table(df$NameOfColumn)

## Have a look at the first few rows
head(df)

#################################
###  Working with subsamples  ###
#################################

## The functionalities in R here go beyond what we can cover. 
## Arguably the most useful and simple command is “subset”

Subset1 = subset(df, select = c(NameOfColumn1, NameOfColumn2, NameOfColumn3))

#################################################
###  Multiple Correspondence Analysis  in ca  ###
#################################################

# The ca library has a range of different algorithms for doing the analysis,
# some of these allow us to have numerical output for correspondence analysis

# load the MASS library /  package (this only needs to be done once per session)
library(ca)

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

# check the data
summary(df)

# Burt MCA
MCA_Burt <- mjca(df, lambda= "Burt")
plot(MCA_Burt, labels= c(0, 2), col= c("white", "black")) 
summary(MCA_Burt)

# Joint MCA
MCA_Joint <- mjca(df, lambda= "JCA")
plot(MCA_Joint, labels= c(0, 2), col= c("white", "black")) 
summary(MCA_Joint)

# Adjusted MCA
MCA_ADj <- mjca(df, lambda= "adjusted")
plot(MCA_ADj, labels= c(0, 2), col= c("white", "black")) 
summary(MCA_ADj)

##########################################################
###  Multiple Correspondence Analysis in FactoMineR   ####
##########################################################

# Good for very detailed numeric output

library(FactoMineR)
library(explor)

# 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)
explor(MCAFac)

###############################################
###  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)

# You will need to “cut” the tree at certain level.
# R will suggest a number of clusters to “cut” at with a dark horizontal line
# R also gives you a scree plot to help you decide

###################################################
###  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"))
##################################################
###              Chi Squared test              ###
### Checking for independence in two-way table ###
##################################################

# where file is two columns of a data frame
df <- read.table(file.choose(), header=T, sep="\t") 

# check data is in order
summary(df)

# make a cross tabulation of your data
xtab = table(df)

# check data is in order
xtab

# turn the data around if you want to
xtab.t = t(xtab)

# run the test for Chi Squared test for independence
CHI = chisq.test(xtab)
CHI

# examine the residuals
RES= CHI$res
RES

# for explanatory purposes, comparing observed and expected
O = CHI$observed
E = CHI$expected
E
O
(O-E)^2/E

############################################################
###            LogLinear Analysis                        ###
### Checking for significant correlations in k-way table ###
############################################################

# Libraries Needed
library(vcd)
library(MASS)

# Load and examine data from data frame 
# use subset command if needed - selecting only the variables you wish to examine. 
# Note that the order of the columns will determine the order they are treated by the analysis. You may have re-important data from external file.

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

# convert dataframe to table

d = table(df)

# Run loglinear analysis (effective maximum of 4 variables)
LogLin = loglm(formula = ~ V1 + V2 +  V3 , data = d)

# Plot the results in a mosaic plot
mosaic(LogLin, legend_args = list(fontsize = 8 ), labeling_args = list(
gp_labels = gpar(fontsize = 8, fontface = 1,  fontfamily = "Times New Roman"),
gp_varnames = gpar(fontsize = 10, fontface = 2, fontfamily = "Times New Roman")))