class: center, middle, title-slide .title[ # XGBoost ] .subtitle[ ## Em classificação ] .author[ ### Jaime Utria ] .institute[ ### Departamento de Estatística - UFF ] --- <!-- macro comandos matemáticos Latex --> <script type="text/x-mathjax-config"> MathJax.Hub.Config({ TeX: { Macros: { vet: ["{\\mathbf #1}",1], prodint: ["\\langle #1, #2 \\rangle",2], posto: ["{\\mathrm{posto} (#1)}",1], tr: ["{\\mathrm{tr} (#1)}",1] } } }); </script> ## XGBoost (Chen e Guestrin (2016)) - *Extreme gradient boosting* (XGBoost) é um algoritmo (também um pacote) boosting desenhado para ser eficiente, flexivel e portável em múltiplas linguagens. - Embora XGBoost forneça as mesmas opções de boosting baseado em árvores de decisão e os mesmos hiperparâmetros vistos anteriormete, também fornece algumas vantagens sobre o método tradicional de boosting como: - Regularização; - Parada antecipada; - Processamento paralelizado; - Funções de perda customizadas; - Diferentes preditores fracos. --- ## Dados de inadimplência completo (Peter Bruce & Andrew Bruce) ``` r library(xgboost) library(caret) loan_data <- read.csv("loan_data_full.csv", stringsAsFactors = TRUE) loan_data <- loan_data[,-c(1,2)] loan_data$outcome <- ordered(loan_data$outcome, levels = c('paid off', 'default')) ``` --- ``` r dim(loan_data) ``` ``` ## [1] 45342 19 ``` ``` r colnames(loan_data) ``` ``` ## [1] "loan_amnt" "term" "annual_inc" ## [4] "dti" "payment_inc_ratio" "revol_bal" ## [7] "revol_util" "purpose" "home_ownership" ## [10] "delinq_2yrs_zero" "pub_rec_zero" "open_acc" ## [13] "grade" "outcome" "emp_length" ## [16] "purpose_" "home_" "emp_len_" ## [19] "borrower_score" ``` --- ``` r # separar preditoras e resposta predictors <- data.matrix(loan_data[,-which(names(loan_data) %in% 'outcome')]) label <- as.numeric(loan_data$outcome)-1 # separar treino/teste set.seed(12345) train_index <- createDataPartition(y = loan_data$outcome, p =0.8, list = F) X_train <- predictors[train_index,] X_test <- predictors[-train_index,] y_train <- label[train_index] y_test <- label[-train_index] # treinar xgboost xgb <- xgboost(data=X_train, label=y_train, objective='binary:logistic', nrounds=250, verbose=0, eval_metric='error') ``` --- ``` r # fazer predições no conjunto teste pred <- predict(xgb, X_test) # calcular o erro de predição error <- abs(y_test - pred) > 0.5 # erro nos dados de treino xgb$evaluation_log[250,] ``` ``` ## iter train_error ## <num> <num> ## 1: 250 0.1307548 ``` ``` r # erro nos dados de teste mean(error) ``` ``` ## [1] 0.3493604 ``` --- - Pelos resultados observados, nota-se que o erro de treino foi ao redor de 13%, ao passo que o erro de teste foi ao redor de 35%, o que sugere a presença de sobreajuste (*overfitting*). - Para contornar o problema de sobreajuste podemos fazer o ajuste dos hiperparâmetros, por exemplo, fixando a taxa de aprendizado em um valor "pequeno". - Outra alternativa é aplicar técnicas de **regularização**, isto é, modificamos a função de perda a fim de **penalizar** a complexidade do modelo. - Em XGBoost é possível modificar a função de perda adicionando um termo que "mede" a complexidade do modelo. --- ## Regularização L2 ``` r xgb_penalty <- xgboost(data = X_train,label = y_train, params = list(eta = 0.1, subsample = 0.65, lambda = 1000), objective = "binary:logistic",nrounds = 250,verbose = 0, eval_metric='error') pred_penalty <- predict(xgb_penalty,X_test) error_penalty <- abs(y_test - pred_penalty) > 0.5 xgb_penalty$evaluation_log[250,]; mean(error_penalty) ``` ``` ## iter train_error ## <num> <num> ## 1: 250 0.3125103 ``` ``` ## [1] 0.3248787 ``` --- ``` r error <- rep(0,250) error_penalty <- rep(0,250) for (i in 1:250) { pred <- predict(xgb,X_test,iterationrange = c(1,i)) error[i] <- mean(abs(y_test - pred) > 0.5) pred_penalty <- predict(xgb_penalty, X_test, iterationrange = c(1,i)) error_penalty[i] <- mean(abs(y_test - pred_penalty) > 0.5) } errors <- rbind(xgb$evaluation_log, xgb_penalty$evaluation_log, data.frame(iter = 1:250, train_error = error), data.frame(iter = 1:250, train_error = error_penalty)) errors$tipo <- rep(c("sem penalty train", "penalty train", "sem penalty test", "penalty test"),rep(250,4)) ``` --- <img src="data:image/png;base64,#Aula19_XGBoost_files/figure-html/unnamed-chunk-7-1.svg" style="display: block; margin: auto;" /> --- ## XGBoost para classificação: Ideias Principais - A **impureza** de cada nó é substituida pelo *Score* de similaridade (SS): `$$SS(a) = \frac{\left(\sum_{i: \mathbf{x}_i \in a} r_{ib}\right)^2}{\sum_{i: \mathbf{x}_i \in a} \hat{p}_i (1-\hat{p}_i) + \lambda}$$` - `\(r_{ib}\)` é o i-ésimo pseudo-resíduo obtido quando usamos o preditor `\(\hat f_b\)`; - `\(\hat p_i\)` probabilidade prevista da classe positiva para amostra `\(i\)` quando usamos o preditor `\(\hat f_{b-1}\)`; - `\(\lambda \ge 0\)` é um parâmetro de regularização. Definimos, um novo ganho de informação na divisão `\(d\)` baseado em `\(SS\)` por: `$$IG^{(SS)}(a_p,d)= SS(a_e) + SS(a_d) - SS(a_p).$$` --- ## Algoritmo XGBoost (Classificação) 1. Iniciamos `\(\hat f_0 = \log(1) \iff \hat p_i^{(0)} = 0,5\)` para todo `\(i=1,\ldots,n\)`. 2. Para `\(b = 1, \ldots, B\)`: (a) Calculamos os resíduos `\(r_{ib} = y_i - \hat p_i^{(b-1)} = y_i - \hat p_i\)` (b) Ajustamos uma árvore de regressão aos resíduos `\(\{r_{ib}\}_{i=1}^n\)` escolhendo aquela divisão que tenha o maior ganho de informação e geramos `\(J_b\)` regiões terminais. (c) O valor constante de cada nó folha é dado por: `\(\gamma_{jb} = \frac{\sum_{i: \mathbf{x}_i \in R_{jb}} r_{ib}}{\sum_{i: \mathbf{x}_i \in R_{jb}} \hat{p}_i (1-\hat{p}_i) + \color{blue}{\lambda}}.\)` (d) Atualizamos `\(\hat f_b(\mathbf x) = \hat f_{b-1} + \eta \sum_{j} \gamma_{jb} I(\mathbf x \in R_{jb})\)` 3. Preditor final: `\(\hat f_{\mathrm{XGB}}(\mathbf{x}) = \hat f_B(\mathbf{x}).\)` --- ## Poda - A poda pode ser realizada durante o treinamento de cada árvore do *ensemble*, inserindo um parâmetro `\(\gamma\)` na função de perda que controla o tamanho da árvore. - Outra abordagem é deixar crescer a árvore até um certo critério de parada, fixando o número de resíduos em cada nó folha (no contexto de XGBoost é denominado de *cover*). O parâmetro de complexidade `\(\gamma\)` é então usado como segue: `$$IG^{SS}(a) - \gamma = \begin{cases} > 0, &\text{não podamos no} a.\\ < 0, & \text{podamos no } a\end{cases}$$` - `\(\gamma\)` é interpretado como o ganho mínimo necessário para manter o ramo com `\(a\)` como raiz. - Vale a pena destacar que `\(\gamma = 0\)`, não necessariamente implica em poda no nó `\(a\)`, pois ainda podemos ter um ganho negativo. --- ## XGBoost: Detalhes **Função de perda regularizada** (https://arxiv.org/pdf/1603.02754) `$$\mathcal{L} = \sum_{i=1}^n L(y_i,\hat y_i) + \sum_{b = 1}^B \Omega(f_b),$$` em que `\(\Omega(f) = \gamma T + \frac{1}{2}\lambda||w||^2\)`, `\(w = (w_1, \ldots,w_T)\)`. - Como antes, `\(L\)` é a perda (custo) de mal classificação de `\(y_i\)` por `\(\hat y_i\)`. - `\(\Omega\)` é um termo de penalização. Via de regra: - quanto maior `\(\lambda\)` maior a penalização nos *pesos* das folhas `\(w\)` o que previne o sobreajuste. - quanto maior `\(\gamma\)` maior a penalização no tamanho da árvore, pois aqui `\(T_b\)` representa o número de folhas. --- - Note que `\(w\)` depende de `\(\hat f_b\)`, entretanto omiteremos sua dependência. - Considere `\(\gamma = 0\)`. Logo, na iteração `\(b\)`: `$$\begin{align}\mathcal{L}^{(b)} &= \sum_{i=1}^n L(y_i,\hat y_i^{(b-1)} + \hat f_{b}(\mathbf{x}_i)) + \frac{1}{2}\lambda||w||^2\\ &=\sum_{i=1}^n L(y_i,\hat y_i^{(b-1)} + \hat f_{b}(\mathbf{x}_i)) + \frac{1}{2}\lambda \sum_{j=1}^{J_b} w_j^2\end{align}$$` - Lembre que em clasificação, usamos: `$$\begin{align}L(y_i,\hat p_i) &= - [y_i \log(\hat p_i) + (1-y_i)\log (1-\hat p_i)]\\ &=-y_i\log(\hat q) + \log(1+\exp(\log(\hat q_i)))\\ &=L(y_i,\log(\hat q))\end{align}$$` --- - Como antes, usamos a aproximação de Taylor de segunda ordem: `$$\mathcal{L}^{(b)} \approx \sum_{i=1}^n [L(y_i,\hat f_{b-1}(\mathbf x_i)) + g_i \hat f_{b}(\mathbf x_i) + \frac{1}{2} h_i\hat{f}_b^2(\mathbf{x}_i)] + \frac{1}{2}\lambda \sum_{j=1}^{J_b} w_j^2,$$` em que `\(g_i = \frac{d}{d\hat f_{b-1}(\mathbf x_i)}L(\hat y_i, \hat f_{b-1}(\mathbf x_i))\)` e `\(h_i = \frac{d^2}{d^2\hat f_{b-1}(\mathbf x_i)}L(\hat y_i, \hat f_{b-1}(\mathbf x_i))\)` - Além disso, `\(\hat{f}_{b}(\mathbf x_i)\)` é a log chance prevista para a amostra `\(\mathbf{x}_i\)` pelo previsor `\(\hat f_b\)`. - Eliminando o termo constante em `\(\hat f_b(\mathbf x_i)\)` (que não depende desse valor), resulta em: `$$\begin{align}\tilde{\mathcal{L}}^{(b)} &\approx \sum_{i=1}^n [g_i \hat f_{b}(\mathbf x_i) + \frac{1}{2} h_i\hat{f}_b^2(\mathbf{x}_i)] + \frac{1}{2}\lambda \sum_{j=1}^{J_b} w_j^2\\ &=\sum_{j=1}^{J_b} \left[\left(\sum_{i: \mathbf x_i \in R_{jb}} g_i\right) w_j + \frac{1}{2}\left(\sum_{i: \mathbf x_i \in R_{jb}} h_i + \lambda\right)w_j^2\right]\end{align}$$` --- - Minimizando em `\(w_j \equiv w_{jb}\)` encontramos (após "ALGUNS" cálculos) que o valor ótimo é: `$$\begin{align}\gamma_{jb} &= -\frac{\sum_{i: \mathbf{x}_i \in R_{jb}} g_i}{\sum_{i: \mathbf{x}_i \in R_{jb}} h_i + \color{blue}{\lambda}}\\ &=\frac{\sum_{i: \mathbf{x}_i \in R_{jb}} r_{ib}}{\sum_{i: \mathbf{x}_i \in R_{jb}} \hat{p}_i (1-\hat{p}_i) + \color{blue}{\lambda}}.\end{align}$$` - Logo, a **impureza** (*similarity score*) para a árvore construída na iteração `\(b\)` é: `$$\begin{align}\mathcal L^{(b)}(T_b) = -\color{red}{\frac{1}{2}}\sum_{j=1}^{J_b}\frac{\left(\sum_{i: \mathbf{x}_i \in R_{jb}} g_i\right)^2}{\sum_{i: \mathbf{x}_i \in R_{jb}} h_i + \color{blue}{\lambda}} &\overset{*}{=}\sum_{j=1}^{J_b}\frac{\left(\sum_{i: \mathbf{x}_i \in R_{jb}} r_{ib}\right)^2}{\sum_{i: \mathbf{x}_i \in R_{jb}} \hat{p}_i (1-\hat{p}_i) + \color{blue}{\lambda}}\\ &= \sum_{j=1}^{J_b} SS(\tau_j).\end{align}$$` --- ## Ajuste de hiperparâmetros ``` r N <- nrow(loan_data) set.seed(12345) fold_number <- sample(1:10, N, replace=TRUE) # Grid de parâmetros para testar # 9 configuraçoes params <- data.frame(eta = rep(c(0.1, 0.5, 0.9), 3), max_depth = rep(c(3, 6, 12), rep(3,3))) head(params) ``` ``` ## eta max_depth ## 1 0.1 3 ## 2 0.5 3 ## 3 0.9 3 ## 4 0.1 6 ## 5 0.5 6 ## 6 0.9 6 ``` --- ``` r rf_list <- vector('list', 9) error <- matrix(0, nrow=9, ncol=5) for(i in 1:nrow(params)){ for(k in 1:5){ cat('Fold', k, 'for model', i, '\n') fold_idx <- (1:N)[fold_number == k] xgb <- xgboost(data=predictors[-fold_idx,], label=label[-fold_idx], params=list(eta=params[i, 'eta'], max_depth=params[i, 'max_depth']), objective='binary:logistic', nrounds=100, verbose=0, eval_metric='error') pred <- predict(xgb, predictors[fold_idx,]) error[i, k] <- mean(abs(label[fold_idx] - pred) >= 0.5) } } ``` ``` ## Fold 1 for model 1 ## Fold 2 for model 1 ## Fold 3 for model 1 ## Fold 4 for model 1 ## Fold 5 for model 1 ## Fold 1 for model 2 ## Fold 2 for model 2 ## Fold 3 for model 2 ## Fold 4 for model 2 ## Fold 5 for model 2 ## Fold 1 for model 3 ## Fold 2 for model 3 ## Fold 3 for model 3 ## Fold 4 for model 3 ## Fold 5 for model 3 ## Fold 1 for model 4 ## Fold 2 for model 4 ## Fold 3 for model 4 ## Fold 4 for model 4 ## Fold 5 for model 4 ## Fold 1 for model 5 ## Fold 2 for model 5 ## Fold 3 for model 5 ## Fold 4 for model 5 ## Fold 5 for model 5 ## Fold 1 for model 6 ## Fold 2 for model 6 ## Fold 3 for model 6 ## Fold 4 for model 6 ## Fold 5 for model 6 ## Fold 1 for model 7 ## Fold 2 for model 7 ## Fold 3 for model 7 ## Fold 4 for model 7 ## Fold 5 for model 7 ## Fold 1 for model 8 ## Fold 2 for model 8 ## Fold 3 for model 8 ## Fold 4 for model 8 ## Fold 5 for model 8 ## Fold 1 for model 9 ## Fold 2 for model 9 ## Fold 3 for model 9 ## Fold 4 for model 9 ## Fold 5 for model 9 ``` ``` r avg_error <- 100 * round(rowMeans(error), 4) ``` --- ## Melhor configuração (modelo) ``` r results <- data.frame(params, avg_error) results[which.min(results$avg_error),] ``` ``` ## eta max_depth avg_error ## 1 0.1 3 32.8 ``` --- ``` r xgb_melhor <- xgboost(data=X_train, label=y_train, params=list(eta=0.1, max_depth=3), objective='binary:logistic', nrounds=100, verbose=0, eval_metric='error') pred <- predict(xgb_melhor, X_test) error <- abs(y_test - pred) > 0.5 xgb_melhor$evaluation_log[100,]; mean(error) ``` ``` ## iter train_error ## <num> <num> ## 1: 100 0.320698 ``` ``` ## [1] 0.3298412 ```