NMST539 | Cvičenie 6 Metóda hlavných komponent | Teoretická časť LS 2020/2021 | 06/04/21 | (online výuka) zdrojový Rmd súbor (kódovanie UTF8) Outline šiesteho cvičenia: metóda hlavných komponent (PCA- Principal Component Analysis); teoretické aspekty metódy hlavných komponent; súvislosť hlavných komponent a lineárnej regresie; 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. Metóda hlavných komponent / teoretické aspekty Metóda hlavných komponent (PCA) je mnohorozmerná (prevažne exploratívna) štatistická metóda založena na teoretickej variančnej-kovariančnej matici náhodného vektoru \(\boldsymbol{X} \in \mathbb{R}^{p}\). V praxi sa samozrejme metóda hlavných komponent aplikuje na základe empricikého odhadu teoretickej variančnej-kovariančnej matice, ktorý získame na základe náhodného výberu \(\boldsymbol{X}_1, \dots, \boldsymbol{X}_n\) z rovnakého rozdelenia, ako pôvodný náhodný vektor \(\boldsymbol{X} \in \mathbb{R}^{p}\). Jednoducho môžeme hlavné komponenty definovať ako ‘vhodné’ lineárne kombinácie pôvodných zložiek náhodneho vektoru \(\boldsymbol{X} = (X_1, \dots, X_p)^\top\). Hlavná myšlienka spočíva v definovaní/určení \(\ell < p\) hlavných komponent (lineárnych kombinácii zložiek \(X_1, \dots, X_p\)), ktoré využijeme namiesto \(p \in \mathbb{N}\) pôvodných zložiek (tzv. redukcia dimenzionality). Nové zložky – \(\ell\) hlavných kompoment – sú v určitom zmysle najlepšie možné. Hlavné kompomenty sú orgogonálne a navyše sú usporiadané vzhľadom k postupne klesajúcej variabilite jednotlivých kompoment. V praktických úlohach často požadujeme, aby \(\ell \ll p\), čo predstavuje výraznú redukciu pôvodnej dimenzie. Na druhej strane sa ale snažíme zachovať čo najväčšiu mieru z pôvodnej informácie. Z formálneho hľadiska je možné hlavne komponenty definovať rôzne: ortogonálne lineárne kombinácie zložiek náhodného vektoru s postupne klesajúcou variabilitou; hlavné komponenty definujú najlepší lineárny podpriestor pre aproximáciu dat; z geometrického hľadiska sa jedná o EIV rozklad teoretickej (alebo empirickej) variančnej-kovariančnej matice; z pohľadu dat sa jedná o SVD rozklad \(n \times p\) rozmernej datovej matice; Z pohľadu štatistika sú hlavné komponenty definované ortogonálne lineárne kombinácie pôvodných premenných \((X_1, \dots, X_p)^\top\) v náhodnom vektore \(\boldsymbol{X} \in \mathbb{R}^{p}\), tak aby tieto lineárne kombinácie boli vzájomne nekorelované, s nulovou strednou hodnotou a diagonálnou variančnou-kovariančnou maticou, ktorá ma na hlavnej diagonále vlastné čísla variančnej-kovariančnej matice náhodného vektoru \(\boldsymbol{X} \in \mathbb{R}^{p}\). V tomto zmysle sa jedná o maximalizačnú úlohu, kde prvá hlavná komponenta \(\boldsymbol{\delta}_1^\top \boldsymbol{X}\) (t.j. lineárná kombinácia zložiek vektora \(\boldsymbol{X} \in \mathbb{R}^{p}\)) je definovaná pre \(\boldsymbol{\delta}_1 \in \mathbb{R}^p\) také, že \[ \boldsymbol{\delta}_1 = Argmax_{\boldsymbol{\{\delta};~\|\boldsymbol{\delta}\| = 1 \}} ~~Var (\boldsymbol{\delta}^\top \boldsymbol{X}). \] Druhá hlavná komponenta \(\boldsymbol{\delta}_2^\top \boldsymbol{X}\) je definovaná ako analogická maximalizačná uloha, ale pri dodatočnej podmienke, že sa jedná o ortogonálne lineárne kombinácie, t.j. \(\boldsymbol{\delta}_1^\top\boldsymbol{\delta}_2 = 0\), resp. náhodné veličiny \(\boldsymbol{\delta}_1^\top \boldsymbol{X}\) a \(\boldsymbol{\delta}_2^\top \boldsymbol{X}\) sú nekorelované, t.j., že platí \[ Cov(\boldsymbol{\delta}_1^\top \boldsymbol{X}, \boldsymbol{\delta}_2^\top \boldsymbol{X}) = 0. \] Na hlavné komponenty sa ale môžeme pozerať aj z geometrického pohľadu projekcií do lineárných podpriestorov a hlavne konceptu lineárnej regresie (prípadne tzv. totálne najmenších štvorocov – TLS / Total Least Squares). 2. Hlavné komponenty, ortogonálne projekcie a regresia Nech \(\boldsymbol{Y} \in \mathbb{R}^n\) je nejaký náhodný vektor dĺžky \(n \in \mathbb{N}\) a nech \(\mathbb{X}\) predstavuje (náhodnú) maticu napr. vysvetlujúcich premenných, ktorá je typu \(n \times p\). Vo všeobecnosti predpokládame, že \(p \in \mathbb{N}\) (väčšinou počet rôznych premenných) je menšie, ako \(n\) (väčšinou počet pozorovaní). Klasická úloha lineárnej regresie spočíva vo využití informácie obsiahnutej v matici \(\mathbb{X}\) k vhodnej aproximácii \(n\) rozmerného náhodného vektoru \(\boldsymbol{Y}\) (tzv. vektor závislej premennej) v \(p\) rozmernom lineárnom podpriestore \(\mathcal{L}(\mathbb{X})\) – lineárnom podpriestore \(R^n\) generovanom stĺpcami matice \(\mathbb{X}\). Ak označíme ortogonálnu projekciu náhodného vektoru \(\boldsymbol{Y}\) do \(p\) rozmerného lineárneho podpriestoru \(\mathcal{L}(\mathbb{X})\) ako \(P_{\mathbb{X}}(\boldsymbol{Y})\), tak potom klasický lineárny regresný model pre vektor závislej premennej \(\boldsymbol{Y}\) a nezávisle premenné v stĺpcoch matice \(\mathbb{X}\) vo forme \(P_{\mathbb{X}}(\boldsymbol{Y}) = a_1 \boldsymbol{X}_1 + \dots + a_p \boldsymbol{X}_p\) (lineárna kombinácia stĺpcov matice \(\mathbb{X}\)) dáva najmenšiu možnú štvorcovú chybu v zmysle, že výraz \[ \| \boldsymbol{Y} - P_{\mathbb{X}}(\boldsymbol{Y}) \|_2^2 \] je minimalizovaný (pripomeňte si Gauss-Markovovu vetu pre tzv. BLUE (Best Linear Unbiased Estimate) odhad v lineárnej regresii). Príslušná (idempotentná) projekčná matica je definovaná predpisom \[ \mathbb{H} = \mathbb{X}(\mathbb{X}^\top \mathbb{X})^{-1}\mathbb{X}^\top = \mathbb{Q}\mathbb{Q}^\top, \] kde matica \(\mathbb{X}\) obsahuje vektory (v stĺpcoch) generujúce lineárný podpriestor \(\mathcal{L}(\mathbb{X})\) resp. stĺce matice \(\mathbb{Q}\) predstavujú ortogonálnu bázu lineárneho podpriestoru \(\mathcal{L}(\mathbb{X})\) . Samozrejme platí, že \[ P_{\mathbb{X}}(\boldsymbol{Y}) = \mathbb{H}\boldsymbol{Y} = \mathbb{X}(\mathbb{X}^\top \mathbb{X})^{-1}\mathbb{X}^\top \boldsymbol{Y} = \mathbb{Q}\mathbb{Q}^\top \boldsymbol{Y} = P_{\mathbb{Q}}(\boldsymbol{Y}) \] Koncept lineárnej regresie (t.j. ortogonálnej projekcie do lineárneho podpriestoru generovaného stĺpcami matice \(\mathbb{X}\)) môže byť zobecnený na problém hľadania najlepšej ortogonálnej projekcie \(\boldsymbol{Y}\) do lineárneho podpriestoru, ktorý je generovaný stĺpcami nejakej obecnej matice \(\mathbb{B}\) typu \(n \times p\), tak že platí \(\mathbb{B}^\top\mathbb{B} = \mathbb{I}\). Označme obecne takúto projekciu ako \(P_{\mathbb{B}}(\boldsymbol{Y})\). Potom platí, že \(P_{\mathbb{B}}(\boldsymbol{Y}) = \mathbb{B}\mathbb{B}^\top\boldsymbol{Y}|\) a, navyše, lze ukázať, že zároveň platí aj \[ \|\boldsymbol{Y} - P_{\mathbb{\Gamma}}(\boldsymbol{Y})\|_2^2 \leq \|\boldsymbol{Y} - P_\mathbb{B}(\boldsymbol{Y})\|_2^2 \] kde \(\Gamma\) je matica typu \(n\times p\) ktorá obsahuje (v stĺpcoch) prvých \(p\) vlastných vektorov teoretickej variančnej-kovariančnej matice \(\Sigma\) náhodného vektoru \(\boldsymbol{Y}\) (tzv. Eckart a Young (1936) a Mirsky (1960) approximácia). Pripomeňme, že v súvislosti s konceptom lineárnej regresie – t.j. projekciou \(P_{\mathbb{X}}(\boldsymbol{Y})\) – predpokládame vlasnosť linearity (BLUE). V nasledujúcom príklade použijeme tri konkrétne ortogonálne projekcie náhodného vektoru \(\boldsymbol{Y}\) do lineárneho podpriestoru, ktorý je generovaný stĺpcami vhodnej matice typu \(n \times p\). Rovnaké data – t.j. náhodný vektor \(\boldsymbol{Y}\) – použijeme pre všetky projekcie, ale lineárny podpriestor, do ktorého bude vektor projektovaný, bude zakaždým iný. Výsledné ortogonálne projekcie sa budu preto vájomne líšiť. Ktorá ortogonálna projekcia je najlepšia (a v akom zmysle)? set.seed(1234) x <- seq(0,1, length = 100) y <- 0 + 2 * x + rnorm(100) par(mfrow = c(1,3)) ### regression of Y on X plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection of Y onto L(X)", ylim = c(-1.5, 4)) m <- lm(y ~ x) abline(m, col = "red", lwd = 2) index <- c(10, 60, 80) for (i in 1:100){ lines(c(x[i], x[i]), c(y[i], m$fitted[i]), col = "blue", lty =2) } ### regression of X on Y plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection of X onto L(Y)", ylim = c(-1.5, 4)) m <- lm(x ~ y) lines((x - m$coeff[1])/m$coeff[2] ~ x, col = "red", lwd = 2) index <- c(10, 60, 80) for (i in 1:100){ lines(c(x[i], m$coeff[1] + m$coeff[2] * y[i]), c(y[i], y[i]), col = "blue", lty =2) } ### regression onto the first PC ### data standardization means <- apply(cbind(x,y), 2, mean) sdErs <- apply(cbind(x,y), 2, sd) cd <- cbind((x - means[1])/sdErs[1], (y - means[2])/sdErs[2]) ### first eigen vector firstEigenVector <- eigen(cov(cd))$vectors[,1] ### projection onto the linear space generated by the first eigen vector project2pc1 <- t(firstEigenVector) %*% t(cd) ### rotation back to the original coordinates rotation <- as.numeric(project2pc1) %*% t(firstEigenVector) ### rescaling to the original mean/variance rescaled <- cbind(rotation[,1] * sdErs[1] + means[1], rotation[,2] * sdErs[2] + means[2]) plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection onto 1.PC", ylim = c(-1.5,4), xlim = c(-0, 1)) ### first PC lines(rescaled[,2] ~ rescaled[,1], col = "red", lwd = 2) ### orthogonal projections for (i in 1:100){ lines(c(x[i], rescaled[i,1]), c(y[i], rescaled[i,2]), col = "blue", lty =2) } Na základe predchádzajúceho príkladu môžeme obecne povedať, že pri hľadaní hlavných komponent sa dá namiesto spektrálnej dekompozície (EIV) teoretickej variančnej-kovariančnej matice náhodného vektoru \(\boldsymbol{Y}\) postupovať vrámci konceptu klasickej lineárnej regresie a iteratívnym postupom (alternáciou medzi dvoma podobnými lineárnymi regresnými úlohami) sa dopracovať k ekvivalentnému riešeniu – hlavným komponentám. Pre háhodný vektor závislej premennej \(\boldsymbol{Y} \in \mathbb{R}^n\) múžeme tieto dve úlohy formulovať následovne: \[ \|\boldsymbol{Y} - P_\mathbb{B}(\boldsymbol{Y})\|_2^2 = \|\boldsymbol{Y} - \mathbb{B}\mathbb{B}^\top\boldsymbol{Y}\|_2^2 = \|\boldsymbol{Y} - \mathbb{B} \boldsymbol{v}\|_2^2 = \sum_{i = 1}^{n} [Y_{i} - \boldsymbol{b}_i^\top\boldsymbol{v}]^2, \] kde \(\boldsymbol{v} = \mathbb{B}^\top\boldsymbol{Y}\), pre \(i = 1, \dots, n\), a \(\boldsymbol{b}_i\) je príslušný řádek matice \(\mathbb{B}\). Musíme si uvedomiť, že ani matica \(\mathbb{B}\) ani vektor \(\boldsymbol{v}\) nie sú známe (preto potrebujeme riešiť dve vzálomne alternujúce regresné problémy). Výraz \[ \sum_{i = 1}^{n} [Y_{i} - \boldsymbol{b}_i^\top\boldsymbol{v}]^2 \] budeme preto minimalizovať aj vzhľadom k neznámej ortogonálnej matici \(\mathbb{B}\) (v jednom kroku), tzn., že platí \(\mathbb{B}^\top\mathbb{B} = \mathbb{I}\) a tiež vzhľadom k neznámemu vektoru \(\boldsymbol{v}\) (v následujúcom kroku). Oba kroky dostatočne dlho alternujeme, až kým nie je dosiahnutá konvergencia. Hlavná myšlienka samozrejme spočíva v tom, že sa snazíme nájsť najlepšiu ortogonálnu projekciu do lineárneho podpriestoru, ktorý je generovaný orgogonálnou maticou \(\mathbb{B}\), tak, aby súčet štvorcových chýb bol najmenší možný. Z teórie lineárnej regresie vieme, že za platnosti dodatočného predpokladu linearity je najlepším odhadom práve model lineárnej regresie. V obecnom prípade (bez predpokladu linerity) ale lze získať lepšiu ortogonálnu projekciu. Ak bude matica \(\mathbb{B}\) v jednotlivých stĺpcoch obsahovať prvých \(p \in \mathbb{N}\) hlavných kompoment, tak získame najlepšiu možnú ortogonálnu projekciu náhodného vektoru \(\boldsymbol{Y}\) do \(p\) rozmerného lineárneho podpriestoru (samozrejme všetko pouze v zmysle najmenšej štvorcovej chyby). Jednoduchý príklad alternujúcej regresie od Matíasa Salibiána Barreru. alter.pca.k1 <- function(x, max.it = 500, eps=1e-10) { n2 <- function(a) sum(a^2) p <- dim(x)[2] x <- scale(x, scale=FALSE) it <- 0 old.a <- c(1, rep(0, p-1)) err <- 10*eps while( ((it <- it + 1) < max.it) & (abs(err) > eps) ) { b <- as.vector( x %*% old.a ) / n2(old.a) a <- as.vector( t(x) %*% b ) / n2(b) a <- a / sqrt(n2(a)) err <- sqrt(n2(a - old.a)) old.a <- a } conv <- (it < max.it) return(list(a=a, b=b, conv=conv)) } set.seed(1234) n <- 20 p <- 5 x <- matrix(rt(n*p, df=2), n, p) Prvú hlavnú kompomentu pomocou alternujúcich lineárnych regresných modelov získame ako alter.pca.k1(x)$a ## [1] 0.08363314 -0.95027213 -0.01427383 0.11629502 -0.27615231 Jedná sa o stĺpcový vektor matice \(\mathbb{B}\), ktorý nam dá minimálnú štvorcovú chybu projekcie do jednorozmerného lineárneho podpriestoru, ktorý je generovaný prvým vlastnym vektorom variančnej-kovariančnej matice – čo múžeme priamo porovnať s využitím EIV alebo SVD rozkladu: svd(cov(x))$u[,1] ## [1] -0.08363314 0.95027213 0.01427383 -0.11629502 0.27615231 eigen(cov(x))$vectors[,1] ## [1] -0.08363314 0.95027213 0.01427383 -0.11629502 0.27615231 Je nutné si uvedomiť, že rozdiel v znamienku je v tomto prípade irelevantný. Z hľadiska praktických aplikacii môže byť užitočné porovnať fungovanie a časovú náročnosť oboch postupov napr. pomocou simulačnej štúdie. Využijeme výrazne vyšší počet dimenzii, aby bolo porovnanie viac relevantné a výpovedné. n <- 2000 p <- 500 x <- matrix(rt(n*p, df=2), n, p) system.time( tmp <- alter.pca.k1(x) ) ## user system elapsed ## 0.340 0.024 0.364 system.time( e1 <- svd(cov(x))$u[,1] ) ## user system elapsed ## 0.846 0.005 0.852 Výpočetná časová náročnosť sa zdá byť v prospech alternujúcej regresie. Obecne platí, že v programe R je nutné dbať na maticové operácie a pristupovať k ním efektívne, čo do časovej náročnosti výpočtu. 3. Metóda hlavných komponent v programe R Metóda hlavných komponent je v programe R implementovaná pomocou príkazu princomp() ktorý využíva spektrálný rozklad (EIV) teoretickej (resp. empirickej) variančnej-kovariančnej matice. Pre konkrétne data \(\mathcal{X} = (\boldsymbol{X}_1, \dots, \boldsymbol{X}_p)\) s odhadnutou variančnou-kovariančnou maticou \(\mathcal{S}\) získame spektrálny rozklad – t.j. \(EIV(\mathcal{S})\) (analogicky, ak je známa teoretická matica, získame \(EIV(\mathcal{\Sigma})\)). Alternatívne a pri praktických (t.j. empirických) úlohach aj ekvivalentne je možné v programe R využiť aj príkaz prcomp(), ktorý funguje výhradne na empirickom princípe a využíva singulárný rozklad (SVD) datovej matice \(\mathcal{X}\) – teda \(SVD(\mathcal{X})\)). Označme oba rozklady ako \[ EIV(\mathcal{S}) = \Gamma \Lambda \Gamma^\top \quad \quad \textrm{a analogicky} \quad \quad SVD(\mathcal{X}) = UDV^\top. \] Ekvivalentnosť oboch rozkladov je viac-menej okamžitá: \[ n \mathcal{S} = \mathcal{X}^\top\mathcal{X} = VDU^\top UD V^\top = V D^2 V^\top = \Gamma \widetilde{\Lambda} \Gamma^\top = EIV(n\mathcal{S}), \] pre \(\Gamma \equiv V\) a \(\widetilde{\Lambda} \equiv D^2\). Taktiež platí, že \[EIV(n\mathcal{S}) = \Gamma \widetilde{\Lambda} \Gamma^\top = \Gamma (n \Lambda) \Gamma^\top = n \Gamma \Lambda \Gamma^\top = n EIV(\mathcal{S}).\] Využitie metódy hlavných komponent z praktického hľadiska (prevažne ako exploratívny a dimenziu redukujúci nástroj) bude predmetom následujúceho cvičenia. Samostatne Pre fungovanie hlavných komponent je podstatné porozumenie fungovania vlastných čísel a vlastných vektorov. Interaktívne grafické ilustrácie pre pripomenutie je možné nájsť napr. tu: Vlastné čísla a vlastné vektory matice: http://setosa.io/ev/eigenvectors-and-eigenvalues/ Hlavné komponenty v 2D a v 3D: http://setosa.io/ev/principal-component-analysis/ Domáca (samostatná) úloha (Deadline: Siedme cvičenie / 13.04.2021) Uvažujte náhodný výber \((X_1,Y_1)^\top, \dots (X_n,Y_n)^\top\) generovaný, resp. simulovaný z nejakého vhodného dvojrozmerného rozdelenia. Náhodný vektor \(\boldsymbol{Y} = (Y_1, \dots, Y_n)^\top\) lze chápať aj ako element v \(n\) rozmernom lineárnom priestore \(\mathbb{R}^n\). Budu nás zaujímať tri rôzne ortogonálne projekcie – zakaždým do iného dvojrozmerného lineárneho podpriestoru v \(\mathbb{R}^n\). Pre každú projekciu explicitne spočítajte štvorcové chyby (t.j. empiricky overte kvalitu danej projekcie vhľadom k teoretickým očakávaniam). Projekcie sa pokúste nejak vizualizovať. Zvoľte si ľubovolný dvojrozmerný lineárny podpriestor generovaný stĺpcovými vektormi vhodnej matice typu \(n \times 2\) (ale uvažujte obecne iný lineárny podpierstor, ako v nasledujúcich dvoch úlohach). Pomocou projekčnej matice data ortogonálne projektujte do vami zvoleného lineárneho podpriestoru. Využijte metodológiu lineárnej regresie a aplikujte ortogonálnu projekciu na regresnú přímku, t.j. do lineárneho podpriestoru ktorý je generovaný stĺpcami regresnej matice \(\mathbb{X} = (\boldsymbol{1}, \boldsymbol{X})\), kde prvý stĺpec je stĺpec jednotiek a druhý stĺpec je \(\boldsymbol{X} = (X_1, \dots, X_n)^\top\). Nájdite ortogonálnu projekciu do dvojrozmerného lineárneho podpriestoru, ktorý je daný prvou a druhou principálnou kompomentou. Pokúste sa túto úlohu riešiť bez pomoci príkazov princomp() a prcomp(). Riešenie umiestnite na svoju webovú stránku, najneskôr v utorok, 13.04.2021, do 14:00.