## LDA/QDA examples

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

# 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 there is no provision for costs

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

## This function implements the RELIEF feature selection algorithm
# Usage

# relief(data, nosample, threshold,vnom)

# Arguments
# data 	the dataset for which feature selection will be carried out
# nosample 	number of instances drawn from the original dataset
# threshold 	the cutoff point to select the features
# vnom 	a vector containing the indexes of the nominal features

library (dprep)

data (my.iris)
relief (my.iris,150,-Inf) 

# 150 does not need to be the length of the data, it can be higher
# this will entail more resampling (at a higher computational cost)
# note the target variable has to be a class number, e.g. 1,2,3,...