Podczas zajêæ bêd± wykorzystywane biblioteki klaR
, MASS
i caret
.
Eksploracja danych (data mining) to proces odkrywania znaczacych nowych powi±zañ, wzorców i trendów poprzez przeszukiwanie du¿ych ilo¶ci danych zgromadzonych w bazach danych przy u¿yciu metod matematycznych. Tu zajmiemy siê problemem klasyfikacji, czyli jak na podstawie cech pewnych obiektów przypisaæ je do okre¶lonych klas (inaczej: jak dokonaæ podzia³u na poszczególne klasy), jak równie¿ kwesti± redukcji wymiaru. W ogólno¶ci mo¿emy wyró¿niæ dwa typu klasyfikacji: pod nadzorem (supervised learning) i bez nadzoru (unsupervised learning). W pierwszym przypadku musimy dysponowaæ pewnym zbiorem danych, w którym istniej± ju¿ dane sklasyfikowane (czyli posiadaj±ce przypisan± klasê, tzw. "próba ucz±ca"). Algorytm "uczy siê" cech zwi±zanych z obiektami, a nastêpnie korzystaj±c z takiej wiedzy, klasyfikuje nowe dane. W przypadku klasyfikacji bez nadzoru, rozdzia³ na klasy nastêpuje bez uprzedniego procesu uczenia.
Prawdopodobnie najprostsz± metod± klasyfikacji pod nadzorem jest LDA (linear discriminant analysis - liniowa analiza dyskryminacyjna). U jej podstaw le¿y za³o¿enie, ¿e rozk³ady zmiennych w ka¿dej klasie s± opisane wielowymiarowym rozk³adem Gaussa. Ponadto, macierze kowariancji wszystkich klas s± takie same. Mo¿na wtedy udowodniæ, ¿e klasy najlepiej rozdziela hiperp³aszczyzna (czyli prosta dla dwóch wymiarów). Zaczynamy od prostego przyk³adu dwóch klas, ka¿da z nich jest reprezntowana przez punkty wylosowane z dwuwymiarowego rozk³adu normalnego (funkcja mvrnorm(liczba_danych, wartosc_oczekiwana, macierz_kowariancji)
z biblioteki MASS
)
# Przyk³ad 11.1
library(MASS)
library(ggplot2)
S <- matrix(c(3,0,0,3),2,2)
m1 <- c(2,2)
m2 <- c(-1,-1)
n1 <- 60
n2 <- 20
n <- n1 + n2
x1 <- mvrnorm(n1, m1, S)
x2 <- mvrnorm(n2, m2, S)
klasy <- c(rep("1", n1), rep("2", n2))
wsp1 <- c(x1[,1],x2[,1])
wsp2 <- c(x1[,2],x2[,2])
dane <- data.frame(klasy, wsp1, wsp2)
theme_set(theme_bw())
ggplot(dane) + geom_point(aes(x = wsp1, y = wsp2, color = klasy), shape = 19, size = 3)
W pakiecie R LDA jest realizowane przez funkcjê lda()
, gdzie jako argumenty podajemy zale¿no¶æ pomiêdzy polami ramki danych (np. z ~ x + y
). W efekcie otrzymamy warto¶ci prawdopodobieñstw a priori (czêsto¶ci klas), warto¶ci oczekiwane punktów nale¿±cych do klas a tak¿e kierunek, na który nale¿y rzutowaæ obserwacje tak, aby dokonaæ najlepszego rozdzielenia klas. Je¿eli warto¶æ tak dokonanego rzutowania danej (punktu) jest wiêksza od zera, przypisujemy go do klasy "1", w przeciwnym wypadku do klasy "2".
# Przyk³ad 11.2
dane.lda <- lda(klasy ~ wsp1 + wsp2, data = dane)
dane.lda
## Call:
## lda(klasy ~ wsp1 + wsp2, data = dane)
##
## Prior probabilities of groups:
## 1 2
## 0.75 0.25
##
## Group means:
## wsp1 wsp2
## 1 1.7956722 1.873495
## 2 -0.2881444 -1.081629
##
## Coefficients of linear discriminants:
## LD1
## wsp1 -0.2660515
## wsp2 -0.4867086
proj <- as.matrix(dane[,2:3]) %*% dane.lda$scaling
dane$proj <- proj[,1]
ggplot(dane) + geom_point(aes(x = proj, y = 0, color = klasy), shape=21, size=5) + geom_vline(xintercept=0)
Oczywi¶cie "najlepsze" rozdzielenie klas nie oznacza, ¿e otrzymamy podzia³, w którym wszystkie punkty z klasy "1" faktycznie zostan± sklasyfikowane jako nale¿±ce do tej grupy - widaæ to na powy¿szym rysunku. Za pomoc± biblioteki klaR
i zaimplementowanej w niej funkcji partimat()
mo¿na otrzymaæ explicite prost± rozdzielaj±c± klasy, a tak¿e szybko zweryfikowaæ, które punkty zosta³y b³êdnie sklasyfikowane. Ilo¶ciowej oceny dok³adno¶ci klasyfikacji dokonamy przy u¿yciu macierzy pomy³ek.
# Przyk³ad 11.3
library(klaR)
partimat(klasy ~ wsp1 + wsp2, data = dane, method="lda")
library(caret)
## Loading required package: lattice
dane.pred <- predict(dane.lda, dane)
dane.conf <- confusionMatrix(dane.pred$class,dane$klasy, positive = "1")
dane.conf
## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2
## 1 57 6
## 2 3 14
##
## Accuracy : 0.8875
## 95% CI : (0.7972, 0.9472)
## No Information Rate : 0.75
## P-Value [Acc > NIR] : 0.001839
##
## Kappa : 0.6842
##
## Mcnemar's Test P-Value : 0.504985
##
## Sensitivity : 0.9500
## Specificity : 0.7000
## Pos Pred Value : 0.9048
## Neg Pred Value : 0.8235
## Prevalence : 0.7500
## Detection Rate : 0.7125
## Detection Prevalence : 0.7875
## Balanced Accuracy : 0.8250
##
## 'Positive' Class : 1
##
W sytuacji gdy nie jest spe³nione za³o¿enie o równo¶ci macierzy kowariancji we wszystkich grupach, skuteczno¶æ liniowej analizy dyskryminacyjnej jest du¿o mniejsza i dlatego lepiej skorzystaæ z kwadratowej analizy dyskryminacyjnej QDA (quadratic discriminant analysis). Porównamy obie metody na sztucznym zbiorze zawieraj±cym obserwacje z trzech grup o ró¿nych macierzach kowariancji.
# Przyk³ad 11.4
S1 <- matrix(c(4,0,0,4),2,2)
S2 <- matrix(c(2,0,0,1),2,2)
S3 <- matrix(c(1,0,0,4),2,2)
m1 <- c(3,3)
m2 <- c(-1,-1)
m3 <- c(4,-2)
n1 <- 60
n2 <- 40
n3 <- 50
n <- n1 + n2 + n3
x1 <- mvrnorm(n1, m1, S1)
x2 <- mvrnorm(n2, m2, S2)
x3 <- mvrnorm(n3, m3, S3)
klasy <- c(rep("1", n1), rep("2", n2), rep("3", n3))
wsp1 <- c(x1[,1], x2[,1], x3[,1])
wsp2 <- c(x1[,2], x2[,2], x3[,2])
dane <- data.frame(klasy, wsp1, wsp2)
ggplot(dane) + geom_point(aes(x = wsp1, y = wsp2, color = klasy), shape = 19, size = 3)
partimat(klasy ~ wsp1 + wsp2, data = dane, method="lda")
partimat(klasy ~ wsp1 + wsp2, data = dane, method="qda")
# Trenujemy klasyfikator LDA i testujemy go przy u¿yciu 10-krotnej walidacji krzy¿owej
dane.lda.cv <- train(klasy ~ wsp1 + wsp2, data = dane, method="lda", trControl = trainControl(method = "cv", number=10))
# Trenujemy klasyfikator QDA i testujemy go przy u¿yciu 10-krotnej walidacji krzy¿owej
dane.qda.cv <- train(klasy ~ wsp1 + wsp2, data = dane, method="qda", trControl = trainControl(method = "cv", number=10))
dane.lda.cv
## Linear Discriminant Analysis
##
## 150 samples
## 2 predictors
## 3 classes: '1', '2', '3'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 135, 135, 135, 135, 135, 135, ...
## Resampling results:
##
## Accuracy Kappa
## 0.9066667 0.8591809
dane.qda.cv
## Quadratic Discriminant Analysis
##
## 150 samples
## 2 predictors
## 3 classes: '1', '2', '3'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 135, 135, 135, 135, 135, 135, ...
## Resampling results:
##
## Accuracy Kappa
## 0.9666667 0.9491761
# Obserwacje zrzutowane na proste ortogonalne do prostych rozdzielaj±cych
dane.lda <- lda(klasy ~ wsp1 + wsp2, data = dane)
proj <- as.matrix(dane[,2:3]) %*% dane.lda$scaling
dane.proj <- data.frame(proj, klasy=dane$klasy)
g <- ggplot(dane.proj, aes(x=LD1, y=LD2))
g + geom_point(aes(color=klasy), size=3)
Jedn± z metod klasyfikacji bez nadzoru jest analiza skupieñ (klastrów), realizowana tu przez metodê k-¶rednich (k-means, w R funkcja kmeans()
). Mamy do dyspozycji pewien zestaw danych, które chcieliby¶my pogrupowaæ tak, aby dane podobne do siebie znalazly siê w jednej klasie, natomiast ró¿ne w odrêbnych klasach. Podobieñstwo i ró¿nicê definiujemy poprzez odleg³o¶æ - w przypadku przez nas rozwa¿anym jest to odleg³o¶æ euklidesowa. Algorytm ma za zadanie tak rodzieliæ punkty, aby do jednego skupiska trafi³y te, które s± sobie najbli¿sze. Parametrem algorytmu jest liczba skupieñ, na która ma zostaæ podzielony zestaw danych.
# Przyk³ad 11.5
sigma <- matrix(c(1,0,0,1),2,2)
mu1 <- c(5,5)
mu2 <- c(1,1)
mu3 <- c(4,-2)
kolory <- c(rep("orange", 30), rep("violet", 30), rep("green", 30))
clust1 <- mvrnorm(30, mu1, sigma)
clust2 <- mvrnorm(30, mu2, sigma)
clust3 <- mvrnorm(30, mu3, sigma)
all_points <- rbind(clust1, clust2, clust3)
xrange <- range(all_points[,1])
yrange <- range(all_points[,2])
xmin = xrange[1]; xmax = xrange[2]
ymin = yrange[1]; ymax = yrange[2]
par(mfrow = c(2,2))
plot(all_points, col=kolory, pch=19, cex=2, xlab="X", ylab="Y", xlim=c(xmin,xmax), ylim=c(ymin,ymax))
title("Klastry", cex.main=1.4, font=2)
cl <- kmeans(all_points, 3)
plot(all_points, col=kolory, pch=19, cex=2, xlab="X", ylab="Y", xlim=c(xmin,xmax), ylim=c(ymin,ymax))
text(all_points[,1], all_points[,2], cl$cluster, font=2)
title("Metoda 3-srednich", cex.main=1.4, font=2)
cl1 <- kmeans(all_points, 4)
plot(all_points, col=cl1$cluster, pch=19, cex=2, xlab="X", ylab="Y", xlim=c(xmin,xmax), ylim=c(ymin,ymax))
title("Metoda 4-srednich", cex.main=1.4, font=2)
cl2 <- kmeans(all_points, 5)
plot(all_points, col=cl2$cluster, pch=19, cex=2, xlab="X", ylab="Y", xlim=c(xmin,xmax), ylim=c(ymin,ymax))
title("Metoda 5-srednich", cex.main=1.4, font=2)
Jednym z zadañ eksploracji danych jest redukcja wymiaru, czyli okre¶lenie, które ze sk³adowych wektora obserwacji s± nieistotne lub te¿ jakie inne kombinacje sk³adowych mog± siê okazaæ przydatne do dalszej analizy. Standartow± metod± redukcji wymiaru jest analiza sk³adowych g³ównych (PCA - Principal Component Analysis). Polega ona na znalezieniu nowego kierunku, który maksymalizuje wariancjê zrzutowanych na niego obserwacji. Nastêpnie szukamy kolejnego kierunku, równie¿ o jak najwiêkszej wariancji, tyle, ¿e ortogonalnego do poprzedniego itd. Okazuje siê, ¿e takie cechy odpowiadaj± wektorom w³asnym zwi±zanym z kolejnymi warto¶ciami w³asnymi (pocz±wszy od najwiêkszej). W poni¿szym przyk³adzie mamy do czynienia z "zaszumion±" relacj± y=x, alogrytm PCA (funkcja princomp()
) wykrywa jako nowe kierunki wektory [1,1] i [-1,1] (pole loadings
). Warto¶ci oryginalnych danych zrzutowane na nowe kierunki otrzymamy za pomoc± opcji scores
, natomiast przeci±¿ona funkcja plot()
prezentuje warto¶ci wariancji w kolejnych nowych kierunkach.
# Przyk³ad 11.6
x <- seq(-5, 5, by=.1)
y <- x
eta <- runif(101, max = 1)
dzeta <- runif(101, max = 1)
x <- x + eta
y <- y + dzeta
par(mfrow = c(2,2))
plot(x, y, pch=19, xlab="Test 1", ylab="Test 3", font=2, font.lab=2, xlim=c(-5,5), ylim=c(-5,5))
abline(h=0, v=0, lwd=2, col="gray")
abline(0,1,lwd=2,col="red")
abline(0,-1,lwd=2,col="green")
text(4.5,-0.5,expression(x[1]),cex=2)
text(-0.5,4.5,expression(x[2]),cex=2)
text(4.7,4,expression(y[1]),cex=2, col="red")
text(-4.5,4,expression(y[2]),cex=2, col="green")
title("Dane", cex.main=1.4)
test <- data.frame(x, y)
test.pc <- princomp(~., cor=T, data=test)
plot(test.pc, main="")
title("Wariancja", cex.main=1.4)
plot(test.pc$scores, xlim=c(-2,2), ylim=c(-2,2), xlab="Sk³adowa 1", ylab="Sk³adowa 2")
title("Sk³adowe", cex.main=1.4)
Wczytaj zbiór danych iris
z biblioteki datasets
i zapoznaj siê z nim. Przeskaluj wektory numeryczne w tej ramce danych odejmuj±c warto¶æ ¶redni± i dziel±c przez odchylenie standardowe. Wykonaj liniow± analizê dyskryminacyjn± na ca³ym zbiorze, a nastêpnie zrzutuj wszystkie obserwacje (oryginalnie czterowymiarowe) na p³aszczyznê dwuwymiarow± okre¶lon± przez proste wyznaczone w LDA. Dokonaj wizualizacji tak przetransformowanych obserwacji w pakiecie ggplot2
. Sprawd¼ skuteczno¶æ klasyfikacji wykonuj±c dowoln± walidacjê krzy¿ow±.
## Call:
## lda(Species ~ ., data = irysy.trans)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa -1.0111914 0.8504137 -1.3006301 -1.2507035
## versicolor 0.1119073 -0.6592236 0.2843712 0.1661774
## virginica 0.8992841 -0.1911901 1.0162589 1.0845261
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.6867795 0.01995817
## Sepal.Width 0.6688251 0.94344183
## Petal.Length -3.8857950 -1.64511887
## Petal.Width -2.1422387 2.16413593
##
## Proportion of trace:
## LD1 LD2
## 0.9912 0.0088
## Linear Discriminant Analysis
##
## 150 samples
## 4 predictors
## 3 classes: 'setosa', 'versicolor', 'virginica'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 135, 135, 135, 135, 135, 135, ...
## Resampling results:
##
## Accuracy Kappa
## 0.98 0.97