##########################################
##### Multivariate Ratio Analysis (MRA), see Baur H., Leuenberger C. (2011) Systematic Biology 60: 813-825
##### 2. Tutorial: The *LDA Ratio Extractor*
##########################################


# **** Make sure, source and data files are in one and the same folder! ****

rm(list=ls()); ls()  # Tidying up workspace
source("LDA_Ratio_Extractor_v1.01.R", chdir = TRUE) # Load newest verson of the script

#### Load Dataset
filename <- "Anisopteromalus_data.csv"
dat0 <- read.csv(filename)





#### Using the LDA Ratio Extractor

#### Example 1 ####

# Use all data
dat <- dat0
X <- dat[,3:21] # Select data variables. One variable, gst.b, is excluded, since it contains artifacts (see Baur et al. 2014, p. 695-696 & Fig. 1)
species <- factor(dat$species) # Select grouping variable, here species

RE <- ldaRE(U=X, g=species, 
             k.max=4,  # Defaul setting
             rx=1, ry=2,  # Defaul setting
             g1=1, g2=2,  # Defaul setting
             col=c(1:nlevels(species))[species], # Default setting
             pch=c(1:nlevels(species))[species], # Default setting
             placelegend="bottomright", # Default setting
             cexlegend=0.9 # Defaul setting
             )

RE # Print entire function output
RE$ratios # Print only data frame with ratios. May be used for further analysis




#### Example 2 ####

# Using a subset of the data with group colours and symbols adjusted according to entire dataset
dat <- subset(dat0, species == "calandrae" | species == "caryedophagus" | species == "quinarius")
X <- dat[,3:21]; species <- factor(dat$species); factor.levels <- c(2,3,6)

RE <- ldaRE(U=X, g=species, 
             k.max=2, 
             rx=1, ry=2, 
             g1=1, g2=3,  # Print statistics for calandrae versus quinarius
             col=factor.levels[species], 
             pch=factor.levels[species], 
             pchlegend=factor.levels, 
             collegend=factor.levels, 
             placelegend="topright", 
             cexlegend=1.1
             )



#### Example 3 ####

# For finding the very best separation, it is usually best to make pairwise comparisons
dat <- subset(dat0, species == "calandrae" | species == "quinarius")
X <- dat[,3:21]; species <- factor(dat$species); factor.levels <- c(2,6)

RE <- ldaRE(U=X, g=species, 
             k.max=2, 
             rx=1, ry=2, 
             g1=1, g2=2,  # Print statistics for calandrae versus quinarius
             col=factor.levels[species], 
             pch=factor.levels[species], 
             pchlegend=factor.levels, 
             collegend=factor.levels, 
             placelegend="bottomright", # Default setting
             cexlegend=1.1
             )

# Remark:
# Now, the very best ratio for separating calandrae from quinarius is:
# ool.l/eye.h!
# In the 1. Tutorial, Example 2, by using the PCA Ratio Spectrum we found 
# ool.l/eye.b as the most significant ratio for the two species.
# The ratio ool.l/eye.h was only number 4 in the spectrum.
# Can you figure out why this could be the case? 
# It has to do with a basic difference between
# PCA (Principal Component Analysis) and LDA (Linear Discriminant Analysis).




#### Example 4 ####

# Regrouping specimens for finding best separation with pairwise comparison
dat <- dat0
X <- dat[,3:21]; species <- factor(dat$species); levels(species) <- list("apiovorus" = c("apiovorus"), "rest" = c("calandrae", "caryedophagus", "ceylonensis", "cornis", "quinarius"))
factor.levels <- c(1,8)

RE <- ldaRE(U=X, g=species, 
             k.max=2, 
             rx=1, ry=2, 
             g1=1, g2=2,  # Print statistics for apiovorus versus rest
             col=factor.levels[species], 
             pch=factor.levels[species], 
             pchlegend=factor.levels, 
             collegend=factor.levels, 
             placelegend="bottomleft", 
             cexlegend=1.1
             )




#### Example 5 ####

# Regrouping specimens for finding best separation with pairwise comparison
dat <- dat0
X <- dat[,3:21]; species <- factor(dat$species); levels(species) <- list("ceylonensis" = c("ceylonensis"), "rest" = c("apiovorus", "calandrae", "caryedophagus", "cornis", "quinarius"))
factor.levels <- c(4,8)

RE <- ldaRE(U=X, g=species, 
             k.max=2, 
             rx=1, ry=2, 
             g1=1, g2=2,  # Print statistics for ceylonensis versus rest
             col=factor.levels[species], 
             pch=factor.levels[species], 
             pchlegend=factor.levels, 
             collegend=factor.levels, 
             placelegend="topleft", 
             cexlegend=1.1
             )



#### Example 6 ####

# Regrouping specimens for finding best separation with pairwise comparison
dat <- dat0
X <- dat[,3:21]; species <- factor(dat$species); levels(species) <- list("cornis" = c("cornis"), "rest" = c("apiovorus", "calandrae", "caryedophagus", "ceylonensis", "quinarius"))
factor.levels <- c(5,8)

RE <- ldaRE(U=X, g=species, 
             k.max=2, 
             rx=1, ry=2, 
             g1=1, g2=2,  # Print statistics for cornis versus rest
             col=factor.levels[species], 
             pch=factor.levels[species], 
             pchlegend=factor.levels, 
             collegend=factor.levels, 
             placelegend="topleft", 
             cexlegend=1.1
             )