Commit dc44bbf9 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Adiciona tutorial de regularização.

parent 06eebccd
......@@ -27,6 +27,8 @@ navbar:
href: tutorials/01-cross-validation.html
- text: "Gradiente descendente"
href: tutorials/02-gradient-methods.html
- text: "Regularização"
href: tutorials/03-regularization.html
- text: "Scripts"
icon: fa-file-text
href: scripts/
......
---
title: "Regularização"
author: Prof. Walmes M. Zeviani & Prof. Eduardo V. Ferreira
date: 2017-10-19
#bibliography: ../config/Refs.bib
#csl: ../config/ABNT-UFPR-2011-Mendeley.csl
---
```{r, include = FALSE}
source("../config/setup.R")
opts_chunk$set(
cache = FALSE,
message = FALSE,
warning = FALSE)
```
# Métodos de seleção de variáveis
```{r, message = FALSE}
rm(list = objects())
library(lattice)
library(latticeExtra)
```
```{r}
# Endereço das tabelas.
pre <- "https://raw.githubusercontent.com/walmes/EACS/master/data-raw/"
files <- c(hidric = "teca_crapar.csv",
quimic = "teca_qui.csv",
prod = "teca_arv.csv")
urls <- paste0(pre, files)
names(urls) <- names(files)
# Lista com as tabelas.
da <- sapply(urls,
FUN = read.table,
header = TRUE,
sep = ";",
simplify = FALSE)
str(da)
# Manipular as tabelas para fazer o merge.
da$quimic <- subset(da$quimic, cam == "[0, 5)", select = -cam)
da$hidric <- subset(da$hidric, cam == "[0, 5)", select = -c(cam, cad))
da$prod <- subset(da$prod, select = c(loc, prod))
str(da)
# Aplica o merge recursivamente.
db <- Reduce(f = merge, x = da)
# Elimina a variável identificadora (agora desnecessária).
db$loc <- NULL
# Estrutura da tabela.
str(db)
# Padronizando as variáveis.
db[, -ncol(db)] <- sapply(db[, -ncol(db)],
FUN = scale,
center = TRUE,
scale = TRUE)
summary(db)
# Gráfico de pares de diagrama de dispersão.
splom(db, as.matrix = TRUE, type = c("p", "r"))
```
## Melhor subconjunto de variáveis
```{r}
# Ajuste do modelo com todas as variáveis.
m0 <- lm(prod ~ ., data = db)
summary(m0)
# Quantos subconjuntos podem ser formados de tamanho 1 até k?
k <- ncol(db) - 1
# É a quantidade total de modelos!
sum(sapply(1:k, choose, n = k))
# Recursos para avaliar todas as combinações possíveis.
library(leaps)
# help(regsubsets, h = "html")
# Melhor modelo de cada tamanho.
b0 <- regsubsets(prod ~ .,
data = db,
method = "exhaustive",
nvmax = 12)
sel <- summary(b0)
sel
# Extrai as medidas de ajuste.
mea <- sapply(c("rsq", "rss", "adjr2", "cp", "bic"),
FUN = function(i) sel[[i]])
mea <- cbind(mod = 1:nrow(mea), as.data.frame(mea))
mea
# Exibe as medidas de ajuste.
f <- paste(paste(colnames(mea)[-1], collapse = " + "), " ~ mod")
xyplot(as.formula(f),
data = mea,
outer = TRUE,
scales = "free",
type = "h",
lwd = 3,
as.table = TRUE,
xlab = "Modelo",
ylab = "Medida de ajuste")
str(sel)
# Quais as variáveis do melhor modelo?
v <- names(which(sel$which[2, ])[-1])
v
summary(lm(prod ~ ., data = db[, c("prod", v)]))
```
## Métodos passo a passo (stepwise)
```{r}
summary(m0)
# Critério de seleção AIC: k = 2.
m1 <- step(m0, direction = "both")
summary(m1)
# Critério de seleção BIC: k = log(n).
m2 <- step(m0, direction = "both", k = log(nrow(db)))
summary(m2)
```
# Métodos de regularização de variáveis
## Regressão Ridge
```{r}
#-----------------------------------------------------------------------
# Usando o pacote MASS.
library(MASS)
m3 <- lm.ridge(prod ~ .,
data = db,
lambda = seq(from = 0, to = 0.5, by = 0.01))
# head(m3)
# Traço das estimativas como função de lambda.
plot(m3)
# Seleção do melhor valor.
select(m3)
#-----------------------------------------------------------------------
# Implementando a regressão ridge (vanilla flavor).
rid <- function(lambda, y, X) {
XlX <- crossprod(X)
Xly <- crossprod(X, y)
II <- diag(ncol(X))
II[1, 1] <- 0
beta <- solve(XlX + lambda * II, Xly)
return(beta)
}
names(db)
# Vetor resposta e matriz do modelo.
y <- cbind(db$prod)
X <- model.matrix(~ ph + ctc + mg + I + are, data = db)
# X <- model.matrix(~ . - prod, data = db)
# ATTENTION: tem algum detalhe ainda não resolvido aqui.
# rid(lambda = 0.5, y = y, X = X)
# lm.ridge(prod ~ ph + ctc + ca + mg + S,
# data = db,
# lambda = 0.5)
# Sequencia de valores de lambda e estimativas dos parâmetros.
l <- seq(0, 50, length.out = 31)
r <- matrix(0, nrow = ncol(X), ncol = length(l))
for (i in seq_along(l)) {
r[, i] <- rid(lambda = l[i], y = y, X = X)
}
# Traço.
matplot(l,
t(r[-1, ]),
type = "o",
lty = 1,
pch = 1,
xlab = expression(lambda),
ylab = expression(hat(beta)))
abline(h = 0, lty = 2)
legend("topright",
legend = colnames(X),
col = palette()[1:ncol(X)],
lty = 1,
bty = "n")
```
## Regressão Lasso
```{r}
#-----------------------------------------------------------------------
# Implementação baunilha.
las <- function(beta, lambda, y, X) {
# Valores preditos.
ypred <- X %*% beta
# Sum of squares.
rss <- sum((y - ypred)^2)
# Penalidade no tamanho dos coeficientes.
penalty <- lambda * sum(abs(beta[-1]))
return(rss + penalty)
}
# Valores iniciais.
start <- coef(lm(y ~ 0 + X))
start
# Otimiza nos parâmetros para minimizar a função objetivo.
op <- optim(par = start,
fn = las,
lambda = 0,
y = y,
X = X)
op$par
# Avalia uma sequência de valores de lambda.
l <- seq(from = 0, to = 250, length.out = 50)
r <- sapply(l,
FUN = function(lam) {
op <- optim(start,
fn = las,
lambda = lam,
y = y,
X = X)
start <<- op$par
return(op$par)
})
# Traço.
matplot(l,
t(r[-1, ]),
type = "o",
lty = 1,
pch = 1,
xlab = expression(lambda),
ylab = expression(hat(beta)))
abline(h = 0, lty = 2)
legend("right",
legend = colnames(X),
col = palette()[1:ncol(X)],
lty = 1,
bty = "n")
```
## Usando o pacote `glmnet`
```{r}
library(glmnet)
y <- cbind(db$prod)
X <- model.matrix(~ . - prod, data = db)
# X <- model.matrix(~ ph + ctc + ca + mg + S, data = db)
# X <- model.matrix(~ ph + ctc, data = db)
# Ajuste com penalização lasso, ridge e elastic net.
mlasso <- glmnet(x = X, y = y, family = "gaussian", alpha = 1)
mridge <- glmnet(x = X, y = y, family = "gaussian", alpha = 0)
melnet <- glmnet(x = X, y = y, family = "gaussian", alpha = 0.5)
# Traços.
par(mfrow = c(1, 3))
plot(mlasso, xvar = "lambda", main = "LASSO")
plot(mridge, xvar = "lambda", main = "RIDGE")
plot(melnet, xvar = "lambda", main = "ELNET")
layout(1)
```
```{r}
#-----------------------------------------------------------------------
# Validão cruzada para a escolha do lambda.
cv.ridge <- cv.glmnet(X,
y,
family = "gaussian",
alpha = 0,
parallel = TRUE,
type.measure = "mse")
plot(cv.ridge)
cv.ridge$lambda.min
cv.ridge$lambda.1se
coef(cv.ridge, s = cv.ridge$lambda.min)
cv.lasso <- cv.glmnet(X,
y,
family = "gaussian",
alpha = 1,
parallel = TRUE,
type.measure = "mse")
plot(cv.lasso)
cv.lasso$lambda.min
cv.lasso$lambda.1se
coef(cv.lasso, s = cv.lasso$lambda.min)
par(mfrow = c(1, 2))
plot(mridge, xvar = "lambda", main = "RIDGE")
abline(v = log(cv.ridge$lambda.min), lty = 2)
plot(mlasso, xvar = "lambda", main = "LASSO")
abline(v = log(cv.lasso$lambda.min), lty = 2)
layout(1)
```
# Links úteis
* <https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html>;
* <http://www4.stat.ncsu.edu/~post/josh/LASSO_Ridge_Elastic_Net_-_Examples.html>;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment