Eksploracja danych

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

Klasyfikacja pod nadzorem

Prawdopodobnie najprostszą metodą klasyfikacji pod nadzorem jest LDA (linear discriminant analysis - liniowa analiza dyskryminacji). 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])

dataf <- data.frame(klasy, wsp1, wsp2)

theme_set(theme_bw())

ggplot(dataf) + 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

data.lda <- lda(klasy ~ wsp1 + wsp2, data = dataf)
data.lda

## Call:
## lda(klasy ~ wsp1 + wsp2, data = dataf)
## 
## Prior probabilities of groups:
##    1    2 
## 0.75 0.25 
## 
## Group means:
##        wsp1      wsp2
## 1  2.253803  2.058642
## 2 -1.199486 -1.385207
## 
## Coefficients of linear discriminants:
##             LD1
## wsp1 -0.4044101
## wsp2 -0.4077440

proj <- as.matrix(dataf[,2:3]) %*% data.lda$scaling
dataf$proj <- proj[,1]
ggplot(dataf) + geom_point(aes(x = proj, y = 0, color = klasy), shape=21, size=5)

Oczywiście "najlepsze" rozdzielenie klas nie oznacza, że otrzymamy podział, w którym wszystkie punkty z klasy "1" faktycznie zostaną skalsyfikowane 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.

# PRZYKŁAD 11.3

library(klaR)
partimat(klasy ~ wsp1 + wsp2, data = dataf, method="lda")

Skuteczność klasyfikatora

Kluczową sprawą oprócz samej klasyfikacji jest jest skuteczność, czyli mówiąc inaczej, jak dobrze dany klasyfikator przewiduje przynależność do klas. W tym celu możemy użyć funkcji predict(), której argumentami są klasyfikator wytrenowany na próbie uczącej oraz zbiór danych. Oczywiście, w docelowych zastosowaniach skuteczność klasyfikatora jest sprawdzana na tzw. "próbie walidacyjnej", która powinna być rodzielna z próbą uczącą. Dla naszych potrzeb będziemy testować klasyfikatory LDA i QDA na próbie uczącej, otrzymując w ten sposób bardzo "optymistyczne" oszacowanie błędu. Interesować nas będzie pole class zmiennej otrzymanej z działania funkcji predict, która prównamy za pomocą polecenia table() z oryginalnymi klasami. Dla danych binarnych (tzn. dla dwóch klas) otrzymamy tzw. macierz pomyłek (confusion matrix), której pola mają nazwy jak na rysunku poniżej (numenklatura jest związana z testami medycznymi, mającymi na celu określenie nosicielstwa danej choroby).

# PRZYKŁAD 11.4

pred <- predict(data.lda, dataf)
pred$class

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [39] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 1 2 2 2 2 1 2 2 2 2 2
## [77] 1 2 2 2
## Levels: 1 2

pred.lda <- predict(data.lda, dataf)
pred.lda$class

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [39] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 1 2 2 2 2 1 2 2 2 2 2
## [77] 1 2 2 2
## Levels: 1 2

CM.lda <- table(pred.lda$class, dataf$klasy)
CM.lda

##    
##      1  2
##   1 58  3
##   2  2 17

Aby w łatwy sposób porównać poszczególne klasyfikatory (lub też ich warianty) korzysta się często z następujących trzech wielkości:

• dokładności klasyfikatora (ang. accurracy), zdefiniowanej jako \(ACC = \frac{TN + TP}{TN + TP + FN + FP}\),
• czułości klasyfikatora (ang. recall), zdefiniowanej jako \(R = \frac{TP}{TP + FN}\) oraz
• specyficzności klasyfikatora (ang. specifity), zdefiniowanej jako \(S = \frac{TN}{TN + FP}\),

Pierwsza wielkość mówi po prostu na ile dobrze algorytm przewiduje dowolną klasę. Dwie kolejne są prawdopodobieństwami warunkowymi dobrej klasyfikacji, pod warunkiem, że obiekt faktycznie należy do danej klasy. Mają one szczególne znaczenie, jeśli rozkład liczebności klas jest nierówny.

# PRZYKŁAD 11.5

cat("Accuracy LDA: ",round(sum(diag(CM.lda))/sum(CM.lda),3),"\n")

## Accuracy LDA:  0.938

cat("Recall LDA: ",round(CM.lda[1,1]/sum(CM.lda[1,]),3),"\n")

## Recall LDA:  0.951

cat("Specifity LDA: ",round(CM.lda[2,2]/sum(CM.lda[2,]),3),"\n")

## Specifity LDA:  0.895

Klasyfikacja bez nadzoru

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

sigma <- matrix(c(1,0,0,1),2,2)

mu1 <- c(4,4)
mu2 <- c(1,1)
mu3 <- c(4,-1)

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)

Redukcja wymiaru

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

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="Skladowa 1", ylab="Składowa 2")
title("Skladowe", cex.main=1.4)