###################################################
######     Handling Data and Significance      ####
######               Poznan 2022               ####
###################################################

#####################################
######       Loading Data        ####
#####################################

# Step one - Copy data from spreadsheet into text editor (not Word Pad)
# Step two - Load data into R from txt file
df <- read.table(file.choose(), header=T, sep="\t", stringsAsFactors = T)

# First rows of data
head(df)

# Summary of dataframe
# NB: some recent versions of R, this command is not working
summary(df)

# Summary of one 'column' of a dataframe
table(df$NameOfColumn)

##############################
######       Chi2         ####
##############################

# Step one - Copy data from spreadsheet into text editor (not Word Pad)
# Step two - Load data into R from txt file
df <- read.table(file.choose(), header=T, sep="\t", stringsAsFactors = T)

xtab <- table(df)

# check data, to make sure it has worked
xtab

CHI = chisq.test(xtab)
CHI

# Pearson residuals
# These show which cells are responsible for differences between the sets
RES= CHI$res
RES

# To understand how it works
# Observed vs. Expected outcomes
O = CHI$observed
E = CHI$expected
E
O

######################################
######     Association Plots      ####
######################################

library(vcd)
library(MASS)

# A simple plot that can help visualise differences in frequency 
# in binomial data 

# Load data in to R
df <- read.table(file.choose(), header=T, sep="\t", stringsAsFactors = T)

# Convert data to cross-tabulation
xtab <- table(df)

# Plot the data
assoc(xtab)

# Sometimes it is clearer to transpose (turn around) the table
xtab2=t(xtab)

# Plot the transposed data
assoc(xtab2)

# Examples of commands to make the plot prettier
assoc(xtab2, compress=T, spacing_equal(sp = unit(0.8, "lines")), residuals_type = "pearson", gp_varnames = gpar(fontsize = 10, fontfamily = "Times New Roman"), gp_labels = gpar(fontsize = 8, fontfamily = "Times New Roman", keep_aspect_ratio = F),  xscale = 0.9, yspace = unit(0.4, "lines"))

#######################################
######       LogLinear Analysis  	  ####
#######################################

# Libraries Needed
library(vcd)
library(MASS)

# Load and examine data from data frame - 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.

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

# convert dataframe to table
d = table(df)

# Run loglinear analysis (effective maximum of 4 variables)
# Replace V1, V2 etc with the names of the variables
LL = loglm(formula = ~ V1 + V2 +  V3 , data = d)

# Plot the results in a mosaic plot
mosaic(LL)

# Examples of commands to make the plot prettier
mosaic(LL, 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")))

#####################
### Handling data ###
#####################

# Change the order of features in a variable 
df $ ColumnName = relevel(df$ColumnName, ref= "FeatureName")

# Remove Columns
df $ NameOfColumn=NULL

# convert the data to cross-tabulation in R 
# make sure you have run the script at the end of this file
xtab = stack.fnc(df)

# check data, to make sure it has worked
xtab

# Creating a contingency table from
# For Cluster Analysis and Binary Correspondence Analysis

xtab = stack.fnc(df)

# Turn contingency table around (transpose)
xtab.t = t(xtab)
xtab.t

# Telling R that a given variable is a factor 
# (non numeric variable, often important for dates)
df $ ColumnName = as.factor(df$ColumnName)

# Telling R that a given variable is ordinal, not nominal 
# (essential for ordinal regression)
df $ ColumnName = as.ordered(df$ColumnName)

# Telling R that a given variable is numeric not categorical 
# (not sure why you would need this :)
df $ ColumnName = as.numeric(df$ColumnName)

# Telling R not to ignore NA
df $ ColumnName = addNA(df$ColumnName)

#################################################################
### Converting flat data frame to a stackec contingency table ###
#################################################################

# Step 1: Load the data frame
# Important! Make sure the first column in the data frame 
# (in excel) 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