## LDA/QDA examples of use

## First we separate test data

data(iris3)
Data <- data.frame(rbind(iris3[,,1], iris3[,,2], iris3[,,3]),
                   Sp = rep(c("s","c","v"), rep(50,3)))
Data
learn <- sample(1:150, 75)
learn
table(Data$Sp[learn])

mydata.learn <- Data[learn, ]
mydata.test <- Data[-learn, ]

## We would like first to have a look at the data in 2D:
## These are scatterplots with color coding by class

pairs(Data[c("Sepal.L.","Sepal.W.","Petal.L.","Petal.W.")], main="A nice Title ", pch=22,
   bg=c("red", "yellow", "blue")[unclass(Data$Sp)]) 

## Linear discriminant analysis 

library(MASS) 

## Note the priors can be changed at will; if unspecified, ML estimates
## from the supplied data will be used

mydata.lda <- lda(Sp ~ ., Data, prior = c(1,1,1)/3, subset = learn)

## or alternatively

## mydata.lda <- lda(Sp ~ ., prior = c(1,1,1)/3, data = mydata.learn)

## lda() prints discriminant functions based on centered (not standardized) variables. 
## The "proportion of trace" that is printed is the proportion of between-class variance 
## that is explained by successive discriminant functions.
## the list element 'class' gives the predicted class 
## The list element 'posterior' holds posterior probabilities  

mydata.lda 

## resubstitution error (optimistic!)

predict(mydata.lda, mydata.learn)$class

# or 

summary(predict(mydata.lda, mydata.learn)$class)

## may also have a look at

predict(mydata.lda, mydata.learn)$posterior

## Scatter plot using the first two discriminant dimensions
plot(mydata.lda)

## The partimat() function in the 'klaR' package can display the results 
## of linear or quadratic classification using two variables at a time.

library(klaR)
partimat(Sp ~ ., Data,method="lda") 

## Estimation of true error
## The test error is a much better (unbiased) estimation

lda.predictions <- predict(mydata.lda, mydata.test)
lda.predictions$class

## Store posterior probabilities of predicted classes in a vector

posteriors <- vector(length = dim(mydata.test)[1])

for (i in 1:length(posteriors)) 
{
   posteriors[i] = lda.predictions$posterior[i,lda.predictions$class[i]]
}

## Assessment of predictive accuracy via LOOCV (recommended)
# If CV=TRUE, the method returns results (classes and posterior probabilities) for leave-out-out crossvalidation.
# Note that if the prior is estimated, the proportions in the whole dataset are used.

mydata.lda.cv <- lda(Sp ~ ., prior = c(1,1,1)/3, data = Data, CV=TRUE) 
summary(mydata.lda.cv$class)

tab <- table(Data$Sp, mydata.lda.cv$class)  
sum(tab[row(tab)==col(tab)])/sum(tab)

# If a formula is given as the principal argument the object may be modified using update()
# Example: omit Petal.W from the model and re-evaluate

mydata.lda.cv.noPetal.W <- update(mydata.lda.cv, . ~ . - Petal.W.)

tab <- table(Data$Sp, mydata.lda.cv.noPetal.W$class)  
sum(tab[row(tab)==col(tab)])/sum(tab) 

## NOTE1: na.omit() omits any rows that have one or more missing values 
# usage: data.nomissing = na.omit(data)
# alternatively one can use the method 'ce.impute' in package 'dprep'

## NOTE2: sadly (and apparently) there is no provision for costs

## Quadratic Discriminant Analysis is the same, replacing 'lda' by 'qda'


##########################################################

# Laboratory EXERCISES

##########################################################

# EXERCISE 1: Design a validation protocol. Is the result different every time? Should this be 'fixed'? Why? Do we need a separate test set?

# EXERCISE 2: Design a search method.
# 	a. Forward search
#     b. Backward search
#     c. Hybrid search

# EXERCISE 3: work with a real-world data set: n cannot be increased, the knowledge on variable importance is (completely?) unknown. Devise a wrapper LDA/QDA-based feature selection solution. Compare LDA with QDA and all the search methods you have implemented. In conclusion, what is the best combination?

# suggestion: use a dataset on fragments of glass collected in forensic work.
# This data frame contains the following columns:
# RI refractive index; more precisely the refractive index is 1.518xxxx.
# The next 8 measurements are percentages by weight of oxides.
# Na sodium.
# Mg manganese.
# Al aluminium.
# Si silicon.
# K potassium.
# Ca calcium.
# Ba barium.
# Fe iron.
# type The fragments were originally classed into seven types, one of which was absent in this
# dataset. The categories which occur are window float glass (WinF: 70), window non-float
# glass (WinNF: 76), vehicle window glass (Veh: 17), containers (Con: 13), tableware (Tabl:
# 9) and vehicle headlamps (Head: 29).

data (fgl)

# warning! 6 classes for only 214 observations overall! Difficult modelling problem

# Starting hint:

glass.loocv.model <- lda (x=fgl[,-10],grouping=fgl$type,CV=TRUE)
summary(glass.loocv.model$class)
tab <- table(fgl$type,glass.loocv.model$class)
(glass.loocv.error <- 1-sum(tab[row(tab)==col(tab)])/sum(tab))