<table width = 900 border = 0><tr><td>
# <b>NMST539 | Lab Session 12</b>

## <b>Discriminant Analysis & Classification</b>
### <font color = "#880000"><b>Summer Term 2021/2022 |  03/05/22 </b></font>
###### <a href = "Lab12_2022.Rmd">source file (UTF8 encoding)</a>

**************


The R software is available for download (free of charge) from the website: <a href = "https://www.r-project.org">https://www.r-project.org</a>

A user-friendly interface (one of many): <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Manuals and introduction into <i>R</i>
(in Czech or English):</b><br>
    
<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i> (in Czech) | (<a href = "NMFM301/Rmanual2.pdf">PDF manual</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i> (in Czech) | (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF manual</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (in Czech) | (<a href = "NMFM301/uvodr.pdf" >PDF manual</a>)</li>
<li>Venables, W.N., Smith, D.M., and the R Core Team: <i>An Introduction to R.</i> |  (<a href = "https://cran.r-project.org/doc/manuals/R-intro.pdf">PDF manual</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) |   (<a href = "http://www.ievbras.ru/ecostat/Kiril/R/Biblio/R_eng/R%20dummies.pdf">PDF book</a>)</li>		
</ul>

<br><br>

<b>Theoretical references</b><br>
<ul>
<li>Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015</li>
<li>Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.</li>
</ul>

**************


### <font color = "#880000"><b>1. Classification via discriminant analysis</b></font>

Discriminant analysis is (not necessarily a multivariate) statistical method  designed for allocating various objects into similar (homogenous) groups. The idea is to use some available information (given covariates) in order to create some classification rule which can be used later to assign  new objects into one of the existing groups -- so called subject classification. Geometrically speaking, the discriminant analysis can be seen as some specific partition of the whole sample space -- usually $\mathbb{R}^p$.
Recall the main difference between the clustering algorithms and the classification approaches:<br><br> 

 <ul>
 <li><b>Classification</b> is applied on a dataset in which the groups of subjects (observations) -- so called <b>labels</b> -- are already known. The classification algorithm 
 uses the data to learn  some principle in order to differentiate new observations into the given groups. </li>
 
 <li><b>Clustering</b> algorithms use no prior information regarding the underlying groups (<b>no labels</b>). The groups are unknown and the algorithm tries to 
 identify the groups by considering mutual similarities/differencies in the given observations.</li>
 </ul>

As mentioned above, the classification algorithm tries to learn some underlying principle -- a <b>discimination rule</b> -- to distinguish heterogeneous observations (which should be ideally classified in different groups) and to classify same observations which are rather similar. 
There are of course many different approaches  to define some reasonable descrimination rule. Such rules can be based, for instance, on the prior probabilities (Bayes rule), the likelihood approach, or more complex (sophisticated neural network) algorithms.  

Using a linear regression jargon, the discriminant analysis can be also seen as an alternative/counterpart for a linear regression model (continuous dependent variable and continuous or discrete predictor variables)  used mainly to analyze the  data when the criterion of interest (i.e., dependent variable) is categorical  and the predictors are continuous (or discrete) covariates. The term categorical variable means that the dependent variable is divided into a specific number of different and disjoint categories.

<br><br>


In practical applications, the main focuss is given on a <b>linear discriminant analysis (LDA)</b>  (a discrimination rule based on a maximum likelihood approach for two multivariate normal populations with the same variance-covariance matrix) which is a statistical method meant to find a linear combination of features (usually a set of continuous random variables) that characterizes or separates two (or more) classes of objects. The resulting linear combination is used as a linear classifier. 
The linear discriminant analysis is also closely related to the analysis of variance (ANOVA) and the regression analysis, which both attempt to express one dependent variable as a linear combination of other features or measurements. However, unlinke ANOVA, which uses categorical independent variables and a continuous dependent variable, the discriminant analysis works other-way around -- it uses continuous independent variables and a categorical dependent variable instead. 

Alternatively, the logistic regression approach is more similar to LDA than ANOVA as it also explains a categorical variable by using a set of continuous (discrete) independent variables. However,  GLM methods are preferred  in applications where it is not reasonable to assume that the data variables are normally distributed, which is a fundamental assumption of the LDA method. A straightforward extention of LDA is a quadratic discrimination rule -- a discrimination rule again based on a maximum likelihood approach for two multivariate  normal populations, however, with different variance-covariance matrices). 

If the normality assumptions is too much restrictive, the Fisher's approach is usually adopted instead. It finds the linear discrimination rule in a way that the <b>between class variability</b> and the <b>within class variability</b> are optimized. This method does not even require homoscedastic errors. 

Suppose that two groups of observations have their means denoted as $\boldsymbol{\mu}_1$ and $\boldsymbol{\mu_2}$ and the corresponding variance-covariances matrices are $\Sigma_1$ and $\Sigma_2$. Thus, there are  observations $\boldsymbol{X}_1, \dots, \boldsymbol{X}_n \sim (\boldsymbol{\mu}_1, \Sigma_1)$ from the first group (only the mean and variance-covariance are specified) and other observations $\boldsymbol{Y}_1, \dots, \boldsymbol{Y}_m \sim (\boldsymbol{\mu}_2, \Sigma_2)$ (the dimension observations from both groups is the same, thus $\boldsymbol{\mu}_1, \boldsymbol{\mu}_2 \in \mathbb{R}^p$). 

The Fisher's linear separation measure between these two groups (the corresponding distributions respectively) is defined as a ratio of the variance between the groups to the variance within the groupds, which is 
$$
S = \frac{\sigma_{between}^2}{\sigma_{within}^2} = \frac{(\boldsymbol{w}^\top \boldsymbol{\mu}_2 - \boldsymbol{w}^\top \boldsymbol{\mu}_1)^2}{\boldsymbol{w}^\top (\Sigma_{1} + \Sigma_2) \boldsymbol{w} }.
$$

The idea is to find an appropriate vector of a linear combination $\boldsymbol{w} \in \mathbb{R}^p$ such that the expression above is maximized. The quantity measures, in some sense,  the signal-to-noise ratio for the class labeling and it can be shown that the maximum separation occurs for 
$$
\boldsymbol{w} = (\Sigma_1 + \Sigma_2)^{-1} (\boldsymbol{\mu}_2 - \boldsymbol{\mu}_1).
$$

<br><br>

### <font color = "#880000"><b>2. Discriminant analysis in R</b></font>

The linear discriminant analysis (LDA) is implemented in the R software in the function `lda()` which is available in the packages 'MASS'.
Similarly, the quadratic discriminant analysis (QDA) is implemented by the command  `qda()` (also in the package 'MASS').

In the following, let us again consider the `mtcars` data set. We use the subset of the continuous covariates only and we search for a linear classification rule (a specific partition of the sample space $\mathbb{R}^5$) to classify the given observations (cars) by the number of cylinders (thus, we have three groups -- cars with 4 cylinders, 6 cylinders, and 8 cylinders). 

Considering a probabilistic model, there is a random vector $\boldsymbol{X} = (mpg, hp, qsec, drat, disp)^\top \in \mathbb{R}^5$ and the random sample $\boldsymbol{X}_1, \dots, \boldsymbol{X}_{32}$ drawn from the distribution of the generic random vector $\boldsymbol{X}$. 

```{r, fig.width = 10}
library(MASS)
lda1 <- lda(cyl ~ mpg + hp + qsec + drat + disp, data = mtcars)

pred <- predict(lda1)
ldahist(data = pred$x[,1], g=mtcars$cyl)
```

From the geometrical point of view the linear discriminant analysis is a linear partition of the $5$ dimensional sample space. 
A graphical illustration of the particion can be obtained by using  the R function `partimat()` from the library `klaR`. 

```{r, fig.width=10, fig.height=8}
library("klaR")
partimat(as.factor(cyl) ~ mpg + hp + qsec + drat + disp, data = mtcars, method = "lda")
```

In the figure above, the partition of the sample space is provided with respect to all possible pairs of two-dimensional (xy) scatterplots. 
Alternatively, we can visualize the same data (without the sample space partition) by using the standard command `pairs()`. 

```{r, fig.width=10, fig.height=8}
pairs(mtcars[c("mpg", "hp", "qsec", "drat", "disp")], pch = 22, bg = c("red", "yellow", "green")[unclass(as.factor(mtcars$cyl))])
```

There are three groups, thus, we have two directions maximazing the between and withing variablity. 
It can be all vizualised as follows: 

```{r, fig.width=10}
tData <- cbind(as.matrix(mtcars[, c("mpg", "hp", "qsec", "drat", "disp")]) %*% lda1$scaling, mtcars$cyl)
m4 <- apply(tData[tData[,3] == 4, 1:2], 2, mean)
m6 <- apply(tData[tData[,3] == 6, 1:2], 2, mean)
m8 <- apply(tData[tData[,3] == 8, 1:2], 2, mean)
plot(tData[,2] ~ tData[,1], pch = 22, bg = c("red", "yellow", "green")[unclass(as.factor(mtcars$cyl))], xlab = "First LDA direction", ylab = "Second LDA direction")
points(m4[1], m4[2], cex = 2, pch = 21, bg = "red")
points(m6[1], m6[2], cex = 2, pch = 21, bg = "yellow")
points(m8[1], m8[2], cex = 2, pch = 21, bg = "green")
```

<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

#### <font color = "#880000"><b>Note</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>For a better graphical insight into the LDA results we can consider a simplification where we only use one continuous covariate as a predictor. 
```{r}
lda2 <- lda(cyl ~ mpg, data = mtcars)
lda2cv <- lda(cyl ~ mpg, data = mtcars, CV = T)
```
Thus, the linear discriminant rule is given in terms of one scaling factor only. The corresponding transformations can be obtained as follows: 
```{r, fig.width = 10}
lda_scores <- as.numeric(lda2$scaling) * mtcars$mpg
par(mfrow = c(1,2))
plot(lda_scores ~ mtcars$mpg, pch = 21, bg = mtcars$cyl, ylab ="LDA Scores", xlab = "Miles per Gallon")

lda2cv <- cbind(lda2cv$posterior, mtcars$mpg, mtcars$cyl)
lda2cv <- lda2cv[order(mtcars$mpg), ]

plot(lda2cv[,1] ~ lda2cv[,4], pch = 21, bg = lda2cv[,5], ylab = "Posterior Probability", xlab = "Miles Per Gallon")
lines(lda2cv[,1] ~ lda2cv[,4], col = "blue")
points(lda2cv[,2] ~ lda2cv[,4], pch = 22, bg = lda2cv[,5])
lines(lda2cv[,2] ~ lda2cv[,4], col = "violet")
points(lda2cv[,3] ~ lda2cv[,4], pch = 23, bg = lda2cv[,5])
lines(lda2cv[,3] ~ lda2cv[,4], col = "gray")

for (i in 1:nrow(lda2cv)){
  lines(rep(lda2cv[i,4], 2), c(0,1), col = "gray", lty = 3)
}
```
 </li>
 
 <li>Consider an analogous situation for more than just one continuous predictor. Use some graphical tools to visualize the results.</li>
 
 <li>Use some additional datapoints (feel free to make up some) and try to use the discrimination 
 rules to assign the new points to the existing groups. </li>
 
 <li>Consider the heteroscedastic version of the model and try to apply the quadratic discriminant function (command 'qda()' in package 'MASS').</li>
 
 <li>Consider, for instance, an alternative classification perfomed within the framework of a logistic model.</li>
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

#### <font color = "#880000"><b>Mowers data: Comparison of various methods</b></font>

Let us consider an artificial and balanced dataset where only three covariates are given for each subject (house holder):  the income of the house holder; the lot landscape; the type of the grass mower. Using the information about the income and the landscape we want to find a classification rule for the mower type. Four different methods are used: generalized linear model; linear discriminant analysis; quadratic discriminant analysis; linear regression model;

The data can be loaded in by the command

```{r}
Mowers <- read.table("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/mowers.txt", sep = "", header = T)

### 2D grid expansion
mowers.gr = expand.grid(seq(20,120,len=1000),seq(14,24,len=1000))
names(mowers.gr) = c("income","lot")
```

And the results of the four methods are the following: 

```{r, fig.width = 10, fig.height = 10}
par(mfrow=c(2,2))
cols <- c("red", "blue")
attach(Mowers)
mowe.ii = unique(mowers.gr$income)
mowe.ll = unique(mowers.gr$lot)

plot(income,lot,pch=21 + riding, bg = cols[riding], main = "Generalized Linear Model")
mow.nn = glm(factor(riding)~income+lot,family=binomial,data=Mowers)

contour(mowe.ii,mowe.ll,
        matrix(as.numeric(predict(mow.nn,new=mowers.gr,type="response")>0.5),1000,1000),
        nlev=1,drawlabels=FALSE,add=T)

plot(income,lot,pch=21 + riding, bg = cols[riding], main = "Linear Discriminant Analysis")
mow.nn = lda(factor(riding)~income+lot,data=Mowers)
contour(mowe.ii,mowe.ll,
        matrix(as.numeric(predict(mow.nn,new=mowers.gr)$class),1000,1000),
        nlev=1,drawlabels=FALSE,add=T)

plot(income,lot,pch=21 + riding, bg = cols[riding], main = "Quadratic Discriminant Analysis")
mow.nn = qda(factor(riding)~income+lot,data=Mowers)
contour(mowe.ii,mowe.ll,
        matrix(as.numeric(predict(mow.nn,new=mowers.gr)$class),1000,1000),
        nlev=1,drawlabels=FALSE,add=T)

plot(income,lot,pch=21 + riding, bg = cols[riding], main = "Linear Regression Model")
mow.nn = lm(I(2*(riding-1)-1)~income+lot,data=Mowers)
contour(mowe.ii,mowe.ll,
        matrix(as.numeric(predict(mow.nn,new=mowers.gr,type="response")>0),1000,1000),
        nlev=1,drawlabels=FALSE,add=T)
```

Note, that the LDA approach and the standard linear regression performs equivalently (under the assumption that there are two groups only, the groups are labeled 
with plus/minus one and, finaly, the data are balanced -- the same number of observations in each group -- 12 in this particular example). <br><br>

Another discrimination approach gaining more and more popularirty is based on deep neural networks which were originally designed to perform such and similar tasks efficiently. 
The deep neural network is a highly non-linear non-convex optimization method. Unfortunately, the results heavily depends on an initial configuration: compare, for instance, the following four outputs (the same neural network is applied for the same dataset for all four cases):

To run the command, one needs to install the corresponding package `install.library(nnet)`. 

```{r, fig.width = 10, fig.height = 10}
require(nnet)
par(mfrow=c(2,2))

for (k in 1:4) {
  mow.nn = nnet(factor(riding)~income+lot,size=50,data=Mowers)
  plot(income,lot,pch=21 + riding, bg = cols[riding], main = "Generalized Linear Model")
  contour(mowe.ii,mowe.ll,
          matrix(as.numeric(predict(mow.nn,new=mowers.gr,type="class")),1000,1000),
          nlev=1,drawlabels=FALSE,add=T)
}
```


In order to get some meaningful results one usually runs the deep neural network for a few times and the final result is averaged over all runs. <br><br>


<br><br>

Other very popularclassification methods (however, not that must transparent from the mathematical point of view) are, for instance, <b>decision trees</b> (command `ctree()` in 
the R library 'install.packages("party")'), <b>support vector machines</b> (command `ksvm()` in the R library 'install.packages("kernlab")'), or the so called <b>K – nearest neighbor</b> approach (command `knn()` in the R library 'install.packages("class")'). 


<br><br>

<b>Decision trees</b>

```{r, fig.width = 10,  message = FALSE, warnings = FALSE}
library("party")

plot(ctree(factor(riding)~income+lot,data=Mowers))
```

<br><br>

<b>Decision trees</b>

```{r, fig.width = 10,  message = FALSE, warnings = FALSE}
library("class")

train <- mtcars[1:25,c("mpg", "hp", "qsec", "drat", "disp")]
test <- mtcars[26:32,c("mpg", "hp", "qsec", "drat", "disp")]

c1 <- as.factor(mtcars$cyl[1:25])

knn(train, test, c1, k = 3, prob=TRUE)
```



<br><br>



## <font color = "#880000"><b>Homework Assignment</b></font>
### (Deadline: Lab Session 13 / 10.05.2022)

<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<ul>
<li>Use some appropriate data set (ideally the same you used for your previous homework assignments) and apply 
the linear discriminat analysis classification. Consider two groups -- either given within some covariate which is already given 
in the data or define some arbitrary binary covariate based on some preliminary exploratory analysis (performed in 
the previous homework assignments). </li>

<li>Use at least one additional alternative classification approach (e.g., linear regression, glm, neural networks, or qda) 
and compare the results. Provide some graphical outputs and some brief interpretation/explanation for both results. </li>

<li>Add your solution to your personal web page, on <font color = "red">Monday, 08.05.2022, 23:59 (CET)</font> at latest. 
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

<br><br>

</td></tr></table>
<br><br><br>