XGBoost

class: center, middle, title-slide

.title[
# XGBoost
]
.subtitle[
## Em regressão
]
.author[
### Jaime Utria
]
.institute[
### Departamento de Estatística - UFF
]

---

<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

- XGBoost (*eXtreme Gradient Boosting*) é um algoritmo da família boosting que estende o tradicional gradiente boosting.

- As principais vantagens do XGBoost incluem técnicas de:

- Regularização;
  - Parada antecipada;
  - Processamento paralelizado;
  - Funções de perda customizadas;
  - Diferentes preditores fracos.

---

## Escore de similaridade para regressão

Seja `$a$` um nó da árvore `$T_b$`, definimos o escore de similaridade por:

`$$SS_{reg}(a) = \frac{\left(\sum_{i : \mathbf{x}_i \in a} y_i - \hat y_i^{(b-1)}\right)^2}{n_a + \lambda} = \frac{(\text{soma de resíduos em }a)^2}{\text{# de resíduos em }a + \lambda}$$`
- `$\hat y_i^{(b-1)}$` é o valor predito na iteração `$b-1$`.

- `$n_a$` é o número de resíduos no nó `$a$`.

- Como antes, `$\lambda \geq 0$` é um parâmtero de regularização.

**OBS**. Note que `$SS_{reg}(a)$` grande significa que os resíduos em `$a$` são *similares*, ao passo que `$SS_{reg}(a)$` pequeno significa que os resíduos em `$a$` são *dissimilares*.

---

## Árvore XGBoost

Definimos, o ganho de informação na divisão `$d$` do nó pai `$a_p$` por:

`$$IG^{(SS_{reg})}(a_p,d) = SS_{reg}(a_e) + SS_{reg}(a_d) - SS_{reg}(a_p)$$`

Uma *árvore de regressão XGBoost* é simplesmente uma árvore de regressão construída usando como critério de divisão o índice `$IG^{(SS_{reg})}$`, ao inves do índice `$IG$` usado na construção das árvores de regressão. tradicionais.

---

## Poda

- A poda de uma árvore XGBoost é realizada com base nos ganhos de informação.

- Introduzimos o parâmetro de complexidade `$\gamma \geq 0$`, e podamos o ramo com raiz no nó `$a$`, de acordo a:

`$$IG^{(SS_{reg})}(a) - \gamma = \begin{cases} > 0, & \text{não podamos no }a,\\<0, & \text{podamos no }a.\end{cases}$$`
- Como antes, `$\gamma = 0$`, não necessariamente implica em poda no nó `$a$`, pois ainda podemos ter um ganho negativo.

---

## Algoritmo XGBoost (Regressão)

1. Iniciamos `$\hat f_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 y_i^{(b-1)}$`.
  
  (b) Ajustamos uma árvore de regressão *XGBoost* aos resíduos `$\{r_{ib}\}_{i=1}^n$` e geramos `$R_{jb}, j = 1, \ldots, 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}}{|R_{jb}| + \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}).$`

---

## Sobre o parâmetro de regularização

- O parâmetro `$\lambda$` reduz a sensibilidade de previsões individuais. 
`$$\gamma_{jb} = \frac{\sum_{i: \mathbf{x}_i \in R_{jb}} r_{ib}}{|R_{jb}| + \color{blue}{\lambda}}.$$`
    - Se `$\lambda = 0$`, então `$\gamma_{jb}$` é simplesmente a média dos resíduos em `$R_{jb}$`. Consequentemente, se temos um único resíduo na região `$R_{jb}$`, então `$\gamma_{jb}$` será aquele resíduo.
  
    - Porém, se `$\lambda > 0$`, então a previsão em `$R_{jb}$`, `$\gamma_{jb}$`, será reduzida (obtendo uma estimativa mais *suavizada*).

- Digamos, por exemplo, que `$R_{jb}$` tem só o resíduo igual a 10. Se `$\lambda = 0$`, então `$\gamma_{jb} = 10$` e se tomamos `$\lambda = 1$`, então `$\gamma_{jb} = 5$` (a contribuição dessa folha à previsão final foi reduzida em 50%).

---

- O parâmetro `$\lambda$` reduz a complexidade do modelo (prevenindo o sobreajuste)

`$$SS_{reg}(a) = \frac{\left(\sum_{i: \mathbf{x}_i \in a} r_{ib}\right)^2}{n_a + \color{blue}{\lambda}},$$`
em que `$a$` é um nó interno de `$T_b$`.

- Quando `$\lambda > 0$` o escore de similaridade de `$a$` é reduzido o que leva a um ganho menor ao dividí-lo, portanto tende a ser podado.

- Considere por exemplo, 
  
  `$IG^{(SS)}(a) = 140$` quando `$\lambda = 0$` e 
  
  `$IG^{(SS)}(a) = 80$` quando `$\lambda = 1$`. 
  
  Tomando, por exemplo, `$\gamma = 120$`, então sem regularização não podamos `$a$` e com regularização podamos `$a$`.  
  
---

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

- O vetor `$w$` contém os valores constantes (saída/previsões) em cada nó folha da árvore.  
  
---

## Valores de saída e escores de similaridade

- Seja `$\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_b||^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_{jb}^2\end{align}$$`

- Os valores de saída `$w_b$` são obtidos minimizando a função de perda `$\mathcal{L}^{(b)}.$`

---

- Como antes, usamos a aproximação de Taylor de segunda ordem:

`$$\mathcal{L}^{(b)} = \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_{jb}^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))$`

- 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)} &= \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_{jb}^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_{jb}^2\right]\end{align}$$`

---

- Finalmente, minimizando em `$w_{jb}$`, obtemos o valor ótimo:

`$$\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}}$$`

- A **impureza** ótima da árvore `$T_b$` é:

`$$\tilde{\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}},$$`

- O escore de similaridade do `$j$`-ésimo nó terminal, `$\tau_{jb}$`, de `$T_b$` é:

`$$SS(\tau_{jb}) = -\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}}$$`
---

- Lembre que em regressão, usamos a perda quadrática,

`$$L(y_i, \hat y_i) = \frac{1}{2}(y_i - \hat y_i)^2,$$`

então,

`$$-g_i = (y_i - \hat y_i^{(b-1)}) = r_{ib}, \; \text{e} \; h_i = 1.$$`

- Substituindo na fórmula geral, obtemos:

`$$\gamma_{jb} = \frac{\sum_{i: \mathbf{x}_i \in R_{jb}} r_{ib}}{|R_{jb}| + \color{blue}{\lambda}}.$$`

`$$SS_{reg}(\tau_{jb}) = \frac{\sum_{i: \mathbf{x}_i \in R_{jb}} r_{ib}}{|R_{jb}| + \color{blue}{\lambda}}.$$`

---

``` r
library(dplyr)
library(caret)
library(xgboost)
```

``` r
# Carregando os dados
mlb <- read.csv("mlb.txt")
# Removendo os nomes dos jogadores
mlb <- mlb[,-18]
# Renomeando as vars indicadoras
names(mlb)[c(14,15,16,17)] <- c("FAE","FA","AE","A") 
colnames(mlb)
```

```
##  [1] "Salary"  "AVG"     "OBP"     "Runs"    "Hits"    "Doubles" "Triples"
##  [8] "HR"      "RBI"     "Walks"   "SO"      "SB"      "Errs"    "FAE"    
## [15] "FA"      "AE"      "A"
```

---

``` r
# separar preditoras e resposta
predictors <- data.matrix(mlb[,-which(names(mlb) %in% "Salary")])
target <- as.numeric(mlb$Salary)
```

---

``` r
# xgboost com validação cruzada
set.seed(12345)
xgb <- xgb.cv(data = predictors,
              label = target,
              nfold = 10,
              objective = "reg:squarederror",
              nrounds = 250,
              verbose = 0,
              params = list(
               eta = 0.1, 
               max_depth = 3,
               subsample = 0.8)
              )
xgb$evaluation_log[250,]
```

```
##     iter train_rmse_mean train_rmse_std test_rmse_mean test_rmse_std
##    <num>           <num>          <num>          <num>         <num>
## 1:   250        106.8174        4.55324       700.3487      157.8793
```

---

``` r
# Busca por grid
hyper_grid <- expand.grid(
  eta = 0.01,
  max_depth = 3,
  subsample = 0.8,
  min_child_weight = 3,
  gamma = c(0,1,10),
  lambda = c(0,0.1,1,100,1000),
  nrounds = 0, # objeto p/ salvar o número de árvores
  rmse = 0   # objeto p/ salvar os rmse
)
```

---

``` r
# treinando os modelos
for(i in seq_len(nrow(hyper_grid))) {
  set.seed(12345)
  m <- xgb.cv(
    data = predictors,
    label = target,
    nfold = 10,
    objective = "reg:squarederror",
    early_stopping_rounds = 50,
    nrounds = 1000,
    verbose = 0,
    params = list( 
      eta = hyper_grid$eta[i], 
      max_depth = hyper_grid$max_depth[i],
      min_child_weight = hyper_grid$min_child_weight[i],
      subsample = hyper_grid$subsample[i],
      gamma = hyper_grid$gamma[i], 
      lambda = hyper_grid$lambda[i]) 
  )
  hyper_grid$rmse[i] <- min(m$evaluation_log$test_rmse_mean)
  hyper_grid$nrounds[i] <- m$best_iteration
}
```

---

``` r
# resultados (modelo ganhador)
hyper_grid[which.min(hyper_grid$rmse),]
```

```
##    eta max_depth subsample min_child_weight gamma lambda nrounds     rmse
## 4 0.01         3       0.8                3     0    0.1     676 677.2596
```

---

## Parâmetros ótimos

``` r
# lista de hiperparâmetros ótimos
params <- list(
  eta = 0.01,
  max_depth = 3,
  min_child_weight = 3,
  subsample = 0.8,
  lambda = 0.1
)
```

---

## Treinamento do mehor modelo

``` r
xgb_melhor <- xgb.cv(
    data = predictors,
    label = target,
    nfold = 10,
    objective = "reg:squarederror",
    nrounds = 67,
    verbose = 0,
    params = params
)
xgb_melhor$evaluation_log[67,]   
```

```
##     iter train_rmse_mean train_rmse_std test_rmse_mean test_rmse_std
##    <num>           <num>          <num>          <num>         <num>
## 1:    67        1047.232       13.94487       1078.453       168.516
```