Zastosowanie pakietu R w statystyce medycznej

LABORATORIUM 11

Podczas zajêæ bêd± wykorzystywane biblioteki klaR, MASS i caret.


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


Klasyfikacja pod nadzorem

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)


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


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 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)


Zadanie punktowane

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