<table width = 900 border = 0><tr><td>
# <b>NMST539 | Cvičenie 11</b>

## <b>Diskriminačná analýza</b>
### <font color = "#880000"><b>LS 2020/2021 |  11/05/21 | (online výuka)</b></font>
###### <a href = "cvicenie2021_11.Rmd">zdrojový Rmd súbor (kódovanie UTF8)</a>

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


<b>Outline desiateho cvičenia:</b> 

<ul>
<li>základy diskriminačnej analýzy v programe R;</li>
<li>alternatívne prístupy ku klasifikácii;</li>
</ul><br>


<b>Štatisticky program R</b> je k dispozícii (GNU public licence) na adrese <a href = "https://www.r-project.org">https://www.r-project.org</a><br><br>

<b>RStudio</b> -- "user-friendly" interface (jeden z mnohých, ktoré na internete nájdete): <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Užitočné návody a manuály pre prácu s programom R:</b><br>

<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i>  (<a href = "NMFM301/Rmanual2.pdf">PDF súbor</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i>  (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF súbor</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (<a href = "NMFM301/uvodr.pdf" >PDF súbor</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) </li>		
</ul><br>


<b>Odborná literatúra:</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>

<br><br>


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


### <font color = "#880000"><b>1. Diskriminačná analýza v R</b></font>

Diskriminačná analýza je štatistická metóda (nie nutne pouze jedna konkrétna a nie nutne pouze určená  pre mnohorozmerné data), 
ktorej cieľom je klasifikovať objekty do konkrétnych (predom daných) skupín -- t.j. metóda, ktorej cieľom je "naučiť" sa princíp diskriminácie (resp. rozlíšovania)
jednotlivých objektov na základe informácie o príslušnosti do skupín a dodatočných premenných (v prípade, že je dodatočných premenných veľké množstvo, 
mohli by sme hovoriť o mnohorozmernej štatistickej metóde). Výstupom diskriminačnej analýzy je tzv. "diskriminačné pravidlo", ktoré s využitím 
dodatočných premenných umožňuje klasifikáciu nových pozorovaní a to bez znalosti primárnej informácie o príslušnosti nových pozorovaní do už existujúcich skupín. 

 <ul>
 <li><b>Klasifikácia</b> pracuje s datami, kde informácia o príslušnosti do skupín je apriórne známa. Cieľom diskriminačnej analýzy 
 je vytvoriť klasifikačné pravidlo, ktoré lze následne aplikovať na budúce pozorovania bez znalosti príslušnosti týchto pozorovaní do daných skupín. <br> </li>
 
 <li><b>Zhlukovanie</b> pracuje s datami, u ktorých informácia o príslušnosti do skupín nie je apriórne známa a jednotlivé skupiny pozorovaní 
 sú výhradne vytvorené na základe konkrétnej miery podobnosti/nepodobnosti medzi pozorovaniami (resp. zhlukmi).</li>
 </ul>
 
Existuju rôzne postupy a algoritmy pre vytvorenie vhodného diskriminačného pravidla. Jednotlivé metódy môžu byť založené napr. na vierohodnosti (napr. klasický zobecněný lineárny model), apriórnom a aposteriórnom rozdelení (napr. Bayesovo diskriminačné pravidlo), súčte štvorcov medzi skupinami a vrámci skupín (napr. Fisherove LDA diskriminačné pravidlo), alebo rózne ďalšíe metódy zahrňujúce komplexné algoritmy, alebo rôzne expertné zásahy (napr. neuronové siete). 
 
V prípade normálneho modelu s homoskedastickými chybami je jednoduché diskriminačné pravidlo založené na lineárnej forme (LDA -- linear discriminant analysis), t.j. lineárnej kombinácii 
dostupných premenných. V prípade heteroskedastických chýb je diskriminačné pravidlo často postavené na kvadratickej forme. Lineárne diskriminačné pravidlo je v programe R implementováne 
vo funkcii `lda()`, ktorá je dostupná v knižnici 'MASS'. Kvadratické diskriminačné pravidlo je implementované vo funkcii `qda()` (knižnica 'MASS').

Lineárna diskriminačná metóda je úzko previazaná z metódou analýzy rozptylu (ANOVA) a klasickým lineárnym regresným modelom. Metóda ANOVA (resp. lineárny regresný model) pracuje so 
spojitou závislou premennou a diskretnými (resp. spojitými) nezávislými premennými, zatiaľ do diskriminačná analýza pracuje s diskrétnou závislou premennou a spojitími (prípadne kategorickými) nezávislými premennými. V tomto smere je možné lineárne diskriminačné pravidlo porovnať aj s logistickým regresným modelom, ktorý taktiež vysvetľuje závislú kategorickú premennú. Na druhej strane, logistický regresný model sa využíva pre data, ktoré nie sú normálne rozdelené -- čo je vpodstate fundamentálny predpoklad pre LDA. 

V prípade, že predpoklad normality je príliš obmedzujúci, je možné aplikovať tzv. Fisherov postup, ktorý je založený na maximalizácii pomeru súčtu štvorcov medzi skupinami 
a súčtom štvorcov vrámci skupín. Tento postup nevyžaduje ani homoskedastické chyby. Ak budeme predpokladať, že dve skupiny pozorovaní maju príslušné stredné hodnoty $\boldsymbol{\mu}_1$ a $\boldsymbol{\mu_2}$ a variančné-kovariančné matice sú  $\Sigma_1$ a $\Sigma_2$, tak Fisherovo lineárne diskriminačné pravidlo bude definované ako pomer 

$$
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} }.
$$
čo je vpodstate tzv. "signal-to-noise" pomer, ktorý je maximalizovaný pre 

$$
\boldsymbol{w} = (\Sigma_1 + \Sigma_2)^{-1} (\boldsymbol{\mu}_2 - \boldsymbol{\mu}_1).
$$

<br><br>

Pre ilustráciu využijeme datový súbor  `mtcars` a apriórna informácia o príslušnosti jednotlivých pozorovaní do skupín bude obsiahnutá vrámci premennej `cyl` -- čo je počet válcov v danom aute (k dispozícii sú pozorovania so 4, 6 a 8 válcami). Pre vytvorenie lineárneho diskriminačného pravidla využijeme niektoré spojité premenné. <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)
```

Z teoretického hľadiska predstavuje lineárne diskriminačné pravidlo rozdelenie $n$ rozmerného výberového prietoru pomocou nadroviny určenej lineárnou kombináciou daných premenných. V tomto konkrétnom prípade uvažujeme 5 rozmerný priestor $\mathbb{R}^5$ a náhodné vektory definované ako

$$
\boldsymbol{X} = (X_1, \dots, X_5)^\top = (mpg, hp, qsec, drat, disp)^\top.
$$

Lineárne diskriminačné pravidlo rozdeľuje lineárny priestor $\mathbb{R}^5$ pomocou dvoch nadrovin (pretože kategorizujeme pozorovania do troch skupín) 
a toto delenie si môžeme vizualizovať pomocou príkazu  `partimat()` z knižnice `klaR`. Výsledne rozdelenie priestoru $\mathbb{R}^5$ je znázornené
pomocou desiatich dvojrozmerných scatterplotov, kde každý zobrazuje konkrétnu dvojicu pôvodných premenných. 


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

Prípadne analogický graf bez explicitného vykreslenia nadrovín rozdeľujúcich výberový priestor: 

```{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"><b>Samostatne</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>Pre lepšiu predstavu ako funguje lineárne diskriminačné pravidlo, môžeme využiť pouze jednú premennú. 
 
 ```{r}
lda2 <- lda(cyl ~ mpg, data = mtcars)
lda2cv <- lda(cyl ~ mpg, data = mtcars, CV = T)

 ```

Diskriminačné pravidlo je definované pomocou tzv. škálovacieho faktoru -- scaling factor. Príslušná transformácia teda je 

```{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>V grafoch výššie je vidieť jednak rozdelenie pozorovaní do skupín a tiež samotné diskriminačné pravidlo. Jednotlivé pozorovania 
 patria do jednej z troch skupin vizualizovanej pomocou rôznych farieb. Budúce pozorovanie, ktoré neobsahuje informáciu o počte válcov, ale 
 obsahuje informáciu o spotrebe (t.j. premenná `mpg`) je možné zaradiť do jednej z troch existujúcich skupín na základe odhadnutých pravdepodobnosti príslušnosti do danej skupiny. 
 </li>
 
 <li>Uvažujte analogický príklad s využitím dvoch rôznych premenných pre vytvorenie lineárneho klasifikačného pravidla.</li>
 
 <li>Využijte nové pozorovania, ktoré neboli k dispozícii pri vytvorení klasifikačného pravidla (napr. využijte pouze určitú proporciu 
 dat pre vytvorenie klasifikačného pravidla, alebo si nové pozorovania vytvorte manuálne) a klasifikujte ich. </li>
 
 <li>Ako by celý postup vyzeral v prípade využitia kvadratického diskriminačného pravidla (príkaz `qda()`)?
 

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

### <font color = "#880000"><b>2. Alternatívne metódy pre klasifikáciu</b></font>

#### <font color = "#880000"><b>Mowers data: </b></font>

Pre ilustráciu alternatívnych prístupov ku klasifikácii využijeme umelo vytvorený (balancovaný) datový súbor, ktorý obsahuje pouze tri rôzne premenné: 
pre každe pozorovanie (majiteľ pozemku) je k dispozícii informácia o veľkosti, resp. rozlohe pozemku (spojitá premenná), informácia o celkovom ročnom príjime majiteľa (opäť spojitá premenná) a informácia o príslušnosti pozorovania do jednej z dvoch skupín -- kategorická premenná rozlíšujúca majiťeľov, ktorý majú manuálnu kosačku a majiťeľov, ktorý maju 
samohybnú kosačku (tzv. mower tractor). <br><br>

Idea je nasledujúca: bohatý majitelia, prípadne majitelia s veľkou rozlohou pozemku by teoreticky mohli mať záujem zakúpiť si motorizovanú samohybnú kosačku. 

Na základe týchto informácii sa pokusíme nájsť diskriminačné previdlo (založené na rôznych štatistických aj neštatistických postupoch) pre typ kosačky. 
Použijeme štyri rôzne metódy: zobecnený lineárny model; linárnu diskriminačnú analýzu; kvadratickú diskriminačnú analýzu a lineárny regresný model. 



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

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

Všimnite si, že LDA prístup a klasická lineárna regresia dávajú ekvivalentné (t.j. stejné) výsledky. To platí za predpokladu, že máme k dispozícii presne dve skupiny, 
pre účely lineárnej regresie sú skupiny označené hodnotami mínus jedna a jedna a navyše, data sú vyvážené (t.j. stejný počet pozorovaní v každej skupine -- v tomto konkrétnom 
prípade 12 majiteľov s manuálnou kosačkou a 12 majiteľov s motorizovanou samohybnou kosačkou). Klasifikácia pomocou GLM je mierne odlišná. Výrazne iné výsledky dáva kvadratické 
diskriminačné pravidlo, čo je ale v tomto prípade celkom očakávané (zvyšné tri postupy su lineárne nielen v zmysle parametrov lineárnej kombinácie, ale aj v zmysle geometrie). 

<br><br>

Iný postup, nie celkom založený na štatistickej teórii, by bolo využiť neuronovú sieť. Neuronová sieť je sofistikovaný, ale často tzv. 'black-box' algoritmus, ktorý 
je navýše výrazne nelineárny a nekonvexný. Okrem toho výsledok výrazne závisí od počiatočnej konfigurácii algoritmu (vstupné, resp. počiatočné hodnoty). Pre ilustráciu 
porovnajte nasledujúcich šesť obrázkov, ktoré predstavujú klasifikáciu pomocou neuronovej siete. 

<br><br>

Pre spustenie príkazov a výpočet pomocou neuronovej siete, je nutná inštalácia knižnice `install.library(nnet)`. 

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

for (k in 1:6) {
  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)
}
```

Základna nevýhoda neuronej siete je zrejmá z obrázkov vyššsie -- výrazna nestabilita výsledkov. Aplikácia na stejné data dáva zakaždým iné (výrazne iné) výsledky, čo je dané faktom, že neuronové siete výrazne závisia na voľbe počiatočného nastavenia -- počiatočných hodnôt parametrov siete -- keďže sa jedna o nekonvexnú funkciu. 
V praxi sa preto postupuje spôsobom opakovaného aplikovania neuronovej siete a následného priemerovania získaných výsledkov (z ktorých mnohé vyzerajú nerealistické, alebo 
aspoň hodne nepraktické). <br><br>




<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>Domáca (samostatná) úloha </b></font>
### (Deadline: Cvičenie č.12 / 18.05.2021)

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

<ul>
<li>Vyberte si vhodný datový súbor (napr. data, ktoré ste uvažovali v prvej, štvrtej, alebo predchádzajúcich štyroch úlohach); </li>

<li>Použijte LDA prístup pre klasifikáciu pozorovaní do skupín. Skupiny definujte buď pomocou premennej, ktorá je v datach už k dispozícíí, 
alebo pomocou informácie, ktorú o datach viete napr. z predchádzajúcich úloh (t.j. napr. vytvorte novú premennú, 
ktorá bude reflektovať výsledky faktorovej analýzy).</li>

<li>Analogickú klasifikáciu vytvorte využitím niektorého alternativného postupu (napr. GLM, QDA, linárny model, neuronová sieť) a získane výsledky 
porovnajte a aspoň stručne okomentujte. </li>

<li> Riešenie umiestnite na svoju webovú stránku, <font color = "red">najneskôr v utorok, 18.05.2021, do 14:00</font>.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

<br><br>

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