Regularização
Quando utilizamos métodos de regressão, podemos nos deparar com bancos de dados que possuem muitas covariáveis em relação à quantidade de observações. Uma suposição que pode ser feita é a de que algumas dessas variáveis tem um alto poder de predição. Assim, um bom método para trabalhar com esse problema é regular a entrada de vairiáveis em um modelo, eliminando-as ou impondo pesos/penalidade nas entradas.
A ideia de utilizar penalidades é aplicada através e Multiplicadores de Lagrange:
\[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
- Neste caso, o objetivo é otimizar:
\[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\]
- A ideia é somar uma constante \(\lambda\) à diagonal da matriz se \(p > n\), solucionando o problema da singularidade:
\[ \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
- Agora, o objetivo é otimizar:
\[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
- A Elastic Net é basicamente a combinação entre os métodos Lasso e 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} \Big(\alpha |\beta_j| + (1 - \alpha) \beta^{2}_j \Big)\]
Selecionando \(\lambda\) ótimo
- Assim como em outros métodos de Machine Learning, utiliza-se validação cruzada para encontrar um parâmetro ótimo, que neste caso é o \(\lambda\). Os passos para isso são:
- Avaliar a taxa de erro de classificação para um certo grid de \(\lambda\);
- Verificar qual dos \(\lambda\) testados produz a menor taxa de erro;
- Ajustar o novo modelo com todas as observações e o \(\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")
sse <- sum((pred.l - real[ind])^2)
# R quadrado
rsq <- 1 - sse / sst
rsq[1] 0.1896636#-------------------------------------------------------------
# Gráfico de perfil dos lambdas
#-------------------------------------------------------------
lambdas <- seq(0, 300, l = 50)
sse <- vector(mode = "numeric", length = length(lambdas))
rsq <- vector(mode = "numeric", length = length(lambdas))
for(i in 1:length(lambdas)){
m0 <- glmnetUtils::glmnet(vlr_renda_total_fam ~ .,
data = db.t,
alpha = 0,
lambda = lambdas[i])
pred <- predict(m0, db.te, type = "response")
ind <- intersect(rownames(pred), 1:length(real))
sse[i] <- sum((pred - real[ind])^2)
# R quadrado
rsq[i] <- 1 - sse[i] / sst
}
# R-quadrado x Lambdas
data.frame(lambdas = lambdas, rsq = rsq) %>%
ggplot(aes(x = lambdas, y = rsq)) +
geom_line( colour = "green2") +
geom_vline(aes(xintercept = lambdas[which.max(rsq)]),
linetype = "dashed", size = 2, color = "gray") +
annotate("text", x = lambdas[which.max(rsq)] + 10, y = max(rsq),
label = paste("lambda == ", round(lambdas[which.max(rsq)], 1)),
parse = TRUE, size = 6)
# Lambdas x sse
data.frame(lambdas = lambdas, sse = sse) %>%
ggplot(aes(x = lambdas, y = sse)) +
geom_line( colour = "yellow") +
geom_vline(aes(xintercept = lambdas[which.min(sse)]),
linetype = "dashed", size = 2, color = "gray") +
annotate("text", x = lambdas[which.min(sse)] + 10, y = min(sse),
label = paste("lambda == ", round(lambdas[which.min(sse)], 1)),
parse = TRUE, size = 6)