# początek podobny jak poprzednio
library(mlbench)
data(HouseVotes84)
library(rpart)
# podzielmy dane na trenujące i testowe losowo w stosunku ok. 2:1
rhv <- runif(nrow(HouseVotes84))
hv.train <- HouseVotes84[rhv>=0.33,]
hv.test <- HouseVotes84[rhv<0.33,]
# drzewo dla domyślnych parametrów
hv.tree0 <- rpart(Class ~ ., hv.train)
# drzewo z maksymalnie rozluźnionym kryterium stopu
hv.tree1 <- rpart(Class ~ ., hv.train, minsplit=2, cp=0)
# wykaz możliwości przycięcia dla różnych wartości cp
hv.tree1$cptable
# przykładowe drzewa uzyskiwane w wyniku przycięcia
hv.tree2 <- prune(hv.tree1, cp=0.001)
hv.tree3 <- prune(hv.tree1, cp=0.005)
hv.tree4 <- prune(hv.tree1, cp=0.01)
# testowanie
err <- function(y.true, y.pred) { sum(y.pred!=y.true)/length(y.true) }
err(hv.test$Class, predict(hv.tree1, hv.test, type="c"))
err(hv.test$Class, predict(hv.tree2, hv.test, type="c"))
err(hv.test$Class, predict(hv.tree3, hv.test, type="c"))
err(hv.test$Class, predict(hv.tree4, hv.test, type="c"))
# sprawdzamy dostępność pakietu e1071 i instalujemy w razie potrzeby
if (! "e1071" %in% row.names(installed.packages()))
install.packages("e1071")
# ładujemy pakiet w celu użycia funkcji tune
library(e1071)
# poniższa linijka potrzebna do ominięcia błędu w funkcji tune
attr(hv.train, "na.action") <- na.rpart
# potrzebne jest też opakowanie na funkcję predict
# niewymagające przekazywania argumentu type
predc <- function(...) { predict(..., type="c") }
# 10-krotna walidacja krzyżowa na zbiorze trenującym
# przy domyślnych parametrach
hv.cv10 <- tune(rpart, Class ~ ., data=hv.train, na.action=na.rpart,
predict.func=predc,
tunecontrol=tune.control(sampling="cross", cross=10))
# sprawdźmy błąd uzyskany dla walidacji krzyżowej
hv.cv10
# końcowy model zbudowanego na całości przekazanych danych
hv.cv10$best.model
# to powinno być to samo co wcześniej przy bezpośrednim wywołaniu rpart
hv.tree0
# zobaczmy jak błąd na zbiorze testowym ma się do błędu z walidacji krzyżowej
err(hv.test$Class, predict(hv.cv10$best.model, hv.test, type="c"))
# strojenie parametrów
# badane wartości dla minsplit, cp, usesurrogate
par.minsplit <- c(2, 5, 10, 20, 50)
par.cp <- c(0, 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5)
par.usesurrogate <- c(0, 1, 2)
# pominiemy argument tunecontrol, gdyż i tak
# 10-krotna walidacja krzyżowa używana jest domyślnie
# uwaga -- to może trochę potrwać
hv.tune <- tune(rpart, Class ~ ., data=hv.train, na.action=na.rpart,
ranges=list(minsplit=par.minsplit, cp=par.cp,
usesurrogate=par.usesurrogate),
predict.func=predc)
# wynik strojenia
hv.tune
# najlepszy zestaw parametrów
hv.tune$best.parameters
# model dla najlepszego zestawu parametrów zbudowany
# na całym przekazanym zbiorze danych
hv.tune$best.model
# czy i ewentualnie czym to się różni od drzewa
# dla parametrów domyślnych?
summary(hv.tree0)
summary(hv.tree.best)
# w szczególności, zobaczmy, które przykłady trafiają
# do innych liści
hv.train[hv.tree0$where!=hv.tree.best$where,]
# wyniki dla wszystkich badanych zestawów parametrów
# uporządkowane wg błędu
hv.tune$performances[order(hv.tune$performances$error),]
# ponownie zbudujmy model dla najlepszego zestawu
hv.tree.best <- rpart(Class ~ ., hv.train,
control=rpart.control(hv.tune$best.parameters))
# czy wyszło to samo drzewo, które jest zachowane jako hv.tune$best.model?
hv.tree.best
# sprawdźmy błąd na zbiorze testowym
err(hv.test$Class, predict(hv.tree.best, hv.test, type="c"))