Commit 2dab6d15 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Aprimora regressão e regularização parte 1.

parent c4b136bf
...@@ -16,15 +16,31 @@ opts_chunk$set( ...@@ -16,15 +16,31 @@ opts_chunk$set(
# Métodos de seleção de variáveis # Métodos de seleção de variáveis
Nesse tutorial serão apresentadas técnicas para seleção de variáveis em
modelos de regressão múltipla. Serão considerados os métodos de todos
os subconjuntos possíveis, seleção passo a passo de variáveis (stepwise)
e regularização.
```{r, message = FALSE} ```{r, message = FALSE}
rm(list = objects()) rm(list = objects())
library(lattice) library(lattice)
library(latticeExtra) library(latticeExtra)
``` ```
* `teca_crapar.csv` contém os parâmetros estimados da curva de
retenção de água do solo.
* `teca_qui.csv` contém valores das variáveis qúmicas do solo.
* `teca_arv.csv` contém os valores de produção de madeira dos nas
localidades onde as variáveis preditoras foram determinadas.
O objetivo destes dados é fazer um modelo para predição da produção de
madeira em função do conjunto de preditoras. A tabela de dados não é
grande pois contém apenas 50 registros.
```{r} ```{r}
# Endereço das tabelas. # Endereço das tabelas.
pre <- "https://raw.githubusercontent.com/walmes/EACS/master/data-raw/" pre <- "https://raw.githubusercontent.com/walmes/EACS/master/data-raw/"
files <- c(hidric = "teca_crapar.csv", files <- c(hidric = "teca_crapar.csv",
quimic = "teca_qui.csv", quimic = "teca_qui.csv",
prod = "teca_arv.csv") prod = "teca_arv.csv")
...@@ -65,8 +81,14 @@ summary(db) ...@@ -65,8 +81,14 @@ summary(db)
splom(db, as.matrix = TRUE, type = c("p", "r")) splom(db, as.matrix = TRUE, type = c("p", "r"))
``` ```
Ao todo são `r nrow(db)` registros completos e `r ncol(db)` variáveis.
## Melhor subconjunto de variáveis ## Melhor subconjunto de variáveis
A técnica do melhor subconjunto de variáveis visa determinar, para uma
quantidade fixa de variáveis, o melhor subconjunto para a predição da
resposta.
```{r} ```{r}
# Ajuste do modelo com todas as variáveis. # Ajuste do modelo com todas as variáveis.
m0 <- lm(prod ~ ., data = db) m0 <- lm(prod ~ ., data = db)
...@@ -74,16 +96,25 @@ summary(m0) ...@@ -74,16 +96,25 @@ summary(m0)
# Quantos subconjuntos podem ser formados de tamanho 1 até k? # Quantos subconjuntos podem ser formados de tamanho 1 até k?
k <- ncol(db) - 1 k <- ncol(db) - 1
k
# É a quantidade total de modelos! # É a quantidade total de modelos! Mais de 1 milhão com 20 variáveis.
sum(sapply(1:k, choose, n = k)) sum(sapply(1:k, choose, n = k))
```
### Ajuste com o pacote `leaps`
```{r}
#-----------------------------------------------------------------------
# Recursos para avaliar todas as combinações possíveis. # Recursos para avaliar todas as combinações possíveis.
library(leaps) library(leaps)
ls("package:leaps")
# help(regsubsets, h = "html") # help(regsubsets, h = "html")
# Melhor modelo de cada tamanho. # Melhor modelo de cada tamanho.
# `nvmax` é o limite superior para o tamanho do modelo.
b0 <- regsubsets(prod ~ ., b0 <- regsubsets(prod ~ .,
data = db, data = db,
method = "exhaustive", method = "exhaustive",
...@@ -118,6 +149,99 @@ v ...@@ -118,6 +149,99 @@ v
summary(lm(prod ~ ., data = db[, c("prod", v)])) summary(lm(prod ~ ., data = db[, c("prod", v)]))
``` ```
### Aplicação de validação cruzada
A validação cruzada será usada para selecionar o número ideal de
variáveis para a predição. Dessa forma, a função custo será o erro
quadrático médio (MSE) da predição no conjunto de teste.
IMPORTANTE: Em cada fold gerado pelas partições o conjunto das variáveis
selecionadas não necessariamente é o mesmo. Ou seja, no fold 1, as
melhores 3 variáveis podem ser `x1`, `x4` e `x8` e no fold 2 podem ser
`x4`, `x7` e `x10`. Isso deve acontecer principalmente entre as
variáveis preditoras mais correlacionadas. No entanto, isso não é um
problema porque a validação cruzada será feita para selecionar o número
ideal de variáveis. Depois será feita a determinação de quais são essas
variáveis.
Aqui será usando validação cruzada 5-fold com 10 repetições
independentes.
```{r}
# Para fazer a predição no conjunto de teste com o ajuste do conjunto de
# treino que considera as variáveis selecionadas.
mse_fit_reg <- function(variables, dt_train, dt_test) {
# print(variables)
# print(str(train[, c("prod", variables)]))
# Ajusta o modelo com os dados de treino.
tra <- dt_train[, c("prod", variables)]
m0 <- lm(prod ~ ., data = tra)
# Predição com os dados de teste.
fitted_vals <- predict(m0, newdata = dt_test[, c("prod", variables)])
# Mean esquare error.
crossprod((dt_test$prod - fitted_vals))/nrow(dt_test)
}
# Ajusta o modelo em cada fold e retorna a medida de ajuste.
fit_fold <- function(fold, indexes) {
# Vetor lógico para particionar treino/teste.
is_test <- indexes == fold
dt_test <- db[is_test, ] # Teste.
dt_train <- db[!is_test, ] # Treino.
# Os melhores ajuste com tamanho máximo de 12 variáveis.
fit <- regsubsets(prod ~ .,
data = dt_train,
method = "exhaustive",
nvmax = 12)
sel <- summary(fit)
# Variáveis mantidas em cada ajuste.
variables <- lapply(apply(sel$which[, -1],
MARGIN = 1,
FUN = which),
FUN = names)
# Medida da função custo nos dados de teste.
mse_test <- sapply(variables,
FUN = mse_fit_reg,
dt_train = dt_train,
dt_test = dt_test)
return(mse_test)
}
n <- nrow(db) # Número de observações.
k <- 5 # Número de folds
i <- ceiling((1:n)/(n/k)) # Indicador de fold
res <- replicate(n = 10, simplify = FALSE, {
ii <- sample(i) # Embaralha os índices.
res <- lapply(1:k, fit_fold, indexes = ii)
res <- as.data.frame(do.call(cbind, res))
res$vars <- 1:nrow(res)
res
})
res <- do.call(rbind, res)
xyplot(V1 + V2 + V3 + V4 + V5 ~ vars,
data = res,
xlab = "Número de variáveis no modelo",
ylab = "Função custo no conjunto de teste",
type = c("p", "a")) +
layer({
m <- aggregate(y ~ x, FUN = mean)
panel.lines(m$x, m$y, col = 1, lwd = 2)
})
# Quais são as melhores variáveis então?
fit <- regsubsets(prod ~ .,
data = db,
method = "exhaustive",
nvmax = 2)
sel <- summary(fit)
v <- names(which(sel$which[2, ])[-1])
summary(lm(prod ~ .,
data = db[, c("prod", v)]))
```
## Métodos passo a passo (stepwise) ## Métodos passo a passo (stepwise)
```{r} ```{r}
...@@ -136,6 +260,9 @@ summary(m2) ...@@ -136,6 +260,9 @@ summary(m2)
## Regressão Ridge ## Regressão Ridge
O pacote `MASS` tem a função `lm.ridge()` para ajuste de modelos de
regressão com regularização de norma $L_2$ ou regressão ridge.
```{r} ```{r}
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Usando o pacote MASS. # Usando o pacote MASS.
...@@ -144,15 +271,51 @@ library(MASS) ...@@ -144,15 +271,51 @@ library(MASS)
m3 <- lm.ridge(prod ~ ., m3 <- lm.ridge(prod ~ .,
data = db, data = db,
lambda = seq(from = 0, to = 0.5, by = 0.01)) lambda = seq(from = 0,
to = 100,
by = 0.1))
# head(m3) # head(m3)
# Traço das estimativas como função de lambda. # Traço das estimativas como função de lambda.
plot(m3) plot(m3)
mtext(side = 1, line = 2, "Valor do parâmetro de regularização")
mtext(side = 2, line = 2, "Valores estimados para os parâmetros")
# Seleção do melhor valor. # Seleção do melhor valor.
select(m3) select(m3)
```
ATENÇÃO! As variáveis preditoras não foram transformadas para a mesma
escala. Os resultados podem ser diferentes.
```{r}
#-----------------------------------------------------------------------
db_norm <- db
i <- names(db) != "prod"
db_norm[, i] <- as.data.frame(scale(db[, i]))
m3 <- lm.ridge(prod ~ .,
data = db_norm,
lambda = seq(from = 0,
to = 100,
by = 0.1))
# Traço das estimativas como função de lambda.
plot(m3)
mtext(side = 1, line = 2, "Valor do parâmetro de regularização")
mtext(side = 2, line = 2, "Valores estimados para os parâmetros")
# Seleção do melhor valor.
select(m3)
```
Pelos testes, o parâmetro de regularização escolhido por GCV foi o mesmo
com a normalização ou não das preditoras. Pode-se acreditar então a
função `lm.ridge()` faça isso internamente. No entanto, a respectiva
documentação não diz nada a respeito.
```{r}
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Implementando a regressão ridge (vanilla flavor). # Implementando a regressão ridge (vanilla flavor).
...@@ -165,44 +328,47 @@ rid <- function(lambda, y, X) { ...@@ -165,44 +328,47 @@ rid <- function(lambda, y, X) {
return(beta) return(beta)
} }
names(db)
# Vetor resposta e matriz do modelo. # Vetor resposta e matriz do modelo.
y <- cbind(db$prod) y <- cbind(db_norm$prod)
X <- model.matrix(~ ph + ctc + mg + I + are, data = db) X <- model.matrix(~ ph + ctc + mg + I + are, data = db_norm)
# X <- model.matrix(~ . - prod, data = db) # X <- model.matrix(~ . - prod, data = db_norm)
# ATTENTION: tem algum detalhe ainda não resolvido aqui. # Valores bem próximos. Diferenças são de sofisticação de implementação.
# rid(lambda = 0.5, y = y, X = X) rid(lambda = 0.5, y = y, X = X)
# lm.ridge(prod ~ ph + ctc + ca + mg + S, lm.ridge(prod ~ .,
# data = db, data = db_norm[, c("prod", colnames(X)[-1])],
# lambda = 0.5) lambda = 0.5)
# Sequencia de valores de lambda e estimativas dos parâmetros. # Sequencia de valores de lambda e estimativas dos parâmetros.
l <- seq(0, 50, length.out = 31) l <- 10^seq(-3, 3, length.out = 200)
r <- matrix(0, nrow = ncol(X), ncol = length(l)) r <- matrix(0, nrow = ncol(X), ncol = length(l))
for (i in seq_along(l)) { for (i in seq_along(l)) {
r[, i] <- rid(lambda = l[i], y = y, X = X) r[, i] <- rid(lambda = l[i], y = y, X = X)
} }
# Traço. # Traço.
matplot(l, matplot(log10(l),
t(r[-1, ]), t(r[-1, ]),
type = "o", type = "o",
lty = 1, lty = 1,
pch = 1, pch = 1,
xlab = expression(lambda), cex = 0.5,
xlab = expression(log(lambda)),
ylab = expression(hat(beta))) ylab = expression(hat(beta)))
abline(h = 0, lty = 2) abline(h = 0, lty = 2)
legend("topright", legend("topright",
legend = colnames(X), legend = colnames(X)[-1],
col = palette()[1:ncol(X)], col = palette()[1:(ncol(X) - 1)],
lty = 1, lty = 1,
bty = "n") bty = "n")
``` ```
## Regressão Lasso ## Regressão Lasso
Para regressão LASSO será utilizado um otimizador geral da função custo
que contém a penalidade de norma $L_1$. O pacote `penalized` tem funções
para o ajuste de regressão LASSO.
```{r} ```{r}
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Implementação baunilha. # Implementação baunilha.
...@@ -225,17 +391,19 @@ start ...@@ -225,17 +391,19 @@ start
op <- optim(par = start, op <- optim(par = start,
fn = las, fn = las,
lambda = 0, lambda = 0,
method = "BFGS",
y = y, y = y,
X = X) X = X)
op$par op$par
# Avalia uma sequência de valores de lambda. # Avalia uma sequência de valores de lambda.
l <- seq(from = 0, to = 250, length.out = 50) l <- 10^seq(-3, 3, length.out = 200)
r <- sapply(l, r <- sapply(l,
FUN = function(lam) { FUN = function(lam) {
op <- optim(start, op <- optim(start,
fn = las, fn = las,
lambda = lam, lambda = lam,
method = "BFGS",
y = y, y = y,
X = X) X = X)
start <<- op$par start <<- op$par
...@@ -243,77 +411,80 @@ r <- sapply(l, ...@@ -243,77 +411,80 @@ r <- sapply(l,
}) })
# Traço. # Traço.
matplot(l, matplot(log10(l),
t(r[-1, ]), t(r[-1, ]),
type = "o", type = "o",
lty = 1, lty = 1,
pch = 1, pch = 1,
xlab = expression(lambda), cex = 0.5,
xlab = expression(log(lambda)),
ylab = expression(hat(beta))) ylab = expression(hat(beta)))
abline(h = 0, lty = 2) abline(h = 0, lty = 2)
legend("right", legend("right",
legend = colnames(X), legend = colnames(X)[-1],
col = palette()[1:ncol(X)], col = palette()[1:(ncol(X) - 1)],
lty = 1, lty = 1,
bty = "n") bty = "n")
``` ```
## Usando o pacote `glmnet` ## Usando o pacote `glmnet`
O pacote `glmnet` é o recomendado para a aplicação de regularização em
modelos de regressão. O pacote permite fazer regularização com
elastic-net da qual LASSO e Ridge são casos extremos quando o parâmetro
$\alpha$ é 1 ou 0, respectivamete.
```{r} ```{r}
library(glmnet) library(glmnet)
y <- cbind(db$prod) y <- cbind(db$prod)
X <- model.matrix(~ . - prod, data = db) X <- model.matrix(~ . - prod, data = db_norm)
X <- X[, -1] # Remove a coluna do intercepto.
# X <- model.matrix(~ ph + ctc + ca + mg + S, data = db) # X <- model.matrix(~ ph + ctc + ca + mg + S, data = db)
# X <- model.matrix(~ ph + ctc, data = db) # X <- model.matrix(~ ph + ctc, data = db)
# Ajuste com penalização lasso, ridge e elastic net. # Ajuste com penalização lasso, ridge e elastic net de mistura 1:1.
mlasso <- glmnet(x = X, y = y, family = "gaussian", alpha = 1)
mridge <- glmnet(x = X, y = y, family = "gaussian", alpha = 0) mridge <- glmnet(x = X, y = y, family = "gaussian", alpha = 0)
melnet <- glmnet(x = X, y = y, family = "gaussian", alpha = 0.5) melnet <- glmnet(x = X, y = y, family = "gaussian", alpha = 0.5)
mlasso <- glmnet(x = X, y = y, family = "gaussian", alpha = 1)
# Traços. # Traços.
par(mfrow = c(1, 3)) par(mfrow = c(1, 3))
plot(mlasso, xvar = "lambda", main = "LASSO")
plot(mridge, xvar = "lambda", main = "RIDGE") plot(mridge, xvar = "lambda", main = "RIDGE")
abline(h = 0, lty = 2)
plot(melnet, xvar = "lambda", main = "ELNET") plot(melnet, xvar = "lambda", main = "ELNET")
abline(h = 0, lty = 2)
plot(mlasso, xvar = "lambda", main = "LASSO")
abline(h = 0, lty = 2)
layout(1) layout(1)
```
```{r}
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Validão cruzada para a escolha do lambda. # Obtendo o valor do parâmetro de penalização.
cv.ridge <- cv.glmnet(X, cvfit_ridge <- cv.glmnet(x = X, y = y, alpha = 0, nfolds = 10, type.measure = "mse")
y, cvfit_lasso <- cv.glmnet(x = X, y = y, alpha = 1, nfolds = 10, type.measure = "mse")
family = "gaussian",
alpha = 0, # Perfil do MSE de validação em função de lambda.
parallel = TRUE, par(mfrow = c(2, 1))
type.measure = "mse") plot(cvfit_ridge, main = "RIDGE")
abline(v = log(cvfit_ridge$lambda.min), col = 2)
plot(cv.ridge) plot(cvfit_lasso, main = "LASSO")
cv.ridge$lambda.min abline(v = log(cvfit_lasso$lambda.min), col = 2)
cv.ridge$lambda.1se layout(1)
coef(cv.ridge, s = cv.ridge$lambda.min)
# Parâmetros de penalização ótimos.
cv.lasso <- cv.glmnet(X, c(Ridge = cvfit_ridge$lambda.min,
y, LASSO = cvfit_lasso$lambda.min)
family = "gaussian",
alpha = 1, # Coeficientes estimados (variáveis selecionadas).
parallel = TRUE, cbind(round(coef(cvfit_ridge, s = "lambda.min"), digits = 4),
type.measure = "mse") round(coef(cvfit_lasso, s = "lambda.min"), digits = 4))
plot(cv.lasso)
cv.lasso$lambda.min
cv.lasso$lambda.1se
coef(cv.lasso, s = cv.lasso$lambda.min)
par(mfrow = c(1, 2)) par(mfrow = c(1, 2))
plot(mridge, xvar = "lambda", main = "RIDGE") plot(mridge, xvar = "lambda", main = "RIDGE")
abline(v = log(cv.ridge$lambda.min), lty = 2) abline(v = log(cvfit_ridge$lambda.min), lty = 2)
plot(mlasso, xvar = "lambda", main = "LASSO") plot(mlasso, xvar = "lambda", main = "LASSO")
abline(v = log(cv.lasso$lambda.min), lty = 2) abline(v = log(cvfit_lasso$lambda.min), lty = 2)
layout(1) layout(1)
``` ```
......
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