Dane

Skorzystamy z danych Boston. Zmienna medv to mediana wartości domu i to będziemy modelować w zależności od innych zmiennych (przekodowałam na 1 - drogie domy, 0 - nie). Od razu robię podział na trening i zbiór testowy.

library(gbm)
library(MASS)
df <- Boston
df$medv <- ifelse(df$medv < 25, 0, 1)
str(df)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  0 0 1 1 1 1 0 1 0 0 ...
set.seed(123)
ind.train <- sample(1:nrow(Boston), size=350, replace=FALSE)
ind.test  <- setdiff(1:nrow(Boston), ind.train)

Bagging, boosting

Pakiet adabag

bagging

mfinal - liczba drzew

parametry drzew można kontrolować dodając coś takiego: control = rpart.control(minsplit = 2, cp = 0).

library(adabag)
df$medv <- factor(df$medv)
model.bagging <- bagging(medv ~ ., data = df[ind.train, ], mfinal = 100)
predykcje.drzewa <- lapply(model.bagging$trees, predict, newdata=df[ind.test, ])
predykcje.drzewa <- data.frame(lapply(predykcje.drzewa, function(x) x[,2]))

predykcje.bagging <- predict(model.bagging, df[ind.test, ])$prob[,2]

kolejnosc <- order(predykcje.bagging)

matplot(predykcje.drzewa[kolejnosc,], pch=16, col="grey", 
        xlab="obserwacje ze zbioru testowego (posortowane wg predykcji z baggingu)",
        ylab="response")
points(as.numeric(as.character(df[ind.test, "medv"]))[kolejnosc], col="red", pch=16)
points(predykcje.bagging[kolejnosc])

model.bagging$importance
##        age      black       chas       crim        dis      indus      lstat 
##  0.7215210  0.2666298  0.0000000  2.2685897  3.6527071  2.7488380 26.1625597 
##        nox    ptratio        rad         rm        tax         zn 
##  0.6171060  1.0015116  0.8123735 58.8123389  2.8737960  0.0620287
predykcje.bagging <- sapply(c(1, 10, 50, 100), 
                            function(i) predict(model.bagging, 
                                                newdata=df[ind.test,], newmfinal=i)$prob[,2])

matplot(predykcje.bagging[kolejnosc,], type="l", col=c("yellow", "orange", "red", "darkred"), lty=1, ylab="response")
points(as.numeric(as.character(df[ind.test, "medv"]))[kolejnosc], col="black", pch=16)
legend("topleft", col=c("yellow", "orange", "red", "darkred"), lty=c(1,1,1,1),
       legend=c(1, 10, 50, 100), title="Liczba drzew")

boosting

model.boosting.adabag <- boosting(medv ~ ., data = df[ind.train, ], mfinal = 50)

Gradient Boosting

Pakiet gbm

df$medv <- as.numeric(as.character(df$medv))
model.boosting <- gbm(medv ~ . , data = df[ind.train, ])
## Distribution not specified, assuming bernoulli ...

Jakie parametry można pozmieniać

Po pierwsze, można zmienić rozkład, a dokładniej jaka funkcja straty będzie optymalizowana.

model.boosting.logist <- gbm(medv ~ . , data = df[ind.train, ], distribution = "bernoulli")
model.boosting.ada <- gbm(medv ~ . , data = df[ind.train, ], distribution = "adaboost")

Po drugie, liczbę drzew.

model.boosting.wiecej.drzew <- gbm(medv ~ . , data = df[ind.train, ],
                                   n.trees = 500)
## Distribution not specified, assuming bernoulli ...

I sporo innych wg potrzeb, np. bag.fraction odnosi się do tego, jak liczna ma być grupa obserwacji, na której uczymy drzewo kolejnej generacji (domuślnie 0.5, czyli połowa wejściowego zbioru).

Jakie informacje można wyciągnąć

Można dowiedzieć się, jakie zmienne są ważne w modelu (jak dokładnie liczone jest rel.inf możecie znaleźć tutaj: ?summary.gbm w details, a dokładniejszy opis w artykule).

summary(model.boosting)

##             var    rel.inf
## rm           rm 42.0987793
## lstat     lstat 40.5420137
## tax         tax  5.3837755
## dis         dis  4.6401381
## ptratio ptratio  1.6294108
## black     black  1.4379879
## age         age  1.2280230
## crim       crim  1.0264490
## nox         nox  0.7907013
## indus     indus  0.6997543
## rad         rad  0.3606352
## chas       chas  0.1623318
## zn           zn  0.0000000

Można dowiedzieć się, jak konkretna zmienna wpływa na zmienną odpowiedzi (czyli tak brzegowo, wszystkie pozostałe obserwacje są ustalone i ruszamy tą jedną zmienną i patrzymy, jaki to ma wpływ na zmienną objaśnianą). Poniżej przykład, że to ma sens (nie szukałam specjalnie, wzięłam najistotniejszą zmienną po prostu). rm oznacza średnią liczbę pokoi. Widać, że gdy liczba pokoi wzrasta, to prawdopodobieństwo, że mamy do czynienia z drogim domem też, i widać też, że istotny skok następujemy przy przejściu z 6 na 7.

plot(model.boosting, i.var="rm", type="response")

Można się dowiedzieć, ile generacji drzew potrzeba. predykcje daje macierz predykcjii w kolumnach jest pstwo 1 dla kolejno 10, 20,…, 100 generacji drzew. Na tej podstawie można stwierdzić, że może warto zakończyć naukę na 80 tzn. ustawić n.trees=80?

predykcje <- predict(model.boosting, df[ind.test,], 
                     n.trees=seq(10, 100, 10), type="response")
y.test <- df[ind.test, "medv"]
logwiarogodnosc <- apply(predykcje, 2, function(x){
  sum(y.test*log(x) + (1 - y.test)*log(1-x))
})

plot(seq(10, 100, 10), logwiarogodnosc, xlab="n.trees")

Pakiet xgboost

Pakiet ma inną składnię. Chociaż podaję mniej przykładów, to zachęcam do stosowania. Twórcy zapewniają, że dobry i szybszy niż gbm.

library(xgboost)

model.xgboost <- xgboost(data=as.matrix(df[ind.train, -14]), label=df[ind.train,"medv"],
                         nrounds=5, objective="binary:logistic")
## [1]  train-error:0.048571 
## [2]  train-error:0.045714 
## [3]  train-error:0.042857 
## [4]  train-error:0.034286 
## [5]  train-error:0.022857
importance_matrix <- xgb.importance(model = model.xgboost)
xgb.plot.importance(importance_matrix, top_n = 10, measure = "Gain")

Lasy losowe

Pakiet randomForest

library(randomForest)

df$medv <- factor(df$medv)

model.randomForest <- randomForest(medv ~ . , data = df[ind.train, ])
model.randomForest$importance
##         MeanDecreaseGini
## crim           5.4948615
## zn             3.7247963
## indus         11.4027399
## chas           0.6606894
## nox            6.6808172
## rm            39.4590815
## age            5.6325015
## dis            7.7673576
## rad            2.7004767
## tax            7.8644570
## ptratio        7.1287258
## black          3.8591295
## lstat         30.4665807
varImpPlot(model.randomForest)

plot(model.randomForest)

par(mfrow=c(1,2))
partialPlot(model.randomForest, df[ind.train, ], x.var="rm", which.class=1)
partialPlot(model.randomForest, df[ind.test, ], x.var="rm", which.class=1)

trzecie.drzewo <- getTree(model.randomForest, k=3, labelVar=TRUE)

Pakiet ranger

Działa szybciej niż randomForest. Wiele paramtrów trzeba odpowiednio ustawić, żeby wyniki zostały zapisane, poniżej przykład z importance ("impurity" związane jest z miarą Giniego).

library(ranger)

model.ranger <- ranger(medv ~ . , data = df[ind.train, ], importance="impurity")
model.ranger$variable.importance
##       crim         zn      indus       chas        nox         rm        age 
##  5.8509189  3.1894777 10.6808521  0.6332659  7.7041253 36.6824968  5.9417027 
##        dis        rad        tax    ptratio      black      lstat 
##  7.7038017  2.6346848  8.9874327  6.8268225  4.1166901 31.5946619

Inne

Są jeszcze np. randomForestSRC, Rborist, Random Jungle.

Polecam

Bardziej przeglądałam niż uważnie czytałam, ale wyglądają sensownie:

http://uc-r.github.io/gbm_regression#xgboost