Métodos de Reamostragem
Autores: Alcides Conte Neto, Bruna Davies Wundervald, Marcelo Adriano Corrêa Maceno.
Data: 19 de Abril de 2016
## install.packages("ISLR")
library(ISLR)
## install.packages("ggplot2")
library(ggplot2)Para 2 grupos
Método de Validação
O método de validação é utilizado para estimar o erro do teste que resulta do ajuste de vários modelos para os dados observados. Esse método é um caso particular de Validação Cruzada quando K=2.
attach(Auto) # Dataset## The following object is masked from package:ggplot2:
##
## mpgcolor <- "#0000EE"Utilizou-se neste exemplo observações sobre carros (392 observações), sendo a variável explicativa o horsepower (potência do motor) e a variável resposta o mpg (milhas por galão).
plot(mpg ~ horsepower, col = color, xlab="Potência do motor",
ylab="Milhas por galão")
Criou-se um subconjunto de treinamento (train), sendo que este subconjunto contém metade das observações.
train <- sample(392, 196)Este subconjunto varia de acordo com os dados selecionados. Após isso, foi feito o ajuste do subconjunto de treino e calculado o Erro Quadrático Médio (EQM) utilizando o subconjunto de validação. O cálculo do EQM foi feito tanto para um modelo utilizando um polinômio de grau 1 quanto para um modelo utilizando grau 2.
lm.fit <- lm(mpg ~ horsepower, data = Auto, subset = train)
mean((mpg - predict(lm.fit, Auto))[-train]^2)## [1] 24.27107lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data = Auto, subset = train)
mean((mpg - predict(lm.fit2, Auto))[-train]^2)## [1] 20.01635Verifica-se que o EQM varia de acordo com o grau do polinômio ajustado. Assim, fez-se um ajuste para diferentes graus de polinômio.
EQMpoly <- rep(0,5)
for (i in 1:5) {
lm.fit <- lm(mpg ~ poly(horsepower, i), data = Auto, subset = train)
EQMpoly[i] <- mean((mpg - predict(lm.fit, Auto))[-train]^2)
}
plot(EQMpoly, type = "l", col = color, xlab = "Grau do Polinômio",
ylab = "EQM")
Neste caso, a partir do polinômio de grau 2 não se verifica uma diminuição do EQM. Assim, calcula-se o EQM para diferentes subconjuntos de treinamento, com um ajuste linear e quadrático.
x <- c(1:10)
EQM1 <- rep(0,10)
EQM2 <- rep(0,10)
for (i in 1:10) {
train <- sample(392,196)
lm.fit1 <- lm(mpg ~ poly(horsepower, 1), data = Auto, subset = train)
lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data = Auto, subset = train)
EQM1[i] <- mean((mpg - predict(lm.fit1, Auto))[-train]^2)
EQM2[i] <- mean((mpg - predict(lm.fit2, Auto))[-train]^2)
}
df <- data.frame(EQM1, EQM2)
ggplot(df, aes(x)) +
geom_line(aes(y=EQM1), colour="blue", linetype="dashed") +
geom_point(aes(y=EQM1), colour="blue", shape=15, size=2) +
geom_line(aes(y=EQM2), colour="red", linetype="dashed")+
geom_point(aes(y=EQM2), colour="red", shape=15, size=2) +
xlab("Subconjunto de treinamento") + ylab("EQM")
Verifica-se que o EQM varia com o subconjunto selecionado para treinamento.
Para K grupos
Validação Cruzada
As funções abaixo servem para obter k grupos independentes de um banco de dados, com o objetivo de separar a amostra de treino e a de validação.
lc <- function(n, k){
stopifnot(k > 0 && n > 0)
aux <- n
nk <- k
i <- vector("numeric", k)
while(nk > 0) {
i[nk] <- round(aux/nk)
aux <- aux-i[nk]
nk <- nk-1
}
return(i)
}
kFold <- function(dados, k = 2){
n <- dim(dados)[1]
interval <- lc(n, k)
res <- vector("raw", k)
for(i in 1:k){
temp <- sample(dim(dados)[1], interval[i],
replace = FALSE)
res[i] <- list(dados[temp, ])
dados <- dados[-temp, ]
}
names(res) <- 1:k
return(res)
}A primeira função somente obtém a quantidade de elementos \(n_k\) que cada intervalo deve ter:
lc(dim(Auto)[1], 5)## [1] 79 78 79 78 78Já a segunda função (kFold) é a que realiza o trabalho de separar os grupos:
k5 <- kFold(Auto, 5)
str(k5)Estimando o Erro de Previsão
Depois de efetuar a separação das classes, uma estimativa do erro de previsão do modelo validado em cada classe contra as restantes é dada pelo \(EQM_k\), que é:
\[ EQM_k = \frac{\sum\limits_{i \in C_k} (y_i - \hat{y}_i) ^ 2}{n_k} \]
Neste caso, como o \(EQM_k\) é calculado k vezes, foi contruida uma função para efetuar as respectivas permutações (entre Treino e Validação) e retornar o \(EQM_k\):
calcEQM <- function(dados, k, g){
dataSet <- kFold(dados, k)
eqm <- c()
eqmt <- c()
for(i in 1:k){
valida <- as.data.frame(dataSet[i])
names(valida) <- names(dados)
treino <- data.frame()
for(elem in dataSet[-i]){
treino <- rbind(treino, elem)
}
names(treino) <- names(dados)
mTreino <- lm(mpg ~ poly(horsepower, g), data = treino)
eqmt <- c(eqmt, mean(mTreino$residuals^2))
eqm <- c(eqm, mean((valida$mpg - predict(mTreino, valida))^2))
}
message("Média do EQM do conjunto de Treino: ", mean(eqmt))
return(eqm)
}É possível conferir abaixo o resultado da aplicação da função calcEQM:
eqm <- calcEQM(Auto, 5, 2)## Média do EQM do conjunto de Treino: 18.9382734313284names(eqm) <- 1:5
eqm## 1 2 3 4 5
## 18.65788 21.87342 17.79582 20.74119 18.02026Uma estimativa do \(EQM\) do modelo é ajustada pela fórmula:
\[ CV_{(k)} = \sum\limits_{k=1}^k \frac{n_k}{n} EQM_k \]
Essa equação faz uma ponderação do \(EQM_k\) pela quantia de observação em cada classe.
cv <- function(dados, eqm){
return(((1/dim(dados)[1]) * lc(dim(dados)[1], length(eqm))) %*% cbind(eqm))
}cv(Auto, eqm)## eqm
## [1,] 19.41164É fácil visualizar que quando a quantia de elementos em cada classe for igual, teremos apenas a média dos \(EQM_k\).
No exemplo acima, como a quantia de elementos em cada classe difere minimamente, o ajuste fica muito próximo da média dos \(EQM_k\):
mean(eqm)## [1] 19.41771Somente para comparação, abaixo segue os gráficos do uso de k-Fold e do Holdout:
for(j in 1:15){
EQMpoly <- 1:5
if (j <= 5){
for (i in 1:5) {
EQMpoly[i] <- cv(Auto, calcEQM(Auto, 2, i))
}
if(j == 1){
plot(EQMpoly, type = "l", col = color,
ylab = "EQM", xlab = "Grau do Polinômio", ylim = c(15, 30),
main = "K-fold")
} else {
lines(EQMpoly, type = "l", col = color)
}
} else if (j <= 10) {
for (i in 1:5) {
EQMpoly[i] <- cv(Auto, calcEQM(Auto, 5, i))
}
lines(EQMpoly, type = "l", col = 2)
} else if (j <= 15) {
for (i in 1:5) {
EQMpoly[i] <- cv(Auto, calcEQM(Auto, 10, i))
}
lines(EQMpoly, type = "l", col = 3)
}
}
legend(4, 30, legend = c(2, 5, 10), col = c(color, 2, 3), lty = 1)
for(j in 1:10){
train <- sample(392,196)
EQMpoly <- rep(0,5)
for (i in 1:5) {
lm.fit <- lm(mpg ~ poly(horsepower, i), data = Auto, subset = train)
EQMpoly[i] <- mean((mpg - predict(lm.fit, Auto))[-train]^2)
}
if(j == 1){
plot(EQMpoly, type = "l", col = color, xlab = "Grau do Polinômio",
ylab = "EQM", ylim = c(15, 30), main = "Holdout")
} else {
lines(EQMpoly, col = sample(colors(), 1))
}
}
Com esses gráficos, é possível perceber que ocorre uma variação maior no uso do Holdout.
Abaixo segue o gráfico do modelo usando o método Leave-One-Out:
EQMpoly <- rep(0,5)
for (i in 1:5) {
EQMpoly[i] <- cv(Auto, calcEQM(Auto, dim(Auto)[1], i))
}## Média do EQM do conjunto de Treino: 23.943297116855## Média do EQM do conjunto de Treino: 18.9844341472522## Média do EQM do conjunto de Treino: 18.9444978682445## Média do EQM do conjunto de Treino: 18.875652211491## Média do EQM do conjunto de Treino: 18.4262167119446plot(EQMpoly, type = "l", col = 2, main = "Leave-One-Out",
xlab = "Grau do Polinômio", ylab = "EQM", ylim = c(15, 30))
Bootstrap
O bootstrap é o método de gerar uma nova amostra do mesmo tamanho da original a partir do sorteio aleatório (c/ reposição) dos seus elementos. Ele pode ser usado para quantificar a incerteza associada a determinado estimador ou método de aprendizagem estatística - o que é o caso.
#install.packages("boot")
library(boot)
# install.packages("scatterplot")
library(car)Fazendo o bootstrap com R = 392 [tamanho do dataset].
- O argumento statistic é uma função que retorna um vetor contendo as estatísticas de interessse.
# install.packages("plyr")
library(plyr)## Warning: package 'plyr' was built under R version 3.2.5f.coef <- function(formula, data, indices){
d <- data[indices,]
fit <- lm(formula, data = d)
return(coef(fit))
}
fitb <- lm(mpg ~ horsepower, data = Auto)
ream.coef.n <- rlply(500, function(x){
lm.b <- boot(data = Auto, statistic = f.coef, R = 392, formula = (mpg ~ horsepower))
# Intervalos de Confiança para o LM
ci1 <- boot.ci(lm.b, type = "norm", index = 1)
ci2 <- boot.ci(lm.b, type = "norm", index = 2)
return(c(ci1$normal[2:3], ci2$normal[2:3]))
})
tab1 <- as.data.frame(do.call(rbind, ream.coef.n))
scatterplot(1:dim(tab1)[1], tab1$V1, col = color, xlab = "Número de Reamostragens BootStrap",
ylab = "Valores dos Coeficientes Obtidos - LI")
scatterplot(1:dim(tab1)[1], tab1$V2, col = color, xlab = "Número de Reamostragens BootStrap",
ylab = "Valores dos Coeficientes Obtidos - LS")
ream.coef.p <- rlply(500, function(x){
lm.b <- boot(data = Auto, statistic = f.coef, R = 392, formula = (mpg ~ horsepower))
# Intervalos de Confiança para o LM
ci1 <- boot.ci(lm.b, type = "perc", index = 1)
ci2 <- boot.ci(lm.b, type = "perc", index = 2)
return(c(ci1$percent[4:5], ci2$percent[4:5]))
})
tab2 <- as.data.frame(do.call(rbind, ream.coef.p))
scatterplot(1:dim(tab2)[1], tab2$V1, col = color, xlab = "Número de Reamostragens BootStrap",
ylab = "Valores dos Coeficientes Obtidos - LI")
scatterplot(1:dim(tab2)[1], tab2$V2, col = color, xlab = "Número de Reamostragens BootStrap",
ylab = "Valores dos Coeficientes Obtidos - LS")
Sabatina
- Pegar um banco de dados relativamente pequeno.
- Realizar o cálculo do erro pelo método de cross-validation k-Fold, com k escolhido.
- Realizar o cálculo de erro pelo método de cross-validation Leave-One-Out para polinômios de grau 1 até grau 5.
- Plotar o gráfico do EQM pelo grau do polinômio do k-Fold.
- Plotar o gráfico do EQM pelo grau do polinômio do Leave-One-Out.
- Fazer o intervalo de confiança bootstrap (percentil).
Referências
JAMES, G.; WITTEN, D.; HASTIE, T. e TIBSHIRANI. An Introduction to Statistical Learning. (livro-texto). Unofficial Solutions, 2013.
KOHAVI, R. A study of cross-validation and bootstrap for accuracy estimation and model selection. In: International joint Conference on artificial intelligence.. [S.l.: s.n.], 1995.