|
Machine Learning
|
Regularização em regressão com implementação
Prof. Walmes M. Zeviani
2018-09-10
Esse tutorial irá explorar computacionalmente conceitos relacionados à regularização. Os exemplos irão explorar a função objetivo acrescida da penalização, será feita estimação dos parâmetros sobre penalização \(L_2\) e \(L_1\) e será empregada validação cruzada para determinação do parâmetro de penalidade.
1 Inspeção da função objetivo com a penalização
Serão utilizados dados simulados e reais para materializar conceitos de regularização. No caso, a resposta será simulada como função de duas variáveis apenas, uma condição na qual raramente se aplica regularização, porém, será útil para produzir gráficos e estudar a função objetivo.
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(lattice)
library(latticeExtra)
library(mvtnorm) # Para simular da normal multivarida.
library(MASS) # lm.ridge()
library(penalized) # LASSO.
library(glmnet) # L1, L2 e elastic net.
library(caret)
sessionInfo()## R version 3.4.4 (2018-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods stats graphics grDevices datasets utils base
##
## other attached packages:
## [1] caret_6.0-77 ggplot2_2.2.1 glmnet_2.0-12
## [4] foreach_1.4.3 Matrix_1.2-14 penalized_0.9-50
## [7] survival_2.41-3 MASS_7.3-50 mvtnorm_1.0-8
## [10] latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-35
## [13] rmarkdown_1.8.10 knitr_1.20.12
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.17 lubridate_1.7.3 class_7.3-14
## [4] assertthat_0.2.0 rprojroot_1.3-2 digest_0.6.14
## [7] ipred_0.9-6 R6_2.2.2 plyr_1.8.4
## [10] backports_1.1.2 stats4_3.4.4 evaluate_0.11
## [13] pillar_1.1.0 rlang_0.2.1 lazyeval_0.2.0
## [16] kernlab_0.9-25 rpart_4.1-13 splines_3.4.4
## [19] CVST_0.2-1 ddalpha_1.3.1 gower_0.1.2
## [22] stringr_1.3.1 munsell_0.4.3 compiler_3.4.4
## [25] xfun_0.3 pkgconfig_2.0.1 dimRed_0.1.0
## [28] htmltools_0.3.6 nnet_7.3-12 tidyselect_0.2.4
## [31] tibble_1.4.2 prodlim_1.6.1 DRR_0.0.2
## [34] codetools_0.2-15 RcppRoll_0.2.2 withr_2.0.0
## [37] dplyr_0.7.5 recipes_0.1.0 ModelMetrics_1.1.0
## [40] grid_3.4.4 nlme_3.1-131 gtable_0.2.0
## [43] magrittr_1.5 scales_0.5.0 stringi_1.2.4
## [46] reshape2_1.4.2 bindrcpp_0.2.2 timeDate_3012.100
## [49] robustbase_0.92-7 lava_1.5.1 iterators_1.0.8
## [52] tools_3.4.4 glue_1.3.0 DEoptimR_1.0-8
## [55] purrr_0.2.4 sfsmisc_1.1-1 yaml_2.2.0
## [58] colorspace_1.3-2 bindr_0.1.1Os dados serão simulados de uma distribição \(3\)-normal. O vetor de médias é 0 e a matriz de covariância estabelece uma alta correlação entre as variáveis indenpendentes x1 e x2.
#-----------------------------------------------------------------------
# Simulando um conjunto de dados.
set.seed(123)
# Criando a matriz de covariância do vetor de v.a. [y, x_1, x_2].
sigma <- matrix(0.98, 3, 3) + 0.02 * diag(3)
sigma[1, 2:3] <- 0.5
sigma[2:3, 1] <- 0.5
sigma## [,1] [,2] [,3]
## [1,] 1.0 0.50 0.50
## [2,] 0.5 1.00 0.98
## [3,] 0.5 0.98 1.00da <- rmvnorm(n = 100, mean = c(10, 5, 5), sigma = sigma)
da <- as.data.frame(da)
names(da) <- c("y", "x1", "x2")
head(da)## y x1 x2
## 1 9.751545 5.665699 5.918686
## 2 10.463534 6.170010 6.394272
## 3 10.019637 3.717290 3.799061
## 4 9.915806 6.053404 5.931178
## 5 10.286147 4.827454 4.733193
## 6 11.386741 4.549304 4.200775# Pares de diagramas de dispersão.
splom(da, type = c("p", "smooth"), as.matrix = TRUE)
# Ajuste do modelo sem penalização.
m0 <- lm(y ~ x1 + x2, data = da)
summary(m0)##
## Call:
## lm(formula = y ~ x1 + x2, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.62112 -0.50295 0.05175 0.41140 1.82286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.8696 0.4197 18.750 <2e-16 ***
## x1 0.6298 0.3800 1.657 0.101
## x2 -0.1881 0.3759 -0.500 0.618
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7615 on 97 degrees of freedom
## Multiple R-squared: 0.234, Adjusted R-squared: 0.2182
## F-statistic: 14.81 on 2 and 97 DF, p-value: 2.434e-06Aspectos interessantes do exercício:
- O modelo de regressão é globalmente significativo.
- Não há efeito de nenhuma das variáveis predidoras
x1ex2. - Embora não significativo, o sinal o efeito da variável
x2é contrário ao que se observa para a relação deycomx2no diagrama de dispersão.
A não significância das preditoras é resultado da alta correlação entre elas, o fato de serem praticamente reduntantes. O quadro de estimativas está testando a hipótese nula do efeito de uma variável ser zero na presença das demais. Conclui-se que x2 não tem efeito em y na presença de x1 e vice-versa.
Do ponto de vista de otimização, a função objetivo é “mal condicionada” conforme pode ser visto nos gráficos a seguir. A função objetivo foi escrita considerando o parâmetro \(\beta_0\) (intercepto) como conhecido para que se possa inspecionar a função nas demais variáveis do problema de otimização, os parâmetros \(\beta_1\) e \(\beta_2\).
# Função para traçar os contornos de nível da função objetivo.
rss <- function(beta1, beta2, beta0 = 7, X, y) {
# crossprod(A) := t(A) %*% A.
crossprod(y - X %*% cbind(c(beta0, beta1, beta2)))
}
RSS <- Vectorize(FUN = rss, vectorize.args = c("beta1", "beta2"))
l2 <- function(beta1, beta2, alpha) {
# crossprod(A) := t(A) %*% A.
alpha * crossprod(c(beta1, beta2))
}
L2 <- Vectorize(FUN = l2, vectorize.args = c("beta1", "beta2"))
l1 <- function(beta1, beta2, alpha) {
# crossprod(A) := t(A) %*% A.
alpha * sum(abs(c(beta1, beta2)))
}
L1 <- Vectorize(FUN = l1, vectorize.args = c("beta1", "beta2"))
# Grid de valores dos parâmetros para a avaliação da função em cada
# ponto.
grid <- expand.grid(beta1 = seq(from = -1, to = 1, length.out = 50),
beta2 = seq(from = -1, to = 1, length.out = 50),
KEEP.OUT.ATTRS = FALSE)
# Avalia a função objetivo em cada ponto.
grid$RSS <- with(grid,
RSS(beta1 = beta1,
beta2 = beta2,
beta0 = coef(m0)[1],
X = model.matrix(m0),
y = m0$model[, 1]))
# Avalia a penalidade L2 em cada ponto.
grid$L2 <- with(grid,
L2(beta1 = beta1,
beta2 = beta2,
alpha = 1000))
# Avalia a penalidade L1 em cada ponto.
grid$L1 <- with(grid,
L1(beta1 = beta1,
beta2 = beta2,
alpha = 1000))
# Função objetivo somada com as penalidades.
grid$RSS_L1 <- with(grid, RSS + L1)
grid$RSS_L2 <- with(grid, RSS + L2)
str(grid)## 'data.frame': 2500 obs. of 7 variables:
## $ beta1 : num -1 -0.959 -0.918 -0.878 -0.837 ...
## $ beta2 : num -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ RSS : num 15640 15124 14616 14118 13627 ...
## $ L2 : num 2000 1920 1843 1770 1700 ...
## $ L1 : num 2000 1959 1918 1878 1837 ...
## $ RSS_L1: num 17640 17083 16535 15995 15464 ...
## $ RSS_L2: num 17640 17044 16460 15888 15327 ...# Valor da função para traçar os contornos.
at <- c(0, 100, 200, 500, 1000, 2000, 5000, 10000, Inf)
# Gráfico de contornos.
contourplot(RSS ~ beta1 + beta2,
data = grid,
at = at,
aspect = "iso",
xlab = expression(beta[1]),
ylab = expression(beta[2]),
key = list(text = list(c("RSS", "L2", "L1")),
lines = list(lty = c(1, 2, 2),
col = c("black", "red", "blue"))),
col.regions = cm.colors(length(at)),
region = TRUE) +
latticeExtra::layer({
# Linhas passando pela origem (0, 0).
panel.abline(v = 0, h = 0, lty = 3)
}) +
latticeExtra::as.layer({
# Contornos da penalização L2 (Ridge).
contourplot(L2 ~ beta1 + beta2,
at = at,
col = "red",
data = grid, lty = 2)
}) +
latticeExtra::as.layer({
# Contornos da penalização L1 (Lasso)
contourplot(L1 ~ beta1 + beta2,
at = at,
col = "blue",
data = grid, lty = 2)
})
p2 <-
contourplot(RSS_L2 ~ beta1 + beta2,
data = grid,
xlab = expression(beta[1]),
ylab = expression(beta[2]),
at = at,
aspect = "iso",
region = TRUE) +
latticeExtra::layer({
panel.abline(v = 0, h = 0, lty = 3)
})
p1 <-
contourplot(RSS_L1 ~ beta1 + beta2,
data = grid,
xlab = expression(beta[1]),
ylab = expression(beta[2]),
at = at,
aspect = "iso",
region = TRUE) +
latticeExtra::layer({
panel.abline(v = 0, h = 0, lty = 3)
})
c("Lasso (L1)" = p1, "Ridge (L2)" = p2)## Warning in formals(fun): argument is not a function
#-----------------------------------------------------------------------
# Variando o parâmetro de penalização para observar o efeito da função
# objetivo.
grid <- grid[, 1:7]
alpha_seq <- c(0, 50, 100, 500)
for (alpha in alpha_seq) {
grid$x <- with(grid, RSS + L2(beta1 = beta1,
beta2 = beta2,
alpha = alpha))
names(grid)[ncol(grid)] <- sprintf("RSS_L2_%d", alpha)
grid$x <- with(grid, RSS + L1(beta1 = beta1,
beta2 = beta2,
alpha = alpha))
names(grid)[ncol(grid)] <- sprintf("RSS_L1_%d", alpha)
}
# Empilha.
gridw <- reshape2::melt(grid,
id.vars = grep("beta\\d", names(grid)),
measure.vars = grep("L\\d_", names(grid)))
str(gridw)## 'data.frame': 20000 obs. of 4 variables:
## $ beta1 : num -1 -0.959 -0.918 -0.878 -0.837 ...
## $ beta2 : num -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ variable: Factor w/ 8 levels "RSS_L2_0","RSS_L1_0",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 15640 15124 14616 14118 13627 ...# Parte em duas variáveis.
aux <- as.data.frame(do.call(rbind,
strsplit(as.character(gridw$variable),
split = "_")))[, -1]
names(aux) <- c("norma", "lambda")
aux$lambda <- as.numeric(as.character(aux$lambda))
# Concatena colunas.
gridw <- cbind(gridw, aux)
# Gráfico de contornos.
useOuterStrips(
contourplot(value ~ beta1 + beta2 | factor(lambda) + norma,
data = gridw,
at = at,
xlab = expression(beta[1]),
ylab = expression(beta[2]),
aspect = "iso",
colorkey = FALSE,
as.table = TRUE,
region = TRUE)) +
latticeExtra::layer({
panel.abline(v = 0, h = 0, lty = 3)
})
- IMPORTATE: Continua sendo uma função convexa para ambas as penalidades. Para a penalidade \(L_2\), continua sendo uma função elíptica.
- Os contornos ficam menos alongados, indo para mais circulares, portanto, um problema de otimização melhor condicionado.
- DETALHE: A solução é mesma porque não foi feita estimação dos parâmetros sujeitos à penalidade. Isso será feito nas sessões seguintes.
2 Carregando dados da web para as aplicações
Os dados a seguir são de medidas de comprimento em folhas de uva com os quais deseja-se obter um modelo preditivo para a área da folha. Veja a descrição na sessão Carregando os dados do tutorial de regressão apresentado no MGEst.
#-----------------------------------------------------------------------
# Dados hospedados na web.
url <- "http://www.leg.ufpr.br/~walmes/data/areafoliarUva.txt"
uva <- read.table(url,
header = TRUE,
sep = "\t",
stringsAsFactors = FALSE)
uva$cult <- factor(uva$cult)
str(uva)## 'data.frame': 300 obs. of 9 variables:
## $ id : chr "malbec_1.jpg" "malbec_10.jpg" "malbec_11.jpg" "malbec_12.jpg" ...
## $ cult: Factor w/ 3 levels "malbec","merlot",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ area: num 100.8 85.8 119.5 137 84.7 ...
## $ mc : num 12 11.5 12.5 15.5 10 12 15.5 17.5 13.5 13.3 ...
## $ nc : num 7.5 9 8.5 10 7 8.5 11 13 10 9.5 ...
## $ ml : num 12.8 10.5 13 14.4 11 12 14 14 12 15 ...
## $ nld : num 6.4 8.5 8.6 9 6.5 9 10 9.5 9 9.2 ...
## $ nle : num 7.5 7 9 10 7 8.2 11 11 8.5 8.3 ...
## $ cll : num 9.5 9.5 10.2 12 7.5 8.9 13.5 10.8 9.7 10.3 ...# Comprimento da nervura lateral: média dos lados direito e esquerdo.
uva$nl <- with(uva, apply(cbind(nld, nle), 1, mean))
uva <- subset(uva, select = -c(nld, nle))
#-----------------------------------------------------------------------
# Ver.
splom(uva[-(1:2)],
groups = uva$cult,
auto.key = TRUE,
cex = 0.2,
type = c("p", "smooth"),
as.matrix = TRUE)
# Existe alta correlação entre as medidas de comprimento.
round(cor(uva[-(1:2)]), digits = 2)## area mc nc ml cll nl
## area 1.00 0.95 0.92 0.93 0.86 0.95
## mc 0.95 1.00 0.96 0.88 0.89 0.94
## nc 0.92 0.96 1.00 0.84 0.89 0.92
## ml 0.93 0.88 0.84 1.00 0.78 0.92
## cll 0.86 0.89 0.89 0.78 1.00 0.87
## nl 0.95 0.94 0.92 0.92 0.87 1.00#-----------------------------------------------------------------------
# Análise apenas para a variedade Malbec.
mal <- subset(uva,
cult == "malbec",
select = -c(cult, id))
# A análise com o log_{10} tem os pressupostos mais atendidos.
mal$larea <- log10(mal$area)
dim(mal)## [1] 100 7# IMPORTANT: Normaliza todas as variáveis: média 0 e variância 1.
sc_vars <- scale(mal)
# Nos atributos estão as médias e desvios usados na normalização.
tail(attributes(sc_vars), n = 2)## $`scaled:center`
## area mc nc ml cll nl
## 116.447477 12.941000 9.410000 12.263000 10.483000 8.354000
## larea
## 2.011195
##
## $`scaled:scale`
## area mc nc ml cll nl
## 46.8878181 3.4047826 2.6483271 2.8235895 2.6922656 1.9981896
## larea
## 0.2617036# Versão original e normalizada.
mal0 <- mal
mal <- as.data.frame(sc_vars)
str(mal)## 'data.frame': 100 obs. of 7 variables:
## $ area : num -0.3347 -0.6528 0.0653 0.4382 -0.678 ...
## $ mc : num -0.276 -0.423 -0.13 0.752 -0.864 ...
## $ nc : num -0.721 -0.155 -0.344 0.223 -0.91 ...
## $ ml : num 0.19 -0.624 0.261 0.757 -0.447 ...
## $ cll : num -0.365 -0.365 -0.105 0.563 -1.108 ...
## $ nl : num -0.703 -0.302 0.223 0.574 -0.803 ...
## $ larea: num -0.0303 -0.2962 0.253 0.4796 -0.3192 ...# Apenas efeitos aditivos dos comprimentos.
m0 <- lm(larea ~ mc + ml + nc + nl + cll, data = mal)
summary(m0)##
## Call:
## lm(formula = larea ~ mc + ml + nc + nl + cll, data = mal)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4436 -0.1128 0.0484 0.1268 0.6935
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.116e-17 2.531e-02 0.000 1.000000
## mc 2.420e-01 1.075e-01 2.252 0.026656 *
## ml 5.842e-01 5.975e-02 9.776 5.5e-16 ***
## nc -1.239e-02 8.302e-02 -0.149 0.881652
## nl 3.524e-01 1.013e-01 3.479 0.000764 ***
## cll -1.767e-01 7.233e-02 -2.443 0.016443 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2531 on 94 degrees of freedom
## Multiple R-squared: 0.9392, Adjusted R-squared: 0.9359
## F-statistic: 290.3 on 5 and 94 DF, p-value: < 2.2e-163 Regularização com norma \(L_2\)
Nos códigos a seguir será feito regularização com norma \(L_2\), conhecida por regressão Ridge.
3.1 Estudo da função objetivo com penalização \(L_2\)
# Para estimar os parâmetros considerando a penalidade de norma L2.
beta_ridge <- function(Xly, XlX, alpha) {
p <- ncol(XlX) - 1
d <- rep(c(0, 1), c(1, p))
beta <- solve(XlX + alpha * diag(d)) %*% Xly
return(beta)
}
# Extrai componentes do modelo declarado acima.
X <- cbind(1, da$x1, da$x2) # Matriz do modelo.
y <- cbind(da$y) # Vetor da resposta.
XlX <- crossprod(X) # X'X.
Xly <- crossprod(X, y) # X'y
# Grid de valores para os parâmetros.
grid <- expand.grid(beta1 = seq(-1, 1, length.out = 50),
beta2 = seq(-1, 1, length.out = 50))
# Grid de valores para o parâmetro de penalidade.
alpha_grid <- c(0, 20, 50, 100, 200, 500, 1000, 2000, 5000)
# Lista vazia para ser preenchida pelo laço.
betas <- vector(mode = "list", length = length(alpha_grid))
i <- 1
for (alpha in alpha_grid) {
# Determina os parâmetros sujeito à penalidade.
beta <- beta_ridge(Xly, XlX, alpha)
# Avalia a função objetivo em cada ponto da malha.
grid$x <- with(grid,
RSS(beta0 = beta[1],
beta1 = beta1,
beta2 = beta2,
X = X,
y = y) +
L2(beta1 = beta1,
beta2 = beta2,
alpha = alpha))
names(grid)[ncol(grid)] <- sprintf("RSS_L2_%d", alpha)
betas[[i]] <- beta
i <- i + 1
}
# Empilha.
gridw <- reshape2::melt(grid,
id.vars = grep("beta\\d", names(grid)),
measure.vars = grep("L2_", names(grid)))
gridw$lambda <- as.integer(gsub(".*_(\\d+)$", "\\1", gridw$variable))
str(gridw)## 'data.frame': 22500 obs. of 5 variables:
## $ beta1 : num -1 -0.959 -0.918 -0.878 -0.837 ...
## $ beta2 : num -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ variable: Factor w/ 9 levels "RSS_L2_0","RSS_L2_20",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 15640 15124 14616 14118 13627 ...
## $ lambda : int 0 0 0 0 0 0 0 0 0 0 ...# Gráfico de contornos.
contourplot(value ~ beta1 + beta2 | factor(lambda),
data = gridw,
cuts = 20,
xlab = expression(beta[1]),
ylab = expression(beta[2]),
layout = c(3, 3),
aspect = "iso",
as.table = TRUE,
region = TRUE) +
latticeExtra::layer({
panel.abline(v = 0, h = 0, lty = 3)
panel.points(x = betas[[which.packet()]][2],
y = betas[[which.packet()]][3],
pch = 19)
})
3.2 Usando o pacote MASS
O pacote MASS tem a função lm.ridge() que ajusta modelo com penalização \(L_2\) para um valor fixado (ou vetor de valores) do parâmetro de penalidade.
X <- model.matrix(m0) # Matriz do modelo.
y <- cbind(m0$model[, 1]) # Vetor da resposta.
XlX <- crossprod(X) # X'X.
Xly <- crossprod(X, y) # X'y
round(beta_ridge(XlX = XlX, Xly = Xly, alpha = 10), digits = 4)## [,1]
## (Intercept) 0.0000
## mc 0.1903
## ml 0.4474
## nc 0.0621
## nl 0.2620
## cll 0.0049library(MASS)
m0r <- lm.ridge(larea ~ mc + ml + nc + nl + cll, data = mal,
lambda = 10)
round(cbind(coef(m0r)), digits = 4)## [,1]
## 0.0000
## mc 0.1905
## ml 0.4483
## nc 0.0617
## nl 0.2624
## cll 0.0040# Avaliando em um grid de valores para o parâmetro de penalidade.
m0r <- lm.ridge(larea ~ mc + ml + nc + nl + cll, data = mal,
lambda = seq(0, 10, 0.01))
# Estimativa baseada em GCV.
select(m0r)## modified HKB estimator is 0.3495133
## modified L-W estimator is 0.2066886
## smallest value of GCV at 1.19# Mostra o traço para os betas associados às variáveis.
plot(m0r)
abline(h = 0, lty = 2)
abline(v = 1.19, lty = 2, col = 2)
# Mostra as estimativas com cada valor de lambda.
round(head(coef(m0r), n = 10), digits = 4)## mc ml nc nl cll
## 0.00 0 0.2420 0.5842 -0.0124 0.3524 -0.1767
## 0.01 0 0.2419 0.5840 -0.0122 0.3522 -0.1763
## 0.02 0 0.2417 0.5838 -0.0121 0.3520 -0.1760
## 0.03 0 0.2416 0.5836 -0.0119 0.3518 -0.1756
## 0.04 0 0.2414 0.5834 -0.0118 0.3516 -0.1752
## 0.05 0 0.2413 0.5832 -0.0116 0.3514 -0.1749
## 0.06 0 0.2411 0.5830 -0.0115 0.3512 -0.1745
## 0.07 0 0.2410 0.5828 -0.0113 0.3510 -0.1741
## 0.08 0 0.2408 0.5827 -0.0111 0.3508 -0.1738
## 0.09 0 0.2407 0.5825 -0.0110 0.3506 -0.1734# Estimativas sem e com penalização.
round(cbind("Sem penal." = coef(m0),
"Com penal." = coef(m0r)[" 0.19", ]),
digits = 4)## Sem penal. Com penal.
## (Intercept) 0.0000 0.0000
## mc 0.2420 0.2393
## ml 0.5842 0.5806
## nc -0.0124 -0.0095
## nl 0.3524 0.3486
## cll -0.1767 -0.16983.3 Validação cruzada para o parâmetro de penalização com norma \(L_2\)
A validação cruzada é geralmente usada para definir o valor apropriado para o parâmetro de penalização. Para alguns modelos ou implementaçõespode já haver um método embutido para determinação do parâmetro de penalização. Foi o que vimos com o pacote MASS para regressão linear múltipla. Mesmo assim, nessa seção será realizada a implementação do zero para determinação do parâmetro de penalização por validação cruzada \(k\)-fold.
# Validação cruzada 8-fold implementada do zero.
# Por que 8? Ué, por que não?
# Embaralha os registros para remover qualquer viés sistemático de
# aquisição de dados, e.g. efeito cronológico ou safra.
set.seed(987)
db <- mal[sample(1:nrow(mal)), ]
n <- nrow(db) # Número de registros.
k <- 8 # Número de folds.
db$i <- ceiling((1:n)/(n/k))
# Número de registros por fold.
table(db$i)##
## 1 2 3 4 5 6 7 8
## 12 13 12 13 12 13 12 13# Cada fold.
by(db, db$i, round, digits = 4)## db$i: 1
## area mc nc ml cll nl larea i
## 48 0.8560 0.4579 0.3738 0.7923 1.0092 0.9238 0.7014 1
## 99 0.6738 0.6047 0.2228 0.7923 0.3777 1.0740 0.6083 1
## 60 1.4734 1.3684 1.1668 0.6152 1.3063 1.4493 0.9829 1
## 56 -0.2797 -0.7169 -0.9100 0.3673 -0.3651 -0.2522 0.0116 1
## 78 0.8990 1.1334 1.2801 0.6152 0.7492 0.8237 0.7226 1
## 23 -0.2922 0.3110 0.4116 0.0131 0.2663 -0.3023 0.0022 1
## 40 -0.3455 0.2523 0.1473 -0.5536 0.4149 -0.1271 -0.0387 1
## 33 -0.4309 -0.5995 -1.0610 0.1548 -0.6251 -0.6526 -0.1063 1
## 80 -1.2003 -1.1575 -1.4764 -1.0848 -1.6651 -1.4783 -0.8859 1
## 65 0.6693 0.8984 1.1668 0.6152 0.5263 0.6986 0.6059 1
## 17 0.2597 0.7222 1.0535 0.2610 1.0092 0.5735 0.3750 1
## 34 1.1319 1.2215 1.0913 1.3589 1.0835 0.8488 0.8331 1
## --------------------------------------------------------
## db$i: 2
## area mc nc ml cll nl larea i
## 15 -0.7267 -0.6876 -0.8345 -0.3056 -1.0337 -0.7026 -0.3645 2
## 28 -0.2030 -0.1589 -0.1171 0.3673 0.0806 -0.2522 0.0684 2
## 90 -1.4874 -1.5980 -1.4386 -1.4389 -0.9223 -1.2782 -1.3062 2
## 55 0.1390 -0.2764 -0.3059 0.0485 0.0063 0.2732 0.3003 2
## 12 -0.8108 -1.0106 -1.2121 -0.4827 -1.2194 -0.9278 -0.4460 2
## 81 0.6931 0.6047 0.2228 0.6506 0.8606 0.5485 0.6184 2
## 9 -0.0674 0.1642 0.2228 -0.0931 -0.2908 0.1982 0.1643 2
## 58 1.9608 1.9264 1.5066 -0.3765 0.4520 0.9489 1.1756 2
## 63 -0.2346 -0.5701 -0.3436 0.0839 -0.5137 -0.1772 0.0452 2
## 94 -2.1751 -2.1855 -2.0428 -3.1035 -2.2223 -2.5543 -3.2516 2
## 88 -1.6246 -1.3925 -0.9855 -1.6869 -1.1080 -1.5534 -1.5521 2
## 83 0.4197 0.6047 0.6004 0.3318 0.0806 0.7737 0.4690 2
## 74 0.1589 0.4579 0.3738 0.0839 0.4520 0.6236 0.3129 2
## --------------------------------------------------------
## db$i: 3
## area mc nc ml cll nl larea i
## 1 -0.3347 -0.2764 -0.7212 0.1902 -0.3651 -0.7026 -0.0303 3
## 37 -0.0698 -0.1295 -0.0793 0.0131 0.4892 0.1481 0.1626 3
## 76 1.1714 1.1921 0.7892 0.9693 0.1920 1.0740 0.8511 3
## 24 -0.6294 -0.4232 -0.1926 -0.4827 -0.3651 -0.1772 -0.2751 3
## 13 -0.2620 0.0173 0.2228 -0.2702 0.4892 -0.0520 0.0249 3
## 32 -0.2121 -0.5701 -0.4947 -0.1994 -0.8480 -0.1021 0.0618 3
## 50 -1.4621 -1.4512 -1.5142 -1.4389 -1.8137 -1.4783 -1.2646 3
## 95 2.4995 1.7796 1.7332 1.4297 2.0492 2.0749 1.3655 3
## 43 0.5507 0.7810 0.8647 0.4735 0.7492 0.7987 0.5423 3
## 8 0.6836 1.3390 1.3556 0.6152 0.1177 0.9489 0.6134 3
## 96 -2.3682 -3.0372 -2.7225 -3.6347 -2.8537 -3.1799 -4.8846 3
## 21 -0.4219 -0.5113 -0.3436 0.0131 0.0063 -0.1521 -0.0991 3
## --------------------------------------------------------
## db$i: 4
## area mc nc ml cll nl larea i
## 89 -1.2442 -0.8932 -0.7212 -1.1556 -0.7737 -1.0279 -0.9436 4
## 53 0.5544 0.7516 0.6759 0.4381 1.0092 0.3483 0.5443 4
## 85 0.1195 0.2523 0.4116 0.0839 0.1920 0.1982 0.2879 4
## 57 -0.0065 -0.5113 -0.6835 0.6506 -0.5508 -0.0771 0.2056 4
## 73 2.0041 1.4859 1.7332 1.9610 1.6778 1.3242 1.1917 4
## 84 2.4522 1.9264 1.5444 1.6777 2.0492 1.5744 1.3497 4
## 14 -0.2713 -0.1002 0.2228 -0.2702 0.2292 -0.1271 0.0180 4
## 4 0.4382 0.7516 0.2228 0.7568 0.5635 0.5735 0.4796 4
## 77 0.2839 0.6047 0.7892 0.2610 0.5635 0.5735 0.3895 4
## 41 -0.7216 -0.8344 -0.9100 -0.1994 -0.4766 -0.6025 -0.3598 4
## 97 1.6336 1.3390 1.3556 1.3235 1.6778 1.4493 1.0487 4
## 87 -1.4796 -0.9813 -0.9100 -1.3681 -0.5508 -1.1781 -1.2931 4
## 25 -1.0343 -1.0987 -0.9100 -0.9431 -0.8851 -0.7527 -0.6839 4
## --------------------------------------------------------
## db$i: 5
## area mc nc ml cll nl larea i
## 7 0.7101 0.7516 0.6004 0.6152 1.1206 1.0740 0.6272 5
## 20 0.2555 0.4579 0.7892 0.0839 0.5635 0.5985 0.3724 5
## 35 0.9918 0.8984 0.6004 1.0756 1.0463 0.7987 0.7675 5
## 27 0.2208 0.1642 0.1850 0.2610 0.4520 0.4734 0.3512 5
## 16 0.9061 1.1921 1.1668 0.4735 1.1206 0.9739 0.7261 5
## 18 0.2133 0.6635 0.7514 0.2256 0.5263 0.9238 0.3466 5
## 3 0.0653 -0.1295 -0.3436 0.2610 -0.1051 0.2232 0.2530 5
## 5 -0.6780 -0.8638 -0.9100 -0.4473 -1.1080 -0.8027 -0.3192 5
## 46 1.4178 1.0453 0.7892 1.2527 1.3806 1.2742 0.9594 5
## 26 -1.4570 -1.5980 -1.6652 -0.9785 -1.5537 -1.6285 -1.2563 5
## 49 0.0427 -0.3351 -0.6457 0.5443 0.7120 -0.0520 0.2382 5
## 68 -0.4916 -0.5113 -0.7967 -0.0223 -0.8851 -0.8528 -0.1562 5
## --------------------------------------------------------
## db$i: 6
## area mc nc ml cll nl larea i
## 31 0.9853 0.8984 0.8269 0.1548 1.1578 0.7487 0.7644 6
## 59 1.5483 1.4859 1.5444 0.8631 1.4921 1.0740 1.0140 6
## 100 -0.8515 -0.8638 -0.9100 -0.8015 -0.6623 -0.4274 -0.4868 6
## 45 -0.4552 -0.0414 0.0717 -0.4119 -0.6994 -0.1271 -0.1261 6
## 93 -2.0820 -2.1855 -1.9673 -2.5723 -1.7766 -2.2290 -2.8140 6
## 62 0.0380 0.1642 0.2228 -0.0931 0.0063 0.4234 0.2351 6
## 72 0.6281 0.5166 0.4871 0.7214 0.3035 0.8738 0.5841 6
## 2 -0.6528 -0.4232 -0.1548 -0.6244 -0.3651 -0.3023 -0.2962 6
## 39 0.3710 0.6928 0.9025 0.3318 0.3406 0.6986 0.4409 6
## 54 0.6267 0.7516 0.4493 0.7214 0.3777 0.6986 0.5833 6
## 11 -1.3986 -1.7449 -1.6652 -1.3327 -1.8509 -1.5534 -1.1645 6
## 22 0.4698 0.5166 0.7137 0.8277 0.7492 0.5985 0.4974 6
## 30 0.5329 -0.1002 -0.2303 0.8985 0.4149 -0.0020 0.5325 6
## --------------------------------------------------------
## db$i: 7
## area mc nc ml cll nl larea i
## 71 0.7921 0.8984 0.7514 0.2256 0.1920 0.9739 0.6693 7
## 47 1.2352 0.8984 0.6004 1.6422 1.4178 0.6986 0.8798 7
## 75 -0.6141 -0.5701 -0.8345 -0.4473 -0.6994 -0.4774 -0.2615 7
## 67 -1.3708 -1.3043 -1.5142 -1.1556 -1.6651 -1.7035 -1.1224 7
## 86 -1.1285 -0.8638 -0.8345 -1.1910 -0.8480 -0.9529 -0.7955 7
## 98 0.8437 -0.1295 2.2996 0.9693 0.5635 0.3233 0.6952 7
## 10 0.4056 0.1054 0.0340 0.9693 -0.0680 0.1982 0.4609 7
## 36 0.1759 0.3698 0.4116 0.2964 -0.7366 -0.4524 0.3235 7
## 29 0.5723 0.6928 0.4493 0.5089 1.3063 0.5735 0.5540 7
## 38 -0.5205 -0.4232 -0.5324 -0.0931 -0.3651 -0.3023 -0.1804 7
## 66 0.5984 0.7810 0.6759 0.6152 1.1206 0.4985 0.5681 7
## 44 0.0705 0.1642 0.1473 0.0839 0.0063 0.3733 0.2564 7
## --------------------------------------------------------
## db$i: 8
## area mc nc ml cll nl larea i
## 42 0.6166 0.9866 1.1290 0.7568 0.6006 0.6236 0.5779 8
## 61 0.6221 0.6047 0.5249 0.7923 0.7492 0.3984 0.5809 8
## 82 0.2902 0.4579 0.7137 0.2610 0.3406 0.3984 0.3933 8
## 6 -0.1613 -0.2764 -0.3436 -0.0931 -0.5880 0.1231 0.0985 8
## 91 -1.7296 -1.5980 -1.2876 -1.8639 -1.2937 -1.5534 -1.7685 8
## 92 -2.1814 -2.1855 -1.9673 -2.9972 -2.0366 -2.6294 -3.2857 8
## 64 -0.0977 -0.2176 -0.2303 0.2610 -0.4394 -0.3523 0.1433 8
## 19 -0.0344 -0.1295 -0.0415 0.4027 0.6378 0.2232 0.1868 8
## 79 -1.6337 -1.8036 -1.8540 -1.7577 -1.8880 -1.8036 -1.5697 8
## 51 0.4642 0.8984 1.0157 0.4381 0.4892 0.4234 0.4943 8
## 69 0.4984 0.4873 0.6004 0.6152 0.5635 1.2742 0.5134 8
## 52 -0.9992 -1.0106 -1.2498 -0.4473 -1.2194 -1.3032 -0.6443 8
## 70 -0.0001 -0.0120 -0.3059 0.2610 -0.1794 -0.0270 0.2098 8# Dado lambda faz: estimação e calcula erro de treino e teste.
myfit_ridge <- function(lambda,
X_train,
y_train,
X_test,
y_test,
XlX = crossprod(X_train),
Xly = crossprod(X_train, y_train)) {
# Estima os parâmetros sujeitos à penalização.
beta <- beta_ridge(Xly, XlX, alpha = lambda)
# Valores ajustados.
fit_train <- X_train %*% beta
fit_test <- X_test %*% beta
# Erro quadrático médio.
mse_train <- crossprod(y_train - fit_train)/nrow(y_train)
mse_test <- crossprod(y_test - fit_test)/nrow(y_test)
return(c(lambda = lambda,
mse_train = mse_train,
mse_test = mse_test))
}
# Vetor de lambda para avaliar o modelo.
lam_seq <- seq(0, 10, by = 0.1)
length(lam_seq)## [1] 101# dput(names(mal))
v <- c("mc", "nc", "ml", "cll", "nl")
cvfit <- lapply(lam_seq,
FUN = function(lambda) {
x <- lapply(1:k,
FUN = function(fold) {
is_test <- db$i == fold
X_train <- model.matrix(~ ., db[!is_test, v])
y_train <- cbind(db[!is_test, "larea"])
X_test <- model.matrix(~ ., db[is_test, v])
y_test <- cbind(db[is_test, "larea"])
XlX = crossprod(X_train)
Xly = crossprod(X_train, y_train)
res <- myfit_ridge(lambda,
X_train,
y_train,
X_test,
y_test,
XlX = XlX,
Xly = Xly)
res <- c(res, fold = fold)
return(res)
})
do.call(rbind, x)
})
cvfit <- data.frame(do.call(rbind, cvfit))
str(cvfit)## 'data.frame': 808 obs. of 4 variables:
## $ lambda : num 0 0 0 0 0 0 0 0 0.1 0.1 ...
## $ mse_train: num 0.0633 0.0584 0.0358 0.0608 0.0658 ...
## $ mse_test : num 0.0399 0.0837 0.267 0.0637 0.0201 ...
## $ fold : num 1 2 3 4 5 6 7 8 1 2 ...xyplot(mse_train + mse_test ~ lambda,
groups = fold,
auto.key = list(title = "Fold",
cex.title = 1,
points = FALSE,
lines = TRUE,
corner = c(0.1, 0.9)),
data = cvfit,
type = "l")
# Obtém o erro médio dos k-folds e o lambda ótimo.
cvfitm <- aggregate(mse_test ~ lambda, data = cvfit, FUN = mean)
lambda_opt <- cvfitm$lambda[which.min(cvfitm$mse_test)]
lambda_opt## [1] 1.2xyplot(mse_test ~ lambda, data = cvfitm, type = "b") +
latticeExtra::layer({
panel.abline(v = lambda_opt, lty = 2)
})
3.4 Ajuste com o pacote glmnet
O pacote glmnet é especializado para o ajuste de modelos de regressão com regualização LASSO, Ridge e elastic net. Ele possui procedimentos extremamente eficientes para isso. Além de modelos de regressão, o pacote acomoda modelos multinomiais, de Poisson, de Cox e regressão logística.
A seguinte vinheta https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html ilustra as principais funcionalidades do pacote. Nesse tutorial, as funcionalidaes serão aplicadas no conjunto de dados de uva. Mais exemplos de uso no glmnet podem ser encontrados nos endereços a seguir.
- https://www4.stat.ncsu.edu/~post/josh/LASSO_Ridge_Elastic_Net_-_Examples.html.
- https://drsimonj.svbtle.com/ridge-regression-with-glmnet.
#-----------------------------------------------------------------------
# glmnet tem recursos para regularização.
y <- m0$model[, 1] # Matriz da resposta.
X <- model.matrix(m0)[, -1] # Matrix das preditoras (sem intercepto).
# Quando alpha = 0, então é Ridge. Se alpha = 1, então é Lasso.
fit <- glmnet(x = X, y = y, alpha = 0, nlambda = 100)
# Traço das estimativas em função do hiperparâmetro de penalização.
plot(fit, xvar = "lambda", label = TRUE)
abline(h = 0, lty = 2)
# str(fit)
# plot(fit$dev.ratio ~ fit$lambda)
# Estimativas dos parâmetros com diferentes níveis de penalidade.
round(coef(fit, s = 0.1), digits = 4)## 6 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.0000
## mc 0.1923
## ml 0.4485
## nc 0.0614
## nl 0.2606
## cll 0.0042round(coef(fit, s = c(0, 0.5, 1, 2, 5)), digits = 4)## 6 x 5 sparse Matrix of class "dgCMatrix"
## 1 2 3 4 5
## (Intercept) 0.0000 0.0000 0.0000 0.0000 0.0000
## mc 0.1934 0.1704 0.1573 0.1352 0.0938
## ml 0.4540 0.2816 0.2198 0.1667 0.1053
## nc 0.0590 0.1204 0.1275 0.1180 0.0858
## nl 0.2627 0.2015 0.1758 0.1453 0.0981
## cll -0.0010 0.1134 0.1262 0.1183 0.0862É importante ressaltar que a diferença no traço das estimativas dos parâmetros bem como no valor que será obtido para o parâmetro de penalização difere do obtido pelo MASS e da implementação manual do método. Isso deve a escolhas na forma de implementar o algorítmo para torná-lo mais eficiente. Existe uma explicação detalhada no Stack Exchange sobre isso: https://stats.stackexchange.com/questions/129179/why-is-glmnet-ridge-regression-giving-me-a-different-answer-than-manual-calculat.
#-----------------------------------------------------------------------
# Validação cruzada para escolha do hiperparâmetro de penalidade.
cvfit <- cv.glmnet(x = X, y = y)
# Perfil do MSE de validação em função de lambda.
plot(cvfit)
abline(v = log(cvfit$lambda.min), col = 2)
# Valor encontrado na escala adotada pelo glmnet.
cvfit$lambda.min## [1] 0.002946199# Estimativas dos parâmetros com o lambda mínimo (cuidado ao
# interpretar).
round(coef(cvfit, s = "lambda.min"), digits = 4)## 6 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.0000
## mc 0.2169
## ml 0.5777
## nc .
## nl 0.3275
## cll -0.1332# Predições.
lam <- cvfit$lambda # Grid de valores de lambda.
Xnew <- X[sample(1:nrow(X), 20), ] # Amostra que pontos observados.
Ypred <- predict(cvfit, s = lam, newx = Xnew) # Predições.
# Traço dos valores preditos.
matplot(log(lam), t(Ypred), type = "l")
abline(h = coef(cvfit, s = 0)[1, ], lty = 2, lwd = 3, col = 2)
O glmnet também permite que sejam passadas partições que definem os folds.
# Mudando o número de folds.
cvfit <- cv.glmnet(x = X, y = y, nfolds = 5)
n <- nrow(X) # Número de registros.
k <- 10 # Número de folds.
i <- ceiling((1:n)/(n/k))
# Passando o identificador de folds.
cvfit <- cv.glmnet(x = X, y = y, foldid = i)3.5 Usando o pacote caret
O pacote caret é um agregador de modelos de machine learning. A principal contribuição do caret, além de agregar as funcionalidades, é padronizar a forma de entrada dos dados e de extração de resultados do ajuste, como predição. Abaixo está ilustrado como o fazer regularização Ridge com o caret que chama internamente o glmnet.
# Uma validação cruzada 5-fold com 3 repetições independentes.
ctrl <- trainControl(method = "repeatedcv",
repeats = 3,
number = 5)
set.seed(123)
# alpha = 0 corresponde à norma L2. lambda é o parâmetro de penalidade.
fit <- train(x = X,
y = y,
method = "glmnet",
trControl = ctrl,
metric = "RMSE",
tuneGrid = expand.grid(alpha = 0,
lambda = seq(0.001, 0.1, by = 0.001)))
fit## glmnet
##
## 100 samples
## 5 predictors
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 3 times)
## Summary of sample sizes: 80, 80, 80, 80, 80, 80, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared MAE
## 0.001 0.2640897 0.9502364 0.1803514
## 0.002 0.2640897 0.9502364 0.1803514
## 0.003 0.2640897 0.9502364 0.1803514
## 0.004 0.2640897 0.9502364 0.1803514
## 0.005 0.2640897 0.9502364 0.1803514
## 0.006 0.2640897 0.9502364 0.1803514
## 0.007 0.2640897 0.9502364 0.1803514
## 0.008 0.2640897 0.9502364 0.1803514
## 0.009 0.2640897 0.9502364 0.1803514
## 0.010 0.2640897 0.9502364 0.1803514
## 0.011 0.2640897 0.9502364 0.1803514
## 0.012 0.2640897 0.9502364 0.1803514
## 0.013 0.2640897 0.9502364 0.1803514
## 0.014 0.2640897 0.9502364 0.1803514
## 0.015 0.2640897 0.9502364 0.1803514
## 0.016 0.2640897 0.9502364 0.1803514
## 0.017 0.2640897 0.9502364 0.1803514
## 0.018 0.2640897 0.9502364 0.1803514
## 0.019 0.2640897 0.9502364 0.1803514
## 0.020 0.2640897 0.9502364 0.1803514
## 0.021 0.2640897 0.9502364 0.1803514
## 0.022 0.2640897 0.9502364 0.1803514
## 0.023 0.2640897 0.9502364 0.1803514
## 0.024 0.2640897 0.9502364 0.1803514
## 0.025 0.2640897 0.9502364 0.1803514
## 0.026 0.2640897 0.9502364 0.1803514
## 0.027 0.2640897 0.9502364 0.1803514
## 0.028 0.2640897 0.9502364 0.1803514
## 0.029 0.2640897 0.9502364 0.1803514
## 0.030 0.2640897 0.9502364 0.1803514
## 0.031 0.2640897 0.9502364 0.1803514
## 0.032 0.2640897 0.9502364 0.1803514
## 0.033 0.2640897 0.9502364 0.1803514
## 0.034 0.2640897 0.9502364 0.1803514
## 0.035 0.2640897 0.9502364 0.1803514
## 0.036 0.2640897 0.9502364 0.1803514
## 0.037 0.2640897 0.9502364 0.1803514
## 0.038 0.2640897 0.9502364 0.1803514
## 0.039 0.2640897 0.9502364 0.1803514
## 0.040 0.2640897 0.9502364 0.1803514
## 0.041 0.2640897 0.9502364 0.1803514
## 0.042 0.2640897 0.9502364 0.1803514
## 0.043 0.2640897 0.9502364 0.1803514
## 0.044 0.2640897 0.9502364 0.1803514
## 0.045 0.2640897 0.9502364 0.1803514
## 0.046 0.2640897 0.9502364 0.1803514
## 0.047 0.2640897 0.9502364 0.1803514
## 0.048 0.2640897 0.9502364 0.1803514
## 0.049 0.2640897 0.9502364 0.1803514
## 0.050 0.2640897 0.9502364 0.1803514
## 0.051 0.2640897 0.9502364 0.1803514
## 0.052 0.2640897 0.9502364 0.1803514
## 0.053 0.2640897 0.9502364 0.1803514
## 0.054 0.2640897 0.9502364 0.1803514
## 0.055 0.2640897 0.9502364 0.1803514
## 0.056 0.2640897 0.9502364 0.1803514
## 0.057 0.2640897 0.9502364 0.1803514
## 0.058 0.2640897 0.9502364 0.1803514
## 0.059 0.2640897 0.9502364 0.1803514
## 0.060 0.2640897 0.9502364 0.1803514
## 0.061 0.2640897 0.9502364 0.1803514
## 0.062 0.2640897 0.9502364 0.1803514
## 0.063 0.2640897 0.9502364 0.1803514
## 0.064 0.2640897 0.9502364 0.1803514
## 0.065 0.2640897 0.9502364 0.1803514
## 0.066 0.2640897 0.9502364 0.1803514
## 0.067 0.2640897 0.9502364 0.1803514
## 0.068 0.2640897 0.9502364 0.1803514
## 0.069 0.2640897 0.9502364 0.1803514
## 0.070 0.2640897 0.9502364 0.1803514
## 0.071 0.2640897 0.9502364 0.1803514
## 0.072 0.2640897 0.9502364 0.1803514
## 0.073 0.2640897 0.9502364 0.1803514
## 0.074 0.2640897 0.9502364 0.1803514
## 0.075 0.2640897 0.9502364 0.1803514
## 0.076 0.2640897 0.9502364 0.1803514
## 0.077 0.2640897 0.9502364 0.1803514
## 0.078 0.2640897 0.9502364 0.1803514
## 0.079 0.2640897 0.9502364 0.1803514
## 0.080 0.2640897 0.9502364 0.1803514
## 0.081 0.2641056 0.9502350 0.1803583
## 0.082 0.2641363 0.9502322 0.1803715
## 0.083 0.2641671 0.9502293 0.1803847
## 0.084 0.2641979 0.9502264 0.1803979
## 0.085 0.2642287 0.9502235 0.1804111
## 0.086 0.2642707 0.9502188 0.1804263
## 0.087 0.2643328 0.9502109 0.1804488
## 0.088 0.2644046 0.9502014 0.1804741
## 0.089 0.2644919 0.9501887 0.1805027
## 0.090 0.2645990 0.9501708 0.1805383
## 0.091 0.2647173 0.9501497 0.1805783
## 0.092 0.2648395 0.9501265 0.1806210
## 0.093 0.2649680 0.9501001 0.1806681
## 0.094 0.2650968 0.9500734 0.1807150
## 0.095 0.2652251 0.9500468 0.1807611
## 0.096 0.2653560 0.9500189 0.1808072
## 0.097 0.2655088 0.9499814 0.1808589
## 0.098 0.2656712 0.9499396 0.1809117
## 0.099 0.2658433 0.9498850 0.1809754
## 0.100 0.2660107 0.9498245 0.1810372
##
## Tuning parameter 'alpha' was held constant at a value of 0
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0 and lambda = 0.08.lam_opt <- fit$finalModel$tuneValue$lambda
lam_opt## [1] 0.08plot(fit) + latticeExtra::layer(panel.abline(v = lam_opt, lty = 2))
# Coeficientes para o lambda ótimo.
round(coef(fit$finalModel, fit$bestTune$lambda), digits = 4)## 6 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.0000
## mc 0.1934
## ml 0.4540
## nc 0.0590
## nl 0.2627
## cll -0.00104 Regularização \(L_1\)
4.1 Usando o pacote penalized
O pacote penalized tem funções primitivas para fazer regressão com regularização LASSO. Então, para ilustrar a implementação do zero, será usado o otimizador deste pacote. Todavia, poderia-se escrever a função objetivo com a penalização de norma L1 e passar para a função optim() estimar os parâmetros sujeito à penalização.
# Exemplo de uso do `penalized` para fazer regressão com regul. LASSO.
m0l <- penalized(response = mal$larea,
penalized = model.matrix(m0)[, -1],
lambda1 = 1,
lambda2 = 0)## # nonzero coefficients: 6
# nonzero coefficients: 6
# nonzero coefficients: 6
# nonzero coefficients: 6
# nonzero coefficients: 6
# nonzero coefficients: 5
# nonzero coefficients: 6
# nonzero coefficients: 6
# nonzero coefficients: 6
# nonzero coefficients: 5
# nonzero coefficients: 5
# nonzero coefficients: 5
# nonzero coefficients: 5
# nonzero coefficients: 5 round(cbind(coef(m0l)), digits = 4)## [,1]
## (Intercept) 0.0000
## mc 0.1705
## ml 0.5561
## nl 0.2848
## cll -0.02554.2 Estudo da função objetivo com penalização \(L_1\)
4.3 Validação cruzada para o parâmetro de penalização com norma \(L_1\)
A seguir estão funções e código para fazer a validação cruzada \(k\)-fold para o parâmetro de penalização com regularização LASSO no conjunto de dados de uva.
#-----------------------------------------------------------------------
# Validação cruzada.
# Dado lambda faz: estimação e calcula erro de treino e teste.
myfit_lasso <- function(lambda,
X_train,
y_train,
X_test,
y_test,
XlX = crossprod(X_train),
Xly = crossprod(X_train, y_train)) {
capture.output({
mod <- penalized(response = y_train,
penalized = X_train,
lambda1 = lambda,
lambda2 = 0)
})
beta <- coef(mod)
n <- names(beta)
fit_train <- X_train[, n] %*% beta
fit_test <- X_test[, n] %*% beta
mse_train <- crossprod(y_train - fit_train)/nrow(y_train)
mse_test <- crossprod(y_test - fit_test)/nrow(y_test)
return(c(lambda = lambda,
mse_train = mse_train,
mse_test = mse_test))
}
# Vetor de lambda para avaliar o modelo.
lam_seq <- seq(0.1, 4, length.out = 100)
length(lam_seq)## [1] 100# dput(names(mal))
v <- c("mc", "nc", "ml", "cll", "nl")
cvfit <- lapply(lam_seq,
FUN = function(lambda) {
x <- lapply(1:k,
FUN = function(fold) {
is_test <- db$i == fold
X_train <- model.matrix(~ ., db[!is_test, v])
y_train <- cbind(db[!is_test, "larea"])
X_test <- model.matrix(~ ., db[is_test, v])
y_test <- cbind(db[is_test, "larea"])
XlX = crossprod(X_train)
Xly = crossprod(X_train, y_train)
res <- myfit_lasso(lambda,
X_train,
y_train,
X_test,
y_test,
XlX = XlX,
Xly = Xly)
res <- c(res, fold = fold)
return(res)
})
do.call(rbind, x)
})
cvfit <- data.frame(do.call(rbind, cvfit))
str(cvfit)## 'data.frame': 800 obs. of 4 variables:
## $ lambda : num 0.1 0.1 0.1 0.1 0.1 ...
## $ mse_train: num 0.0634 0.0585 0.0358 0.0609 0.0659 ...
## $ mse_test : num 0.0389 0.0859 0.2665 0.0629 0.0201 ...
## $ fold : num 1 2 3 4 5 6 7 8 1 2 ...xyplot(mse_train + mse_test ~ lambda,
groups = fold,
auto.key = list(space = "right"),
data = cvfit,
type = "l")
xyplot(mse_train + mse_test ~ lambda,
outer = TRUE,
auto.key = TRUE,
data = cvfit,
cex = 0.4,
type = c("p", "a"))
# Obtém o erro médio dos k-folds e o lambda ótimo.
cvfitm <- aggregate(mse_test ~ lambda, data = cvfit, FUN = mean)
lambda_opt <- cvfitm$lambda[which.min(cvfitm$mse_test)]
lambda_opt## [1] 0.2181818xyplot(mse_test ~ lambda, data = cvfitm, type = "b") +
latticeExtra::layer({
panel.abline(v = lambda_opt)
})
4.4 Ajuste com o pacote glmnet
y <- m0$model[, 1] # Matriz da resposta.
X <- model.matrix(m0)[, -1] # Matrix das preditoras (sem intercepto).
# Quando alpha = 0, então é Ridge. Se alpha = 1, então é Lasso.
fit <- glmnet(x = X, y = y, alpha = 1, nlambda = 100)
# Traço das estimativas em função do hiperparâmetro de penalização.
plot(fit, xvar = "lambda", label = TRUE)
abline(h = 0, lty = 2)
# Estimativas dos parâmetros com diferentes níveis de penalidade.
round(coef(fit, s = 0.1), digits = 4)## 6 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.0000
## mc 0.1218
## ml 0.5075
## nc .
## nl 0.2609
## cll .round(coef(fit, s = c(0, 0.1, 0.25, 0.5, 1, 2, 5)), digits = 4)## 6 x 7 sparse Matrix of class "dgCMatrix"
## 1 2 3 4 5 6 7
## (Intercept) 0.0000 0.0000 0.0000 0.0000 0 0 0
## mc 0.2273 0.1218 0.0515 . . . .
## ml 0.5826 0.5075 0.4315 0.3011 . . .
## nc . . . . . . .
## nl 0.3398 0.2609 0.2454 0.1602 . . .
## cll -0.1601 . . . . . .#-----------------------------------------------------------------------
# Validação cruzada para escolha do hiperparâmetro de penalidade.
cvfit <- cv.glmnet(x = X, y = y)
# Perfil do MSE de validação em função de lambda.
plot(cvfit)
cvfit$lambda.min## [1] 0.002445986abline(v = log(cvfit$lambda.min), col = 2)
# Estimativas dos parâmetros com o lambda mínimo (cuidado ao
# interpretar).
round(coef(cvfit, s = "lambda.min"), digits = 4)## 6 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.0000
## mc 0.2198
## ml 0.5791
## nc .
## nl 0.3309
## cll -0.1407# Predições.
Xnew <- X[sample(1:nrow(X), 20), ] # Amostra que pontos observados.
lam <- cvfit$lambda # Grid de valores de lambda.
Ypred <- predict(cvfit, s = lam, newx = Xnew) # Predições.
# Traço dos valores preditos.
matplot(log(lam), t(Ypred), type = "l")
abline(h = coef(cvfit, s = 0)[1, ], lty = 2, lwd = 3, col = 2)
4.5 Usando o pacote caret
# Uma validação cruzada 5-fold com 3 repetições independentes.
ctrl <- trainControl(method = "repeatedcv",
repeats = 3,
number = 5)
set.seed(123)
# alpha = 0 corresponde à norma L2. lambda é o parâmetro de penalidade.
fit <- train(x = X,
y = y,
method = "glmnet",
trControl = ctrl,
metric = "RMSE",
tuneGrid = expand.grid(alpha = 1,
lambda = seq(0.001, 0.1, by = 0.001)))
fit## glmnet
##
## 100 samples
## 5 predictors
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 3 times)
## Summary of sample sizes: 80, 80, 80, 80, 80, 80, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared MAE
## 0.001 0.2689424 0.9431961 0.1871369
## 0.002 0.2655212 0.9448401 0.1851936
## 0.003 0.2627679 0.9461590 0.1833364
## 0.004 0.2612253 0.9469674 0.1820804
## 0.005 0.2603014 0.9474592 0.1810554
## 0.006 0.2603068 0.9476062 0.1806836
## 0.007 0.2605704 0.9476618 0.1805943
## 0.008 0.2607934 0.9477494 0.1805033
## 0.009 0.2611478 0.9477722 0.1804434
## 0.010 0.2616267 0.9477289 0.1804782
## 0.011 0.2621542 0.9476119 0.1805497
## 0.012 0.2625175 0.9475173 0.1804679
## 0.013 0.2625642 0.9475051 0.1802253
## 0.014 0.2624751 0.9475384 0.1798594
## 0.015 0.2624005 0.9475704 0.1795477
## 0.016 0.2623475 0.9475795 0.1792611
## 0.017 0.2623022 0.9475760 0.1789747
## 0.018 0.2622603 0.9475726 0.1786892
## 0.019 0.2622195 0.9475704 0.1784017
## 0.020 0.2621815 0.9475683 0.1781130
## 0.021 0.2621461 0.9475665 0.1778235
## 0.022 0.2621142 0.9475653 0.1775380
## 0.023 0.2620885 0.9475626 0.1772603
## 0.024 0.2620662 0.9475593 0.1769948
## 0.025 0.2620460 0.9475563 0.1767341
## 0.026 0.2620292 0.9475533 0.1764833
## 0.027 0.2620161 0.9475501 0.1762375
## 0.028 0.2620054 0.9475475 0.1759908
## 0.029 0.2619975 0.9475454 0.1757433
## 0.030 0.2619935 0.9475430 0.1755002
## 0.031 0.2619940 0.9475396 0.1752681
## 0.032 0.2619973 0.9475360 0.1750364
## 0.033 0.2620036 0.9475324 0.1748117
## 0.034 0.2620126 0.9475292 0.1745911
## 0.035 0.2620255 0.9475257 0.1743737
## 0.036 0.2620417 0.9475223 0.1741561
## 0.037 0.2620614 0.9475189 0.1739418
## 0.038 0.2620834 0.9475163 0.1737357
## 0.039 0.2621086 0.9475137 0.1735302
## 0.040 0.2621368 0.9475118 0.1733248
## 0.041 0.2621678 0.9475102 0.1731193
## 0.042 0.2622058 0.9475063 0.1729171
## 0.043 0.2622471 0.9475024 0.1727146
## 0.044 0.2622925 0.9474975 0.1725122
## 0.045 0.2623415 0.9474924 0.1723115
## 0.046 0.2623936 0.9474875 0.1721148
## 0.047 0.2624486 0.9474827 0.1719182
## 0.048 0.2625068 0.9474784 0.1717217
## 0.049 0.2625678 0.9474744 0.1715271
## 0.050 0.2626324 0.9474702 0.1713404
## 0.051 0.2627006 0.9474658 0.1711582
## 0.052 0.2627721 0.9474613 0.1709903
## 0.053 0.2628476 0.9474559 0.1708245
## 0.054 0.2629267 0.9474503 0.1706638
## 0.055 0.2630091 0.9474447 0.1705031
## 0.056 0.2630949 0.9474390 0.1703449
## 0.057 0.2631840 0.9474333 0.1701891
## 0.058 0.2632765 0.9474277 0.1700333
## 0.059 0.2633723 0.9474220 0.1698775
## 0.060 0.2634717 0.9474164 0.1697266
## 0.061 0.2635746 0.9474107 0.1695822
## 0.062 0.2636809 0.9474050 0.1694378
## 0.063 0.2637902 0.9473995 0.1692931
## 0.064 0.2639025 0.9473944 0.1691482
## 0.065 0.2640177 0.9473895 0.1690029
## 0.066 0.2641361 0.9473846 0.1688576
## 0.067 0.2642576 0.9473799 0.1687196
## 0.068 0.2643823 0.9473751 0.1685857
## 0.069 0.2645106 0.9473701 0.1684545
## 0.070 0.2646427 0.9473646 0.1683253
## 0.071 0.2647784 0.9473590 0.1681993
## 0.072 0.2649178 0.9473533 0.1680735
## 0.073 0.2650607 0.9473475 0.1679478
## 0.074 0.2652072 0.9473416 0.1678222
## 0.075 0.2653571 0.9473356 0.1676967
## 0.076 0.2655100 0.9473299 0.1675709
## 0.077 0.2656662 0.9473243 0.1674506
## 0.078 0.2658259 0.9473187 0.1673362
## 0.079 0.2659891 0.9473129 0.1672223
## 0.080 0.2661556 0.9473072 0.1671152
## 0.081 0.2663255 0.9473013 0.1670091
## 0.082 0.2664989 0.9472955 0.1669027
## 0.083 0.2666760 0.9472893 0.1667964
## 0.084 0.2668567 0.9472829 0.1666928
## 0.085 0.2670408 0.9472764 0.1665948
## 0.086 0.2672279 0.9472701 0.1664965
## 0.087 0.2674186 0.9472637 0.1663983
## 0.088 0.2676127 0.9472572 0.1663028
## 0.089 0.2678102 0.9472508 0.1662126
## 0.090 0.2680111 0.9472442 0.1661226
## 0.091 0.2682153 0.9472378 0.1660328
## 0.092 0.2684228 0.9472315 0.1659429
## 0.093 0.2686338 0.9472252 0.1658531
## 0.094 0.2688485 0.9472187 0.1657644
## 0.095 0.2690666 0.9472122 0.1656791
## 0.096 0.2692881 0.9472056 0.1655952
## 0.097 0.2695131 0.9471990 0.1655191
## 0.098 0.2697416 0.9471924 0.1654430
## 0.099 0.2699736 0.9471857 0.1653685
## 0.100 0.2702090 0.9471788 0.1652958
##
## Tuning parameter 'alpha' was held constant at a value of 1
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 1 and lambda = 0.005.lam_opt <- fit$finalModel$tuneValue$lambda
lam_opt## [1] 0.005plot(fit) + latticeExtra::layer(panel.abline(v = lam_opt, lty = 2))
# Coeficientes para o lambda ótimo.
round(coef(fit$finalModel, fit$bestTune$lambda), digits = 4)## 6 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.0000
## mc 0.2049
## ml 0.5720
## nc .
## nl 0.3133
## cll -0.1021| Machine Learning | Prof. Eduardo V. Ferreira & Prof. Walmes M. Zeviani |
