Редукција димензионалности
Претпоставимо да желимо да визуализујемо \(n\) тачака али имамо \(p\) колона: \(X_1\),…,\(X_p\). Једини начин који смо до сада имали за ту визуализацију јесте био да правимо графике свих могућих парова промјенљивих. Међутим, ту се појављују два велика проблема:
- Имамо \(\binom{p}{2}\) графика, за велико \(p\) то је јако велики број графика који је тежак за праћење.
- Губи се јако велика количина информација на овај начин.
Због тога нам треба метода која то све представља информативније. Једна од таквих метода јесте и метода PCA.
Principal Component Analysis
PCA (енг. Principal Component Analysis), или, на српском, Анализа главних компоненти, јесте једна од најпознатијих метода за рад са проблемом димензионалности. Како је она есенцијална за анализу података било ког типа, обрадићемо је детаљно.
Нека имамо \(n × p\) матрицу \(X\) у којој се налазе наши предиктори. Претпоставимо да су подаци центрирани, то јест да је узорачка средња вриједност сваке колоне једнака нули.
Напомена. PCA претпоставља да нема категоричких предиктора. Ако су они присутни треба их привремено склонити, те радити PCA на осталим предикторима, па онда категоричке додати назад.
Идеја је да извршимо трансформације предиктора \(X_1\),…, \(X_𝑝\) облика
\[Z_k = \sum_{i=1}^p a_{ki}X_i = \mathbf{a}_k^T\mathbf{X}\]
такве да је \(DZ_k = \mathbf{a}_k^T \mathrm{Cov}(\mathbf{X})\mathbf{a}_k\) што веће. Да бисмо имали оптимизациони проблем тражења условног екстремума, додаје се и услов \(\sum_{i}a_{ki}^2=1\),а ово заправо значи да је трансформација задата матрицом \(A = [a_{𝑘𝑙}]\) једна ортогонална трансформација.
Рјешавањем горњег оптимизационог проблема за \(𝑘 = 1\) налазимо коефицијенте \(a_{11}\),…,\(a_{1𝑝}\), па самим тим и вектор \(Z_1\) који зовемо прва главна компонента. На наредној слици је илустрација.
Овиме смо пронашли правац највеће распршености наших података. Сљедећи вектор који тражимо, \(Z_2\) треба да одговара сљедећем правцу највеће распршености, послије првог. Овај вектор зовемо друга главна компонента. Њега тражимо тако што стављамо услов \(||\mathbf{a}_2||_2=1\), тражимо да је он ортогоналан на претходну главну компоненту, дакле \(\mathbf{a}_2^T\mathbf{X} \perp Z_1\) и наравно тражимо да је \(\mathbf{a}_2^T\mathrm{Cov}(\mathbf{X})\mathbf{a}_2\) што веће. Поступак понављамо итеративно.
Може се показати (ЛСМ) да је прва главна компонента заправо сопствени вектор који одговара највећој сопственој вриједности матрице \(\mathrm{Cov}(\mathbf{X})\).
Напомена. Подсјетимо се, коваријациона матрица је ненегативно дефинитна, те има ненегативне сопствене вриједности.
Након тога, уместо “одокативног” одбацивања предиктора можемо просто задржати првих неколико главних компоненти. Пракса је да се задржи онај број који објашњава барем \(80\) посто варијабилитета у подацима. Детаљно извођење свега у вези са PCA биће дато на ЛСМ, укључујући и коваријациону матрицу главних компоненти, те и теоријски губитак који имамо одбацивањем последњих \(𝑙\) компоненти.
Основна мана овог алгоритма јесте та што подаци губе на интерпретабилности. Сваки од оригиналних предиктора имао је неко значење: висину, масу, годиште итд, док ови нови, који су њихове линеарне комбинације, не значе ништа. То је цијена коју смо морали да платимо.
PCA - примјер
Декомпозиција матрице у сврху PCA у \(R\)-у може да буде спектрална, и тако је имплементирана у функцији \(princomp\), док је у функцијама \(prcomp\) и \(PCA\) из \(FactoMineR\) имплементирана сингуларна декомпозиција коваријационе матрице. У суштини је свеједно: ова друга даје мало боље резултате у пракси, те ћемо се ми одлучити за \(prcomp\).
Напомена. Ова функција има аргумент scale али је он подразумевано FALSE! Податке треба или ручно скалирати пре позива ове функције, или аргумент поставити на TRUE.
За базу података одлучићемо се за већ познату \(USArrests\).
data("USArrests")
head(USArrests)## Murder Assault UrbanPop Rape
## Alabama 13.2 236 58 21.2
## Alaska 10.0 263 48 44.5
## Arizona 8.1 294 80 31.0
## Arkansas 8.8 190 50 19.5
## California 9.0 276 91 40.6
## Colorado 7.9 204 78 38.7Свака врста је једна америчка савезна држава, а колоне говоре саме за себе.
Сљедећа ствар коју радимо јесте да скалирамо наше податке.
scaled_df <- apply(USArrests, 2, scale)Сада ћемо извршити PCA над овом базом.
pca_result <- prcomp(scaled_df)Да видимо шта нам је функција вратила.
class(pca_result)## [1] "prcomp"Ова класа је заправо листа повратних вриједности, а њихова имена можемо добити као:
names(pca_result)## [1] "sdev" "rotation" "center" "scale" "x"Промјенљиве center и scale одговарају средњим вриједностима и стандардним девијацијама колона које су постојале прије скалирања (код нас 0 и 1, јер смо скалирали ручно). У то се можемо увјерити позивом
pca_result$center## Murder Assault UrbanPop Rape
## -7.667478e-17 1.111611e-16 -4.332645e-16 8.938163e-17Видимо да су ове вриједности приближно нула. У rotation је уписана матрица \(𝐴\). У sdev уписане су стандардне девијације сваке од главних компоненти, а у x саме главне компоненте.
Треба напоменути да \(prcomp\) подразумијевано прави сопствене векторе који су негативно усмерени (“према доље”), тако да, по потреби, можемо узети њихове супротне вриједности.
Дисперзија сваке од главних компоненти је, наравно, њена стандардна девијација на квадрат, а пропорцију дисперзије објашњене њоме добијамо једноставно:
PVE <- (pca_result$sdev^2) / sum(pca_result$sdev^2)
PVE## [1] 0.62006039 0.24744129 0.08914080 0.04335752sum(PVE)## [1] 1Сада треба одлучити колико главних компоненти задржати. За то је згодан Scree plot приказ који даје функција \(fviz\_eig\) из пакета \(factoextra\).
library(factoextra)
fviz_eig(pca_result, addlabels = TRUE)
У овом случају било би оптимално задржати двије главне компоненте. За кумулативан приказ ад-хок рјешење је
plot(cumsum(PVE), type = "b")
Ово је само демонстрација, у било какве практичне сврхе користили бисмо углађеније графике из \(ggplot2\) пакета.
Да сада мало обратимо пажњу на корелациону структуру нових и старих предиктора. Погледајмо корелациону матрицу оригиналних предиктора.
library(corrplot)
corrplot(cor(scaled_df))
Видимо да је присутна јака линеарна зависност међу предикторима. Ако сада погледамо корелациону матрицу главних компоненти:
corrplot(cor(pca_result$x))
Видимо да су главне компоненте некорелисане, што је наравно добро.
PCA - регресија
У претходном одјељку смо показали шта је PCA и како функционише, а сада ћемо демонстрирати како се овај метод може користити заједно са регресијом, а у сврху предвиђања. Користићемо надалеко чувени Boston скуп података који говори о цијенама некретнина, а желимо да предвидимо medv, средњу вриједност цијене некретнине у хиљадама долара.
Прво ћемо учитати податке.
library(MASS)
data("Boston")Сада ћемо да извршимо PCA на свим колонама осим medv, а њу узимамо за зависну.
boston_preds <- subset(Boston, select = -c(medv))
boston_pca <- prcomp(boston_preds, center = TRUE, scale = TRUE)Можемо кратко бацити поглед на дисперзије сваке од компоненти.
plot(boston_pca)
За прављење регресије користићемо се пакетом caret.
library(caret)Генеричка функција за прављење модела из овог пакета јесте функција \(train\), која има јако згодан аргумент preProcess којем може да се каже како да претпроцесира податке прије саме регресије.
set.seed(123)
model_pcr <- train(medv ~.,
data = Boston,
method = "lm", # linearna regresija
preProcess = c("center", "scale", "pca")
)Да видимо шта је урадио.
model_pcr## Linear Regression
##
## 506 samples
## 13 predictor
##
## Pre-processing: centered (13), scaled (13), principal component
## signal extraction (13)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 506, 506, 506, 506, 506, 506, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 5.401168 0.6604614 3.595097
##
## Tuning parameter 'intercept' was held constant at a value of TRUEsummary(model_pcr)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.018 -2.810 -0.888 1.604 31.392
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.53281 0.22382 100.674 < 2e-16 ***
## PC1 -2.27302 0.09051 -25.113 < 2e-16 ***
## PC2 2.19491 0.18714 11.729 < 2e-16 ***
## PC3 3.50099 0.20098 17.419 < 2e-16 ***
## PC4 1.08068 0.24193 4.467 9.84e-06 ***
## PC5 -2.23308 0.24521 -9.107 < 2e-16 ***
## PC6 -0.67063 0.27632 -2.427 0.01558 *
## PC7 -0.09431 0.30620 -0.308 0.75820
## PC8 -1.04018 0.35598 -2.922 0.00364 **
## PC9 0.14945 0.42573 0.351 0.72570
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.035 on 496 degrees of freedom
## Multiple R-squared: 0.7057, Adjusted R-squared: 0.7003
## F-statistic: 132.1 on 9 and 496 DF, p-value: < 2.2e-16Видимо да је користио све главне компоненте за модел, те видимо да има и оних које нису значајне. Сада ћемо подијелити податке на тренинг и на тест скуп да направимо предвиђање, те да провјеримо оверфит за сваки случај.
boston_train_index <- createDataPartition(1:nrow(Boston), p = 0.8, list = FALSE)
boston_train <- Boston[boston_train_index, ]
boston_test <- Boston[-boston_train_index, ]Сада ћемо да направимо модел само на тренинг скупу.
set.seed(123)
model_pcr_train <- train(medv ~.,
data = boston_train,
method = "lm",
preProcess = c("center", "scale", "pca")
)
model_pcr_train## Linear Regression
##
## 406 samples
## 13 predictor
##
## Pre-processing: centered (13), scaled (13), principal component
## signal extraction (13)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 406, 406, 406, 406, 406, 406, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 5.361084 0.6694754 3.59626
##
## Tuning parameter 'intercept' was held constant at a value of TRUEСада вршимо предвиђање.
boston_test_prediktori <- subset(boston_test, select = -c(medv))
boston_test_target <- subset(boston_test, select = c(medv))[, 1]
# kastovanje u "numeric"
predvidjanje <- predict(model_pcr_train, newdata = boston_test_prediktori)Да видимо грешку.
sqrt(mean((boston_test_target - predvidjanje)^2))## [1] 4.241712Грешка на тренинг скупу била је 5.401168, тако да свакако нисмо у оверфиту. Како видимо да нам се грешке доста машавају, сетимо се да имамо и бољи метод процјене грешке на тренинг скупу, а то је кросвалидација.
ctrl <- trainControl(
method = "cv",
number = 35 # radimo 35-fold
)Направимо нови модел. То је заправо исти модел, само другачије оцјењује грешку.
set.seed(123)
model_pcr_train_2 <- train(medv ~.,
data = boston_train,
method = "lm",
preProcess = c("center", "scale", "pca"),
trControl = ctrl
)
model_pcr_train_2## Linear Regression
##
## 406 samples
## 13 predictor
##
## Pre-processing: centered (13), scaled (13), principal component
## signal extraction (13)
## Resampling: Cross-Validated (35 fold)
## Summary of sample sizes: 394, 394, 394, 394, 394, 395, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 4.934065 0.7352477 3.56909
##
## Tuning parameter 'intercept' was held constant at a value of TRUEsummary(model_pcr_train_2)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.686 -2.805 -1.006 1.519 31.387
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.6520 0.2597 87.227 < 2e-16 ***
## PC1 -2.2026 0.1055 -20.887 < 2e-16 ***
## PC2 2.0314 0.2165 9.381 < 2e-16 ***
## PC3 3.8558 0.2355 16.371 < 2e-16 ***
## PC4 1.4997 0.2763 5.427 1.00e-07 ***
## PC5 -1.2728 0.2864 -4.443 1.15e-05 ***
## PC6 -0.9138 0.3183 -2.871 0.00431 **
## PC7 -0.4105 0.3516 -1.167 0.24372
## PC8 -1.0376 0.4060 -2.556 0.01097 *
## PC9 -0.1813 0.4942 -0.367 0.71389
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.233 on 396 degrees of freedom
## Multiple R-squared: 0.6842, Adjusted R-squared: 0.677
## F-statistic: 95.31 on 9 and 396 DF, p-value: < 2.2e-16Видимо да много боље погађа грешку.
Напомена. За штеловање броја фолдова ја сам користио метод “од ока”. Наравно, у реалном животу, треба користити валидациони скуп, као и друге методе које смо раније помињали (“лакат”).
Сада ћемо још да задржимо само оне главне компоненте које објашњавају \(80\%\) варијабилности унутар података. Подразумијевано се задржава \(95\%\).
set.seed(123)
boston_train_for_pp <- subset(boston_train, select=-c(medv, chas))
boston_train_wait <- subset(boston_train, select=c(medv, chas))
# chas je kategoricka, ne valja za pca
pp <- preProcess(boston_train_for_pp, method = "pca", tresh = 0.8)
# Pravimo pretprocesirani boston_train
boston_train_pp <- predict(pp, boston_train_for_pp)
# Dodajemo mu kolone koje smo izbacili
boston_train_pp <- cbind(boston_train_pp, boston_train_wait)
model_pcr_train_3 <- train(medv ~.,
data = boston_train_pp,
method = "lm",
trControl = ctrl
)
model_pcr_train_3## Linear Regression
##
## 406 samples
## 10 predictor
##
## No pre-processing
## Resampling: Cross-Validated (35 fold)
## Summary of sample sizes: 394, 394, 394, 394, 394, 395, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 4.979478 0.7315984 3.618851
##
## Tuning parameter 'intercept' was held constant at a value of TRUEsummary(model_pcr_train_3)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.936 -2.767 -1.022 1.403 29.060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.4651 0.2681 83.798 < 2e-16 ***
## PC1 -2.2147 0.1050 -21.083 < 2e-16 ***
## PC2 0.7731 0.2255 3.428 0.000671 ***
## PC3 4.4397 0.2410 18.425 < 2e-16 ***
## PC4 1.4819 0.2847 5.205 3.13e-07 ***
## PC5 0.8421 0.3148 2.675 0.007784 **
## PC6 0.3434 0.3494 0.983 0.326206
## PC7 1.0372 0.4042 2.566 0.010659 *
## PC8 0.1843 0.4920 0.375 0.708166
## PC9 -1.1232 0.5324 -2.110 0.035506 *
## chas 2.8093 1.0647 2.639 0.008652 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.21 on 395 degrees of freedom
## Multiple R-squared: 0.6877, Adjusted R-squared: 0.6798
## F-statistic: 86.97 on 10 and 395 DF, p-value: < 2.2e-16Видимо да нисмо много покварили модел (можда је чак и мало бољи) у смислу показатеља квалитета, али смо га доста поједноставили. Ово свакако није никакав примјер квалитетног предвиђања, већ више демонстрација како се PCA користи у комбинацији са регресијом.
Ако се некоме чини да PCA некако одузима везу података од зависне промјенљиве, наредни график би требало да га разувјери.
plot(medv ~ PC1, data = boston_train_pp)
Види се јасна зависност. Такође, види се и зашто је модел лош: зависност није баш линеарна. Ако нацртамо дијагностичке графике, видимо да није добро.
par(mfrow=c(2,2))
plot(model_pcr_train_3$finalModel)
Значи да би модел требало даље надограђивати трансформисањем главних компоненти. Ми нећемо, није нам то циљ. Коме је занимљиво нек проба, добро је за вјежбу.
t-SNE
Замислимо да смо добили базу података коју треба да истражимо. Природна је ствар кренути од цртања. У причи о PCA смо видјели да је цртање пар-по-пар графика неинформативно и рачунски захтјевно. Стога се појавила потреба за методом који податке из вишедимензионог простора “пројектовати” на дводимензионални простор који ми можемо у потпуности да перципирамо.
Прва идеја која човјеку пада на памет јесте да се срачунају сва растојања између оригиналних тачака. Нека, рецимо имамо матрицу чије су врсте \(x_𝑘\), и означимо \(𝑑_{𝑖𝑗} = ||x_𝑖 − x_𝑗||\), гдје је норма одабрана по потреби датог случаја.
Сада се чини смисленим да пројекције \(y_𝑖\) тачака \(x_𝑖\) тражимо тако да се чувају растојања, то јест да за све парове тачака важи да је
\[||y_i-y_j||=d_{ij}.\]
То је, међутим, доказано немогуће, и познато је у машинском учењу као “клетва димензионалности” (Curse of Dimensionality). Ми се овдје нећемо теоријски увјеравати у то, али можемо убиједити себе на кратком примјеру. На наредној слици дати су хистограми пар-по-пар растојања из генерисаних 2-димензионих, 50-димензионих и 500-димензионих нормалних случајних величина (ПСУ). Видимо да како расту димензије, нема малих растојања. Стога је немогуће очекивати да се нешто што изгледа као други или трећи хистограм утопи у дводимензион простор у коме очекујемо први хистограм.
Први одговор на овај проблем била је PCA, коју је смислио Карл Пирсон 1901, а независно од њега развио и именовао Харолд Хотелинг 1930их година. Међутим, PCA је у својој сржи линеарна метода и нема велику способност “развлачења података”. Њена основна идеја је да се цртају две главне компоненте, а остале да се занемарују. Јасно је и неутренираном оку колико ту губитака може настати.
t-SNE је скраћено од t-Distributed Stochastic Neighbour Embedding. За почетак ћемо се упознати са тиме шта је Neighbour Embedding, од сада скраћено NE. Прва и основна премиса за NE јесте да се у потпуности одрекнемо идеје о очувању растојања и да пробамо нешто друго. То нешто друго биће покушај да се очувају односи растојања. То би значило да се тачке са лијевог краја зеленог хистограма пресликавају у лијеви крај плавог хистограма, и слично за десни крај и средину. У преводу, биће нам циљ да тачке које су оригинално близу буду близу и при “пројекцији”, а да оне које су биле оригинално далеко тако и остану. Дакле, није идеја да се очува растојање, већ сусједство, или на енглеском to preserve Nearest Neighbours. Отуда и име.
G. Hinton и S.T. Roweis су 2002. у свом раду поставили темеље алгоритма Стохастичко NE, у даљем тексту SNE. Овај рад до сада је цитиран преко 1464 пута.
Постоји и модификација овог рада (Hinton је аутор оба), који је до сада цитиран више од 23200 пута.
За почетак ћемо дефинисати сличности између опсервација као вјероватноће \(𝑝_{𝑖𝑗}\). Што су две опсервације ближе, њихова сличност је већа. На исти начин ћемо дати и сличности пројекција, \(𝑞_{𝑖𝑗}\). Још ћемо видјети како су ове сличности конкретно дефинисане.
Ако претпоставимо да смо њих задали, сада је идеја да нађемо пројекције које минимизују “разлику” \(𝑝_{𝑖𝑗}\)-ова и \(𝑞_{𝑖𝑗}\)-ова. Сада се поставља ново питање: шта узети за функцију губитака, то јест шта ће да нам мјери разлику сличности пројекција и сличности оригинала. Код SNE-а избор је пао на Кулбак-Лајблерово разилажење, које је дефинисано са
\[L = \sum_{i,j} p_{ij}\log \frac{p_{ij}}{q_{ij}}.\]
Ко упише курс Теорија информације упознаће се боље са овим појмом, а овдје је довољно знати само шта мјери. Битно је уочити и да Кулбак-Лајблерово разилажење није симетрично, па није метрика на простору вјероватносних мјера које су дискретног типа. Зато га и зовемо разилажење. Из формуле се види да је акценат на томе да се плаћа велика цијена ако се блиске тачке пројектују далеко. Цијена обрнуте грешке је нешто мања.
Сада је ред да покажемо како се дефинишу сличности међу тачкама. Узимамо
\[p_{j|i} = \frac{\exp (-||x_i-x_j||^2/2\sigma_i^2)}{\sum_{k\neq i}\exp{(-||x_i-x_k||^2/2\sigma_i^2)}}.\]
Сигме се бирају да се постигне жељена Шенонова ентропија. То ћемо радити на Теорији информације, а овдје се може замислити да умјесто сигме стоји било која константа. Интуитивно: бирамо сигме тако да свака опсервација има \(𝐾\) блиских сусједа, а остале далеке, гдје је \(𝐾\) фиксно, подразумевано 30. Како ове вјероватноће нису симетричне, ми их симетризујемо и добијемо
\[p_{ij} = \frac{p_{i|j}+p_{j|i}}{2n}.\] Ове вјероватноће се по конструкцији сабирају на 1.
У многим ситуацијама може се проћи јефтиније, те се може дефинисати да је \(𝑝_{𝑗|𝑖} = 1/𝐾\) за \(𝐾\) најближих, а нула за остале. Често резултати нису много лошији.
Да видимо још како дефинишемo \(𝑞_{𝑖𝑗}\)-ове. Дефинишемо
\[q_{ij} = \frac{k(||y_i-y_j||)}{\sum_{k\neq l}k(||y_k-y_l||)},\]
гдје SNE бира
\[k(d) = e^{-d^2},\]
а t-SNE бира
\[k(d) = \frac{1}{1+d^2}.\]
Напомена. Основна разлика је у тежини репова.
MNIST dataset
Modified National Institute of Standards and Technology скуп података представља скуп од 70 хиљада сличица 28 × 28 пиксела, лабелованих, које представљају ручно исписане цифре.
Да видимо шта уради PCA, а шта t-SNE.
Примјена оквира компетенција у учењу и развоју запослених
http://alas.matf.bg.ac.rs/~mv18004/ss2_faktorska_analiza.html
Задаци за самосталан рад
Задатак 1. У бази \(mtcars\) налазе се подаци о потрошњи горива и 10 аспеката дизајна и перформанси за 32 аутомобила. Спровести анализу главних компоненти у циљу смањења димензионалности података, те извршити предвиђање на основу постављеног модела.
Задатак 2. У бази \(decathlon2\) (пакет \(factoextra\)) налазе се подаци о перформансама атлетичара на два спортска догађаја (Desctar и OlympicG). База садржи информације о 27 атлетичара са 13 предиктора. Спровести анализу главних компоненти у циљу смањења димензионалности података, те извршити предвиђање на основу постављеног модела.