NMST539 | Cvičenie 6

LS 2020/2021 | 06/04/21 | (online výuka)

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.

(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.

---