NMST539 | Cvičenie 11 Diskriminačná analýza LS 2020/2021 | 11/05/21 | (online výuka) zdrojový Rmd súbor (kódovanie UTF8) Outline desiateho cvičenia: základy diskriminačnej analýzy v programe R; alternatívne prístupy ku klasifikácii; Štatisticky program R je k dispozícii (GNU public licence) na adrese https://www.r-project.org RStudio – “user-friendly” interface (jeden z mnohých, ktoré na internete nájdete): RStudio. Užitočné návody a manuály pre prácu s programom R: Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor) Komárek, A.: Základy práce s R. (PDF súbor) Kulich, M.: Velmi stručný úvod do R. (PDF súbor) De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808) Odborná literatúra: Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015 Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979. 1. Diskriminačná analýza v R 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. Klasifikácia 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. Zhlukovanie 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). 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). \] 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é. 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. 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: pairs(mtcars[c("mpg", "hp", "qsec", "drat", "disp")], pch = 22, bg = c("red", "yellow", "green")[unclass(as.factor(mtcars$cyl))]) Samostatne 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 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") 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. Uvažujte analogický príklad s využitím dvoch rôznych premenných pre vytvorenie lineárneho klasifikačného pravidla. 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. Ako by celý postup vyzeral v prípade využitia kvadratického diskriminačného pravidla (príkaz qda())? 2. Alternatívne metódy pre klasifikáciu Mowers data: 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). 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. Mowers <- read.table("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/mowers.txt", sep = "", header = T) head(Mowers) ## income lot riding ## 1 60.0 18.4 1 ## 2 85.5 16.8 1 ## 3 64.8 21.6 1 ## 4 61.5 20.8 1 ## 5 87.0 23.6 1 ## 6 110.1 19.2 1 ### 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: 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). 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. Pre spustenie príkazov a výpočet pomocou neuronovej siete, je nutná inštalácia knižnice install.library(nnet). require(nnet) ## Loading required package: 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) } ## # weights: 201 ## initial value 21.812101 ## iter 10 value 15.060353 ## iter 20 value 12.751390 ## iter 30 value 8.491047 ## iter 40 value 5.519982 ## iter 50 value 2.748038 ## iter 60 value 0.079289 ## iter 70 value 0.001065 ## final value 0.000090 ## converged ## # weights: 201 ## initial value 18.879395 ## iter 10 value 15.771013 ## iter 20 value 13.606618 ## iter 30 value 6.122672 ## iter 40 value 5.958311 ## iter 50 value 5.422826 ## iter 60 value 4.660722 ## iter 70 value 4.620938 ## iter 80 value 4.531911 ## iter 90 value 4.525809 ## final value 4.525783 ## converged ## # weights: 201 ## initial value 22.510072 ## iter 10 value 14.994248 ## iter 20 value 12.824509 ## iter 30 value 9.325663 ## iter 40 value 4.375754 ## iter 50 value 1.217117 ## iter 60 value 0.001801 ## final value 0.000071 ## converged ## # weights: 201 ## initial value 18.560394 ## iter 10 value 14.736960 ## iter 20 value 13.653053 ## iter 30 value 12.633303 ## iter 40 value 12.613236 ## final value 12.613190 ## converged ## # weights: 201 ## initial value 23.384901 ## iter 10 value 16.210271 ## iter 20 value 13.891970 ## iter 30 value 12.600155 ## iter 40 value 9.730337 ## iter 50 value 4.800069 ## iter 60 value 4.121253 ## iter 70 value 2.781828 ## iter 80 value 2.774573 ## iter 90 value 2.689207 ## iter 100 value 1.492312 ## final value 1.492312 ## stopped after 100 iterations ## # weights: 201 ## initial value 17.192270 ## iter 10 value 15.527719 ## iter 20 value 13.765769 ## iter 30 value 13.345368 ## iter 40 value 4.465187 ## iter 50 value 0.631217 ## iter 60 value 0.001096 ## final value 0.000068 ## converged 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é). Domáca (samostatná) úloha (Deadline: Cvičenie č.12 / 18.05.2021) Vyberte si vhodný datový súbor (napr. data, ktoré ste uvažovali v prvej, štvrtej, alebo predchádzajúcich štyroch úlohach); 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). 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. Riešenie umiestnite na svoju webovú stránku, najneskôr v utorok, 18.05.2021, do 14:00.