Drzewa
MŁ
Podstawy
Pakiet do budowy drzewa:
library(rpart)Pakiety do rysowania wykresów:
library(rpart.plot)
library(plotmo)
# library(tree) - wymaga nowszego R niż mam, więc tylko zaznaczam, że istnieje taki pakietDane:
The kyphosis data frame has 81 rows and 4 columns. representing data on children who have had corrective spinal surgery
Kyphosis - a factor with levels absent present indicating if a kyphosis (a type of deformation) was present after the operation.
Age - in months
Start - the number of the first (topmost) vertebra operated on.
df <- kyphosis[, -3]
head(df)## Kyphosis Age Start
## 1 absent 71 5
## 2 absent 158 14
## 3 present 128 5
## 4 absent 2 1
## 5 absent 1 15
## 6 absent 1 16Budowanie drzewa polega na tym, że dokonujemy np. takich podziałów:
wszystkie obserwacje czarne klasyfikujemy jako czarne jeśli Age < 50, pozostałe czerwone
wszystkie obserwacje czarne klasyfikujemy jako czarne jeśli Age < 50 i Start > 7
itd. (w szególności dalej dzielimy Age >= 50, Age < 50 & Start <= 7)
Ale zmienne po których dzielimy i punkt podziału dobieramy na podstawie miary różnorodności klas w węźle i różnicy pomiędzy różnorodnością w węźle rodzicu i węzłach dzieciach.
png('wstepny.png')
par(mfrow=c(1,3))
plot(Start~Age, data=df, col=Kyphosis, pch=16)
plot(Start~Age, data=df, col=Kyphosis, pch=16)
abline(v=50)
plot(Start~Age, data=df, col=Kyphosis, pch=16)
abline(v=50)
lines(c(-2, 50), c(7, 7))
dev.off()
Budowa drzewa
Budujemy drzewo:
set.seed(123)
tree <- rpart(Kyphosis ~ ., data = df)A dzieje się to w takiej kolejności:
png('budowa.png')
par(mfrow = c(3,3))
for(iframe in 1:nrow(tree$frame)) {
cols <- ifelse(1:nrow(tree$frame) <= iframe, "black", "gray")
prp(tree, col = cols, branch.col = cols, split.col = cols)}
dev.off()
Najpierw jesteśmy na górze. Niech \(Q_p\) będzie miarą różnorodności rodzica, a \(Q_{dl}\), \(Q_{dp}\) dziecka prawego i lewego. Chcemy wybrać taki podział, żeby różnica \(Q_p - \hat{p}_{l}Q_{dl} - \hat{p}_{p}Q_{dp}\) była maksymalna (\(\hat{p}_{l}\), \(\hat{p}_{p}\) - estymowane p-stwa bycia lewym i prawym dzieckiem odpowiednio). Jak dokonamy podziału, to schodzimy poziom niżej i od nowa.
Domyślnie \(Q\) to parms = list(split = "gini"), można zmienić na entropię tak: parms = list(split = "information").
Kolejne podziały w drzewie są zapisane tak:
tree$frame## var n wt dev yval complexity ncompete nsurrogate yval2.V1 yval2.V2
## 1 Start 81 81 17 1 0.17647059 1 0 1.00000000 64.00000000
## 2 Start 62 62 6 1 0.01960784 1 1 1.00000000 56.00000000
## 4 <leaf> 29 29 0 1 0.01000000 0 0 1.00000000 29.00000000
## 5 Age 33 33 6 1 0.01960784 1 1 1.00000000 27.00000000
## 10 <leaf> 12 12 0 1 0.01000000 0 0 1.00000000 12.00000000
## 11 Age 21 21 6 1 0.01960784 1 0 1.00000000 15.00000000
## 22 <leaf> 14 14 2 1 0.01000000 0 0 1.00000000 12.00000000
## 23 <leaf> 7 7 3 2 0.01000000 0 0 2.00000000 3.00000000
## 3 <leaf> 19 19 8 2 0.01000000 0 0 2.00000000 8.00000000
## yval2.V3 yval2.V4 yval2.V5 yval2.nodeprob
## 1 17.00000000 0.79012346 0.20987654 1.00000000
## 2 6.00000000 0.90322581 0.09677419 0.76543210
## 4 0.00000000 1.00000000 0.00000000 0.35802469
## 5 6.00000000 0.81818182 0.18181818 0.40740741
## 10 0.00000000 1.00000000 0.00000000 0.14814815
## 11 6.00000000 0.71428571 0.28571429 0.25925926
## 22 2.00000000 0.85714286 0.14285714 0.17283951
## 23 4.00000000 0.42857143 0.57142857 0.08641975
## 3 11.00000000 0.42105263 0.57894737 0.23456790Można narysować to drzewo:
png("ubogie.png")
par(mar=c(2,2,2,2))
par(xpd = NA)
plot(tree)
text(tree, use.n = TRUE)
dev.off()
Na tej podstawie możemy dokonywać klasyfikacji i szacować prawdopodobieństwo przynależności do klasy.
Np. obserwacja, dla której 14.5 > Start >= 8.5, 111 > Age >=55, znajduje się w wężle, w którym mamy 3 obserwacje o klasie “absent” i 4 obserwacje “present”, więc predykcja klasy dla nowej obserwacji to “present”, a szacowane prawdopodobieństwo, że obserwacja jest w klasie “present” to 4/7.
predict(tree, data.frame(Start=10, Age = 90))## absent present
## 1 0.4285714 0.5714286predict(tree, data.frame(Start=10, Age = 90), type="class")## 1
## present
## Levels: absent presentInne wykresy drzewa
png("mniej_ubogie.png")
prp(tree, type=1, extra=106)
dev.off()
png("3d_wykres.png")
plotmo(tree, type = "prob", nresponse = "present")
dev.off()
png("2d_wykres.png")
plotmo(tree, type = "prob", nresponse = "present", type2 = "image", ngrid2 = 200, pt.col = ifelse(kyphosis$Kyphosis == "present", "red", "lightblue"))
dev.off()
Przycinanie drzewa
Drzewo bez przycinania - budujemy je do końca, tzn. aż w węzłach będą obserwacje tylko z jednej klasy.
tree <- rpart(Kyphosis ~ . - Number, data = kyphosis, control = rpart.control(cp = -1, minsplit=2, minbucket = 1))
png("pelne.png")
par(mfrow=c(1, 2))
prp(tree, extra=106)
par(mar=c(1,1,1,1))
par(xpd = NA)
plot(tree)
text(tree, use.n = TRUE)
dev.off()
Takie drzewo jest przeuczone i może nie działać na zbiorze testowym, dlatego ważne jest, by je przyciąć.
Kryterium kosztu-złożoności
\[R_{\alpha}(T) = R(T) + \alpha |T|\]
\(R_{\alpha}(T)\) - miara niedoskonałości drzewa, \(|T|\) - rozmiar drzewa (liczba liści), \(\alpha\) - współczynnik złożoności (cp - complexity parameter).
png("plot_cp_tree.png")
plotcp(tree)
dev.off()
Dla naszych danych słabo wychodzi, niestety.
To jest wykres zależności miary niedopasowania drzewa (relative error, tzn. \(1 – R^2\) w przypadku regresyjnym, a klasyfikacyjnym ułamek błędnych klasyfikacji) od wartości cp. Na górze wypisana jest liczba zmiennych.
To jest mało edukacyjny rysunek, zwykle te wykresy wyglądają raczej tak:
Bardziej edukacyjny przykład :) Tutaj widać, że minimalny błąd jest dla cp=0.048 i od miejsca, gdzie kończy się jego odchylenie jest narysowana pozioma przerywana linia, więc stosując zasadę 1SE możemy wybrać cp=0.14.
tree_przyciete <- prune.rpart(tree, cp=0.068)
tree_przyciete2 <- rpart(Kyphosis ~ ., data = df, control = rpart.control(cp = 0.068))Inne możliwości ograniczania rozrostu drzewa można znaleźć tak: ?rpart.control (są dobrze opisane).