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

## <b>Discriminant Analysis & Canonical Corelation Analysis</b>
### <b><font color = "#880000">LS 2017/2018 | Monday 14/05/18</font></b>
###### <a href = "cvicenie2018_11.Rmd">Rmd file (UTF8 coding)</a>

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


The R-software is available for download from the official 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>Some useful manuals and tutorials for <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>
**************


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

Discriminant analysis is a (not purely a multivariate) statistical method which was 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 some new object into one of the existing groups. From this point of view it can be considered as a classification tool. To remind the difference between clustering and clussification, let us recall the following:<br><br> 

 <ul>
 <li><b>Classification</b> works with  data for which the groups are known and we try to learn 
 what differentiates them in order to be able to assign future observations into the given groups. </li>
 
 <li><b>Clustering</b> has no prior information on what are the underlying groups -- the groups are unknown and we try to identify the groups
 by considering mutual similarities/differencies in the given observations.</li>
 </ul>
 
There are many different approaches how to define some reasonable descrimination rule. Such rules can be based, for instance, on the prior probabilities (Bayes rule) and the likelihood approach. For a normal model and homoscedastic error the discrimination rule simplifies to a linear form (linear discriminant analysis) while for the heteroscedastic errors the discrimination rule has a quadratic form instead.  We mostly focuss on a linear discriminant analysis (LDA)  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 or events. The resulting combination may be used as a linear classifier. This method is implemented in the R software implemented in the function `lda()` which is available in the packages 'MASS'. An straightforward extention is a quadratic discrimination rule wich is performed with the R command 'qda()' (also in package 'MASS').

Linear discriminant analysis is closely related to analysis of variance (ANOVA) and regression analysis, which also attempt to express one dependent variable as a linear combination of other features or measurements. However, ANOVA uses categorical independent variables and a continuous dependent variable, whereas discriminant analysis has continuous independent variables and a categorical dependent variable instead. On the other hand, the logistic regression approach is more similar to LDA than ANOVA is, as they also explain a categorical variable by the values of continuous independent variables. These other methods are preferable in applications where it is not reasonable to assume that the independent variables are normally distributed, which is a fundamental assumption of the LDA method.

If the normality assumptions is too much restrictive, we can use Fisher's approach which finds the linear discrimination in a way that between class and within class variability 
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$. The Fisher's linear separation measure between for 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} \boldsymbol{\mu}_2 - \boldsymbol{w} \boldsymbol{\mu}_1)^2}{\boldsymbol{w}^\top (\Sigma_{1} + \Sigma_2) \boldsymbol{w} }.
$$

This measure is, in some sense, a measure of 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>

In the following, let us consider the `mtcars` data set. We use the subset of the continuous covariates and we search for a linear classification rule to classify the given cars by the number of cylinders (thus, we have three groups -- cars with 4 cylinders, 6 cylinders, and 8 cylinders).<br><br>

```{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 theoretical point of view the linear discriminant analysis is a linear partition of the $n$ dimensional sample space. In the example above we have 
five dimensional sample space $R^5$ with the random vector 
$$
\boldsymbol{X} = (X_1, \dots, X_5)^\top = (mpg, hp, qsec, drat, disp)^\top.
$$

The linear discriminant analysis provides a specific partition, which can be visualized 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")
```

which can be further compared with another plot where we only distinguish for different groups

```{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))])
```


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

#### <font color = "#880000">To Do</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)
 ```
The linear discriminant rule is given by the scaling factor. 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")
```
 </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>Consider, for instance, an alternative classification perfomed within the framework of a logistic model.</li>

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







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

LDA works when the measurements made on independent variables for each observation are continuous quantities. When dealing with categorical independent variables, the equivalent technique is discriminant correspondence analysis.<br><br>



### <font color = "#880000"><b>3. Canonical Correlation Analysis in R</b></font>

The canonical correlation analysis (CCA) is a multidimensional exploratory statistical method
which (roughly) operates on the same principle as the principal component analysis (PCA). However, unlike the PCA where one tries to explore the common variability within the given data set, the canonical correlation is used to explore the mutual 
variability (correlations) between two different sets of quantitative variables observed on the same experimental units. Thus, the PCA method deals with one dataset only and it tries to to reduce the overall dimensionality of the dataset using some linear combination of the initial variables with the maximum variace possible
while the CCA approach searhes for some linear combinations of the initial variables in order to achieve maximum correlations instead. <br><br>

For instance, for two multivariate random vectors $\boldsymbol{X}$ and $\boldsymbol{Y}$ we assume that the common variance-covariance matrix of $(\boldsymbol{X}, \boldsymbol{Y})^\top$ is of the form 
$$
Var \left(\begin{array}{c}\boldsymbol{X}\\ \boldsymbol{Y}\end{array}\right) = \left(\begin{array}{cc} \Sigma_{XX} & \Sigma_{XY}\\ \Sigma_{YX} & \Sigma_{YY} \end{array}\right).
$$

If $\boldsymbol{X}$ is a $p$ dimensional random vector and $\boldsymbol{Y}$ is a $q$ dimensional random vecton, then the corresponding correlation of any linear combination of $\boldsymbol{X}$ and $\boldsymbol{Y}$ in terms
of $\boldsymbol{a}^\top \boldsymbol{X}$ and $\boldsymbol{b}^\top\boldsymbol{Y}$, for $\boldsymbol{a} \in \mathbb{R}^p$ and $\boldsymbol{b} \in \mathbb{R}^q$ can be expressed as 
$$
Cor(\boldsymbol{a}^\top \boldsymbol{X}, \boldsymbol{b}^\top\boldsymbol{Y}) = \boldsymbol{a}^\top \Sigma_{XY} \boldsymbol{b}.
$$

Therefore, the canonical correlations analysis searches for such vectors $\boldsymbol{a} \in \mathbb{R}^p$ and $\boldsymbol{b} \in \mathbb{R}^q$, that the previous correlation is as large as possible. However, 
the canonical correlation analysis does not search for just one pair of the canonical variables $\boldsymbol{a}^\top \boldsymbol{X}$ and $\boldsymbol{b}^\top\boldsymbol{Y}$ but instead, it defines all pairs (up to $min(p, q)$), 
such that all canonical covariates within each dataset are uncorrelated among each other. <br><br>

In the R statistical software there is function `cancor()` available under the standard instalation. 
Some additional packages can be downloaded and installed as well (mainly the R packages 'CCA' and 'vegan'). 

In the following we will consider the same dataset we have already worked with before. The data consists 
of two datasets where each dataset represents measurements in some specific (64-65) river localities in the Czech republic. The first dataset contains measurements of biologial metrics for each locality (17 different metrics and taxons) and the secon dataset contains measurements on chemical concentrations and values at the same localities (7 covariates recored). 

The idea is to somehow relate both dataset in order to explain what bio metrics can correlate which which chemical concentrations. One way of answering this is to apply the canonical correlation approach explained below. 


```{r, warning = F, message = F}
rm(list = ls())
bioData <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/bioData.csv", header = T)
chemData <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/chemData.csv", header = T)

head(bioData)
head(chemData)
```
<br><br>

Again, there is one locality which is missing the dataset of chemical measurements. We identify this locality and we do not consider it for further analysis. Both datasets are alligned with respect to he locality names - the first column in each dataset. 


```{r, warning = F, message = F}
ind <- match(chemData[,1], bioData[,1])
data <- data.frame(bioData[ind, ], chemData[, 2:8])
```


```{r}
X <- data[,2:9]
Y <- data[,19:25]
```

A nice graphical tool to view both datasets with respect to their overall (within and between) correlation structure is avialable within the R package 'CCA' (Canonical Correlation Analysis). The packages needs to be firstly installed (`install.packages("CCA")`) and later one can use functions `matcor()` and `img.matcor()` to graphically display the corration matrix 

$$ cor(\mathcal{X}, \mathcal{Y}) = \left(\begin{align}\Sigma_{XX} & \Sigma_{XY}\\\Sigma_{YX} & \Sigma_{YY}\end{align}\right)$$.

```{r, warning=F, message = F, fig.width=10, fig.height=8}
library("CCA")
correl <- matcor(X, Y )
img.matcor(correl, type = 2)
```


Next we can directly use `cancor()` function to obtain canonical correlations for the datasets $\mathcal{X}$ and $\mathcal{Y}$. Another option is to use function `cc()` (from the R package 'CCA'). Note, the the results are equivalent although not identical. 

```{r}
cc1 <- cancor(X, Y)  ### function from standard R instalation
cc2 <- cc(X, Y)      ### function for the R package 'CCA'
```

The canonical correlations are the same for both approaches (R commands used) 

```{r}
cc1$cor  ### function from standard R instalation
cc2$cor  ### function for the R package 'CCA'
```


The canonical covariates can be easily reconstructred as follows: 
```{r}
CanX <- as.matrix(X) %*% cc1$xcoef
CanY <- as.matrix(Y) %*% cc1$ycoef

cor(CanX, CanY)
```

Which corresponds with the output provided by the function `cancor()` and `cc()`above (see also the picture below
for a visual comparison).<br><br>



```{r, warning=F, message = F, fig.width=10, fig.height=8}
par(mfrow = c(1,2))
barplot(cc1$cor, main = "Canonical correlations for 'cancor()'", col = "gray")
barplot(cc2$cor, main = "Canonical correlations for 'cc()'", col = "gray")
```


However, the actual estimated linear combinations of the original covariates in both datasets $\mathcal{X}$ and $\mathcal{Y}$ are different. Compare the following: 

```{r}
cc1$xcoef  ### function from standard R instalation
cc2$xcoef  ### function for the R package 'CCA'

cc1$ycoef
cc2$ycoef
```

Using the 'CCA' package one can also take an advantage of nice graphical tools which are available for
the canonical correlation analysis. The idea is to display maximized correlations between transformed variables of the dataset $\mathcal{X}$ and the dataset $\mathcal{Y}$. 

```{r, warning=F, message = F, fig.width=10, fig.height=8}
plt.cc(cc2, var.label = TRUE, ind.names = data[,1])
```


#### <font color = "#880000">Questions</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>How would you interpret the left hand side of the figure?</li>
 
 <li>What is the relationship between the subfigure on the left hand side and the subfigure on the 
 right hand side? </li>
 
 <li>What would be your conclussion about the two available dataset and their mutual relationship? </li>
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

A usefull extension to the classical canonical correlation approach is the regularized version - so called regularized canonical correlation analysis. It is meant for scenarios where the number of parameters $p \in \mathbb{N}$ is greater than the sample size $n \in \mathbb{n}$ (usually $p \gg n$). 
There is R function `rcc()` available in the R package 'CCA'. <br><br>

Another option how to perform canonical correlation in R is to install some additional R packages -- for instance, package 'vegan' (use `install.packages("vegan")` for instalation) and to apply function `cca()`.  Use the help session in R for more details about the function itself or the corresponding package respectively.<br><br>




#### <font color = "#880000">To Do</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>Use this markdown and consider (maybe) the whole dataset on biologial measurements. Compare how 
 much the results changes. </li>
 
 <li>Consider a lower dimensional profile and try to refit the whole canonical correlation analysis on
 a lower dimensional datasets (e.g. two covariates for bio data and two covariates for chemical data). </li>
 
 <li>Try to manually reconstruct the results from the above by using the correponding canonical correlation theory form the lecture. </li>
 
 
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

#### <font color = "#880000">Note</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li><b>Canonical Correspondance Analysis</b><br> 
 An analogous method that the canonical correlation analysis, however, mainly meant for discrete data.<br><br> </li>
 
 <li><b>Redundancy Analysis</b><br></li>
 A similar method than the canonical correlation analysis, however, not symetric. The idea is tu use linear combinations of one dataset 
 to maximize the corelation with the original covariates in the second dataset. Thus, it depends on which data set is used as the orginal one. 
 
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>



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




## <font color = "#880000"><b>Homework Assignment</b></font>
### (Final Submission Deadline: Tuesday 28.05.2018)

Use the command below to instal the `SMSdata` packages with various multivariate datasets: 

```{r, eval = F}
install.packages(pkgs="http://www.karlin.mff.cuni.cz/~hlavka/sms2/SMSdata_1.0.tar.gz", repos=NULL, type="source")
```

Chose one dataset of your choise from the list of all available datasets in the package:

```{r, eval = F}
data(package = "SMSdata") 
```
 
 There are 21 different datasets and you can load each of them by typing and running 
`data(dataset_name)`, for instance:

```{r, eval = F}
data(plasma) 
```

<ul>
<li>Choose one dataset and use some judgement to split the data into two disjoint parts. Alternatively, 
you could also consider two different datasets however, they should be recorded on the same units. 
Apply the canonical correlation analysi. </li>

<li>Finilize your report and make it ready for the final submission. The submission deadline 
is <font color = "red">Tuesday, 28/05/2017, 23:59</font></li>. The report should be submitted by emial 
to (at least) one of the following addresses: 
<ul>
<li>maciak AT karlin.mff.cuni.cz</li>
<li>maciak AT ualberta.ca</li>
</ul>

<li>Each student who wants to receive the course credit for this term must submit 
his/her PDF report (and individual evaluation of all questios and tasks assigned during the term) not later than 
the deadline provided above. </li>

<li>If you need to obtain the course credit before, it is necessary to submit your PDF report at leats two days 
before the day you need to have the credit assigned.</li>
</ul>





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