Regularização

Bruna Wundervald

Regularização

\[min_{\beta} \sum_{i = 1}^{n}\Big(y_i - \beta_0 - \sum_{j = 1}^{p} \beta_j x_{ij} \Big)^{2} + \lambda P(\beta)\]

Onde \(P(\beta)\) faz o papel de penalidade, mantendo as estimativas de \(\beta_j\) perto de 0. Utilizando a família das potências, temos:

\[min_{\beta} \sum_{i = 1}^{n}\Big(y_i - \beta_0 - \sum_{j = 1}^{p} \beta_j x_{ij} \Big)^{2} + \lambda \sum_{j = 1}^{p} |\beta_j|^q\]

Ou seja, quanto maior é o valor absoluto do coeficiente, mais penalidade é atribuída à ele. Consideramos \(\lambda\) como o parâmetro de tuning.

Penalização Ridge

\[min_{\beta} \sum_{i = 1}^{n}\Big(y_i - \beta_0 - \sum_{j = 1}^{p} \beta_j x_{ij} \Big)^{2} + \lambda \sum_{j = 1}^{p} \beta_j^2\]

\[ \hat \beta = (\mathbf X^T \mathbf X + \lambda \mathbf I)^{-1} \mathbf X^T \mathbf y\]

Notar que \(\hat \beta_{\lambda}^{R} = \frac{\hat \beta^{OLS}}{ 1 + \lambda}\), caracterizando algo importante na regressão ridge: o “encolhimento” dos parâmetros.

Penalização Lasso

\[min_{\beta} \sum_{i = 1}^{n}\Big(y_i - \beta_0 - \sum_{j = 1}^{p} \beta_j x_{ij} \Big)^{2} + \lambda \sum_{j = 1}^{p} |\beta_j|\]

Elastic Net

\[min_{\beta} \sum_{i = 1}^{n}\Big(y_i - \beta_0 - \sum_{j = 1}^{p} \beta_j x_{ij} \Big)^{2} + \lambda \sum_{j = 1}^{p} \Big(\alpha |\beta_j| + (1 - \alpha) \beta^{2}_j \Big)\]

Selecionando \(\lambda\) ótimo

Exemplos

set.seed(20172)
# Carregamento de pacotes
library(tidyverse)
library(plyr)
library(dplyr)
library(glmnet)
library(glmnetUtils)

# Leitura e organização da base:
# 1. Leitura
# 2. Filtra apenas quem mora no centro e tem renda maior que 0
# 3. Seleciona apenas as colunas de interesse
# 4. Transforma algumas colunas em fator e decide quem fará parte da 
# amostra de treino e da de teste

db <- read.csv("dados2.txt", 
               header = TRUE, 
               sep = "\t",
               encoding = "UTF-8") %>% 
  filter(nom_localidade_fam == "CENTRO",
         vlr_renda_total_fam > 0) %>% 
  dplyr::select(c(vlr_renda_total_fam,
                  qtd_comodos_domic_fam,    qtd_comodos_dormitorio_fam,
                  cod_material_piso_fam,    cod_material_domic_fam,
                  cod_agua_canalizada_fam,  cod_escoa_sanitario_domic_fam,
                  cod_iluminacao_domic_fam,
                  cod_calcamento_domic_fam, 
                  qtd_pessoas_domic_fam,    
                  val_desp_energia_fam, val_desp_agua_esgoto_fam,   
                  val_desp_gas_fam, val_desp_alimentacao_fam,
                  val_desp_transpor_fam,    val_desp_aluguel_fam)) %>% 
  dplyr::mutate(cod_material_domic_fam = factor(cod_material_domic_fam),
                cod_agua_canalizada_fam = factor(cod_agua_canalizada_fam),
                cod_calcamento_domic_fam = factor(cod_calcamento_domic_fam),
                cod_escoa_sanitario_domic_fam = factor(cod_escoa_sanitario_domic_fam),
                part = ifelse(runif(695) > 0.3, "treino", "teste"))
dim(db)
[1] 695  17
#-------------------------------------------------------------
# Regressão Ridge
#-------------------------------------------------------------

# Amostra de treino
db.t <- db %>% 
  filter(part == "treino") %>%
  dplyr::select(-c(part))

db.te <- db %>% 
  filter(part == "teste") %>% 
  dplyr::select(-c(part))

# Modelo linear simples
m0 <- lm(vlr_renda_total_fam ~ ., data = db.t)
summary(m0)        

Call:
lm(formula = vlr_renda_total_fam ~ ., data = db.t)

Residuals:
    Min      1Q  Median      3Q     Max 
-1562.4  -445.1   -62.6   300.9  4945.4 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    -228.51143  222.65542  -1.026 0.305320    
qtd_comodos_domic_fam           184.54097   26.54753   6.951 1.33e-11 ***
qtd_comodos_dormitorio_fam     -199.46653   56.49227  -3.531 0.000458 ***
cod_material_piso_fam            64.44493   32.38802   1.990 0.047241 *  
cod_material_domic_fam2         753.44908  232.85016   3.236 0.001306 ** 
cod_material_domic_fam3          76.72224   85.31339   0.899 0.368991    
cod_material_domic_fam6          26.84502  203.85007   0.132 0.895290    
cod_agua_canalizada_fam2         31.56261  673.89693   0.047 0.962665    
cod_escoa_sanitario_domic_fam2 -225.04918  401.64024  -0.560 0.575546    
cod_escoa_sanitario_domic_fam6 -174.23231  514.77612  -0.338 0.735178    
cod_iluminacao_domic_fam       -195.68158   84.27931  -2.322 0.020704 *  
cod_calcamento_domic_fam2       278.71042  166.66672   1.672 0.095192 .  
cod_calcamento_domic_fam3        63.56178  257.15809   0.247 0.804893    
qtd_pessoas_domic_fam            -0.07274   28.57232  -0.003 0.997970    
val_desp_energia_fam              1.36104    0.65327   2.083 0.037796 *  
val_desp_agua_esgoto_fam         -0.52575    0.83514  -0.630 0.529328    
val_desp_gas_fam                 -0.99182    2.08865  -0.475 0.635125    
val_desp_alimentacao_fam          1.85294    0.21553   8.597  < 2e-16 ***
val_desp_transpor_fam             1.41330    0.66910   2.112 0.035236 *  
val_desp_aluguel_fam              0.02905    0.11402   0.255 0.799007    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 666.7 on 435 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.3373,    Adjusted R-squared:  0.3084 
F-statistic: 11.65 on 19 and 435 DF,  p-value: < 2.2e-16
# Aqui, existem evidências a favor da não nulidade dos coecientes
# para 'qtd_comodos_domic_fam', 'qtd_comodos_dormitorio_fam',
# 'cod_material_domic_fam', 'val_desp_energia_fam' e 
# 'val_desp_alimentacao_fam'

# Ridge
m0.r <- glmnetUtils::cv.glmnet(vlr_renda_total_fam ~ .,
                               data = db.t,
                                   alpha = 0)

plot(m0.r)

# O melhor valor de lambda:
m0.r$lambda.1se
[1] 3179.415
# Coeficientes estimados usando esse lambda
coef(m0.r, s = "lambda.1se")
24 x 1 sparse Matrix of class "dgCMatrix"
                                          1
(Intercept)                    662.28899975
qtd_comodos_domic_fam           32.72997474
qtd_comodos_dormitorio_fam      10.55583452
cod_material_piso_fam           24.15722402
cod_material_domic_fam1          6.33052209
cod_material_domic_fam2        102.96587096
cod_material_domic_fam3        -12.34271292
cod_material_domic_fam6        -47.31151920
cod_agua_canalizada_fam1        19.83967374
cod_agua_canalizada_fam2       -19.83975444
cod_escoa_sanitario_domic_fam1 -15.35781671
cod_escoa_sanitario_domic_fam2  19.69299406
cod_escoa_sanitario_domic_fam6   8.66265356
cod_iluminacao_domic_fam       -50.84883561
cod_calcamento_domic_fam1      -20.00855536
cod_calcamento_domic_fam2       15.92708883
cod_calcamento_domic_fam3       29.09089523
qtd_pessoas_domic_fam           10.80824732
val_desp_energia_fam             0.43364860
val_desp_agua_esgoto_fam         0.30674864
val_desp_gas_fam                 0.70682315
val_desp_alimentacao_fam         0.40572351
val_desp_transpor_fam            0.37462491
val_desp_aluguel_fam             0.02170858
# Predição
pred <- predict(m0.r, db.te, type = "response")
real <- db %>% filter(part == "teste") %>% with(vlr_renda_total_fam)

ind <- intersect(rownames(pred), 1:length(real))

# Predito x Observado
data.frame(pred = pred, r = real[ind]) %>%
  ggplot(aes(x = r, y = pred)) + 
  geom_point( colour = "turquoise") + 
  xlab("Dados reais") +
  ylab("Preditos")

# Soma de quadrados
sst <- sum((real - mean(real))^2)
sse <- sum((pred - real[ind])^2)

# R quadrado
rsq <- 1 - sse / sst
rsq
[1] 0.1176799
#-------------------------------------------------------------
# Regressão Lasso
#-------------------------------------------------------------
m0.l <- glmnetUtils::cv.glmnet(vlr_renda_total_fam ~ .,
                               data = db.t,
                                   alpha = 1)
plot(m0.l)

# O melhor valor de lambda:
m0.l$lambda.1se
[1] 119.7029
# O melhor valor de lambda está indicado como 204.5921
# Coeficientes estimados usando esse lambda
coef(m0.l, s = "lambda.1se")
24 x 1 sparse Matrix of class "dgCMatrix"
                                       1
(Intercept)                    290.91875
qtd_comodos_domic_fam           76.10593
qtd_comodos_dormitorio_fam       .      
cod_material_piso_fam            .      
cod_material_domic_fam1          .      
cod_material_domic_fam2          .      
cod_material_domic_fam3          .      
cod_material_domic_fam6          .      
cod_agua_canalizada_fam1         .      
cod_agua_canalizada_fam2         .      
cod_escoa_sanitario_domic_fam1   .      
cod_escoa_sanitario_domic_fam2   .      
cod_escoa_sanitario_domic_fam6   .      
cod_iluminacao_domic_fam         .      
cod_calcamento_domic_fam1        .      
cod_calcamento_domic_fam2        .      
cod_calcamento_domic_fam3        .      
qtd_pessoas_domic_fam            .      
val_desp_energia_fam             .      
val_desp_agua_esgoto_fam         .      
val_desp_gas_fam                 .      
val_desp_alimentacao_fam         1.37155
val_desp_transpor_fam            .      
val_desp_aluguel_fam             .      
# Predição
pred.l <- predict(m0.l, db.te, type = "response")
ind <- intersect(rownames(pred.l), 1:length(real))

# Predito x Observado
data.frame(pred = pred.l, r = real[ind]) %>%
  ggplot(aes(x = r, y = pred)) + 
  geom_point( colour = "tomato") + 
  xlab("Dados reais") +
  ylab("Preditos")