Métodos de reamostragem
Laboratório de Estatística e Geoinformação - LEG/UFPR
Métodos de reamostragem são ferramentas indispensáveis na estatística moderna. As técnicas envolvem particionar os dados de treino e reajustar os modelos em competição para cada subamostra, a fim de obter informações adicionais sobre o ajuste do modelo (que não seria possível com os dados completos). P. ex., podemos estimar o erro de teste associado a determinado modelo (avaliação do modelo), selecionar o nível de flexibilidade adequado (seleção do modelo) etc. Na prática, treinamos o algoritmo com a amostra de treinamento (curva em azul da figura abaixo) e avaliamos a qualidade do treinamento com a amostra de validação (baseando-se na curva em vermelho). Note que se utilizarmos algoritmos muito simples (pouco flexíveis), teremos um erro de predição alto na amostra de treino (curva em azul), e a medida que aumentamos sua complexidade, o erro de treinamento tende a diminuir. Entretanto, essa melhora aparente vem acompanhada da redução na capacidade de generalizar a informação, ou seja, com a chegada de novos exemplos (amostra de validação), o desempenho não é satisfatório. Frente a essa questão, o desafio é encontrar um meio-termo entre o modelo simples, que aprendeu pouco (viciado), e o complexo que aprendeu demais (alta variabilidade), tal que com a chegada de novos dados, se erre o mínimo possível.
Função Custo
Antes de proceder para a escolha do melhor modelo, devemos definir métricas para dizer o quão bom eles são. Neste material, vamos apresentar algumas (muito embora essas não sejam as únicas opções). Para maiores aprofundamentos, sugere-se por exemplo, A Probabilistic interpretation of precision, recall and F-score.
- Matriz de confusão: Em problemas de classificação, uma matriz de confusão é um layout de tabela que permite a visualização do desempenho do algoritmo. Nas linhas são representadas as instâncias verdadeiras, enquanto nas colunas as preditas (ou vice-versa). Abaixo, temos dois exemplos, a esquerda refere-se ao diagnóstico de pacientes doentes e saudáveis, já à direita trata da questão de classificação de e-mails como spam e não spam.
- Acurácia: trata-se da proporção de resultados corretos que o classificador alcançou. Ou seja, a razão entre as predições corretas pelo total. No exemplo do diagnóstico médico, seria o número de classificações corretas, dentre todos os pacientes:
Já para o exemplo da classificação de e-mails, temos uma acurácia de 80%.
\[ \textrm{Acurácia} = \frac{1000+8000}{10000} = 90\% \]
Agora, imagine um problema de transações de cartão de crédito em que 248.335 das operações, apenas 472 são fraudulentas. Um modelo sugerido, segundo essa métrica, poderia ser aquele que atribuísse todas as transações como legais, acertando mais de 99% das vezes. Entretanto, ele não captaria as fraudes. Por isso, quando temos um grande desequilíbrio de classes, essa medida não é muito adequada, pois o modelo, geralmente, fará previsões segundo a classe majoritária, obtendo bons resultados, mas não será útil no contexto do problema. É o que chamamos de Paradoxo da acurácia.
\[ \textrm{Acurácia} = \frac{100+700}{1000} = 80\% \]
- Precisão: é o número de verdadeiros positivos dividido pelo número de positivos estimados pelo modelo (verdadeiro positivo + falso positivo). É também chamado de Positive Predictive Value (PPV). Para o exemplo do diagnóstico de pacientes, temos:
Uma Precisão de 55,6% pode parecer baixa, entretanto, neste caso, consideramos mais grave o fato de diagnosticar um paciente doente como saudável (falso negativo), do que o contrário. Sendo assim, esse valor pode não ser tão ruim quanto parece. Por outro lado, veja o exemplo da classificação de e-mails:
\[ \textrm{Precisão} = \frac{1000}{1000 + 800} = 55,6\% \]
Agora, o erro que consideramos mais danoso é atribuir acidentalmente um e-mail não spam para a caixa de spam, o que poderia trazer mais problemas ao usuário. Neste caso, a Precisão foi de 76,9%, o que pode ser considerado um valor alto.
\[ \textrm{Precisão} = \frac{100}{100 + 30} = 76,9\% \]
- Recall: dado que o estado verdadeiro é positivo, qual a proporção de verdadeiro positivo. No exemplo do diagnóstico, seria calcular a proporção dos diagnosticados como doente, dado que estamos na população dos pacientes doentes:
Por outro lado, o exemplo da classificação de e-mails, temos um Recall de 37%. Não é demais enfatizar que, neste caso, consideramos ser mais prejudicial classificar um não-spam como spam (falso positivo).
\[ \textrm{Recall} = \frac{1000}{1000 + 200} = 83,3\% \]
\[ \textrm{Recall} = \frac{100}{100 + 170} = 37\% \]
- \(F_1\) score: vimos nos itens anteriores que alguns modelos podem apresentar boa Precisão, mas nem sempre um Recall elevado, e o contrário também pode acontecer. Frente a esse fato, gostaríamos de uma medida que fizesse um compromisso entre essas duas medidas, atribuindo, inclusive, o mesmo grau de importância para ambas. A função \(F_1\) score faz isso, ela considera a média harmônica entre as duas entidades supracitadas, ou seja:
O valor do \(F_{1}\) score pode variar entre 0 e 1, e quanto mais próximo de 1 melhor o modelo. Note que se \(x<y\), então \(2\cdot \frac{x\cdot y}{x + y}< \frac{x + y}{2}\), ou seja, estamos sempre mais perto do menor valor. Esse fato faz com que a média harmônica seja mais robusta a valores extremos. Isto é, ao se observar uma quantidade baixa de Precisão ou Recall já se levanta uma bandeira de atenção, por exemplo:
\[ F_1\ score = 2\times \frac{\textrm{Precisão}\cdot \textrm{Recall}}{\textrm{Precisão} + \textrm{Recall}} \]
- Se Precisão = 1 e Recall = 0 (média aritmética = 0,5): \[ F_1\ score = 2\times \frac{1\cdot 0}{1 +0} = 0; \]
- Se Precisão = 0,2 e Recall = 0,8 (média aritmética = 0,5): \[ F_1\ score = 2\times \frac{0,2\cdot 0,8}{0,2 + 0,8} = 0,32; \]
- \(F_\beta\) score: como dito anteriormente, o \(F_1\) score atribui pesos iguais para Precisão ou Recall (o número “1” no \(F_{1}\) indica exatamente isso). No entanto, considerar o mesmo peso para tais quantidades, muitas vezes, não é adequado. O \(F_\beta\) score é uma generalização desse compromisso, em que o \(\beta\) representa a influência da Precisão no resultado final. Por exemplo, \(F_{0,5}\) significa que estamos considerado mais peso para Recall. Importante constatar que a escolha do \(\beta\) não é algo exato, e exige bastante intuição sobre o problema, além de algumas rodadas de experimento. A função generalizada tem a seguinte forma: \[ F_\beta\ score = (1 +\beta^2)\frac{\textrm{Precisão} \cdot \textrm{Recall}}{\beta^2\cdot \textrm{Precisão} + \textrm{Recall}} \]
- Erro quadrático médio (EQM): agora, apresentaremos algumas funções perda para o caso de regressão. A primeira, trata da soma do quadrado da diferença entre o valor observado, \(y_{obs}\), e o estimado, \(h(\mathbf{x})\). Essa métrica é mais útil quando erros grandes são particularmente indesejáveis (advindos de outliers, por exemplo), uma vez que os penaliza com mais intensidade, quando comparado com o erro absoluto médio (que veremos em seguida), sendo este último mais robusto (menos influenciado) a estas questões de valores atípicos. Definimos o EQM como: \[ J\left[y_i,h(\mathbf{x})\right] = \frac{1}{n}\sum_{i=1}^{n}\left[y_i - h(x_i)\right]^2 \] Quando se quer uma métrica na mesma unidade da variável de interesse, toma-se a raiz quadrada do erro quadrático médio (REQM). É importante constatar que não estamos tratando da medida utilizada para construir (estimar) \(h(\mathbf{x})\).
- Erro absoluto médio (EAM): mede a diferença média entre o valor observado, \(y_{obs}\), e o estimado, \(h(\mathbf{x})\), sem considerar a sua direção. Neste caso, todas as diferenças individuais têm o mesmo peso. Definimos o EAM como: \[ J\left[y_i,h(\mathbf{x})\right] = \frac{1}{n}\sum_{i=1}^{n}|y_i - h(x_i)| \] Se os sinais dos erros não forem removidos, estaremos diante do erro médio do vício, que também nos fornece informações úteis, mas deve ser interpretado com cautela, porque os erros positivos e negativos são cancelados.
- Huber-M cost: é um compromisso entre as funções 7 e 8, em que o parâmetro \(\delta\) é obtido automaticamente para um específico percentil dos erros absolutos: \[ J[y_i,h(\mathbf{x})] = \frac{1}{n}\displaystyle \sum_{i=1}^{n}\begin{cases} \frac{1}{2}[y_i - h(x_i)]^2, & \textrm{para } |y - h(\mathbf{x}_i)| \le \delta, \\ \delta\, |y_i - h(x_i)| - \frac{1}{2}\delta^2, & \textrm{caso contrário.} \end{cases} \]
Observação 1: \([EAM] ≤ [REQM]\), i.e., o resultado da raiz quadrada do erro quadrático médio será sempre maior ou igual ao erro absoluto médio. Se todos os erros tiverem a mesma magnitude, então \([EAM] = [REQM]\).
Observação 2: \([REQM] ≤ [ EAM \times \sqrt{n}]\), em que \(n\) é o número de observações de validação. Isso significa que a raiz quadrada do erro quadrático médio tem uma tendência a ser cada vez maior do que o erro absoluto médio, à medida que o tamanho da amostra de validação aumenta. Isso pode ser problemático quando compararmos resultados da REQM calculados em amostras de validação de tamanhos diferentes, o que é comum no mundo real.
Observação 3: uma grande vantagem do erro quadrático médio é evitar a função módulo, que é bastante indesejável em muitos cálculos matemáticos.
Validação cruzada holdout
Definida a função custo, abordaremos as técnicas de validação propriamente, a começar com a holdout. A atividade envolve dividir aleatoriamente o conjunto de dados em (i) treinamento: utilizado para preparar o modelo (treiná-lo) e (ii) validação: utilizado para avaliar o desempenho do modelo treinado (neste ponto é que consideramos as funções acima descritas). É importante constatar que algumas literaturas não fazem distinção entre validação e teste, muito embora sejam entidades diferentes. Os dados de validação são rotulados, assim como os de treinamento, mas são utilizados para testar o desempenho do modelo. Neste texto, consideramos os dados de teste como informações novas, ainda não rotuladas (veja a figura abaixo). Para entender esse processo, considere a função custo como a perda quadrática:
\[ J\left[y_i,h(\mathbf{x})\right] = \frac{1}{n}\sum_{i=1}^{n}[y_i - h(x_i)]^2. \] Neste caso, \(n\) representa os dados de validação. Já os dados de treinamento foram utilizados na estimação dos parâmetros de \(h(\mathbf{x})\). Não é demais enfatizar que essa partição deve ser disjunta, ou seja, durante a estimação do modelo, não se pode, em momento algum, apresentar observações que serão utilizadas na validação.
Vamos ilustrar os conteúdos apresentados neste capítulo através de um mesmo exemplo simulado. A função geradora dos dados é a seguinte:
\[\begin{eqnarray} f(x) &=& \beta_0 + \beta_1 * x + \beta_2 * sin\left(\frac{x}{10}\right)\nonumber\\[6pt] &=& 10 +0,05* x + 2 * sin\left(\frac{x}{10}\right).\nonumber \end{eqnarray}\]A linha em azul reprenta a curva verdadeira do sistema.
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(lattice)
library(latticeExtra)
# Retorna os valores determinados pelo modelo.
eta <- function(x, beta = c(10, 0.05, 2)) {
beta[1] + beta[2] * x + beta[3] * sin(x/10)
}
# Retorna um rúido normal padrão.
epsilon <- function(x) {
rnorm(n = length(x), mean = 0, sd = 1)
}
# Simulando do modelo.
set.seed(1)
da <- data.frame(x = runif(100, 0, 100))
da$y <- eta(da$x) + epsilon(da$x)
# Exibindo os dados.
plot(y ~ x, data = da)
curve(eta, add = TRUE, col = "dodgerblue3", lwd = 3)
Na figura abaixo, separamos aleatoriamente 100 observações em duas amostras, treinamento (com 80 dados) e validação (20 dados). Em cada painel, temos um polinômio ajustado em cinza, com os graus variando de 1 à 15. A curva verdadeira está representada em azul, e os pontos em vermelho são as observações de validação. Como esperado, a medida que aumentamos a flexibilidade do modelo, observa-se um melhor ajuste aos dados de treinamento.
#-----------------------------------------------------------------------
# Partindo os dados em treino e validação.
# Criando a variável para separar os dados em treino e validação.
ntra <- 80
nval <- 20
i <- sample(rep(c("tra", "val"),
times = c(ntra, nval)))
# Dividindo os dados em treino e validação.
das <- split(da, f = i)
str(das)
# Criando uma lista de fórmulas para repetir os dados nos paineis.
dmax <- 15
form <- replicate(dmax, y ~ x)
names(form) <- 1:dmax
# Diagramas de dispersão dos dados um ajuste variando o grau.
xyplot.list(form,
data = das[["tra"]],
as.table = TRUE) +
layer(panel.smoother(form = y ~ poly(x, degree = panel.number()),
method = lm)) +
layer(panel.curve(eta, col = "dodgerblue3", lwd = 3)) +
layer(panel.points(x = x,
y = y,
pch = 19,
col = "red"),
data = das[["val"]])
Entretanto, o nosso interesse está em treinar o modelo e avaliar a sua capacidade de generalização, isto é, quantificar seu desempenho com os dados de validação. Na figura que segue, apresentamos a raiz quadrada do erro quadrático médio (REQM) para os dados de validação e treinamento, de acordo com a flexibilidade do modelo. A escolha do grau do polinômio se dará pelo menor valor da REQM de validação (curva em azul). Nesse caso, o polinômio escolhido é o de grau 4.
# Calculação a validação cruzada do handout.
Flexibilidade <- 1:dmax
cv <- sapply(Flexibilidade,
FUN = function(d) {
# Ajuste com os dados de treino.
m0 <- lm(y ~ poly(x, degree = d),
data = das[["tra"]])
# Predição para os dados de validação.
yhat <- predict(m0,
newdata = das[["val"]])
# Erro de treino.
# crossprod(fitted(m0) - das[["tra"]]$y)/ntra
Etra <- deviance(m0)/ntra
# Erro de validação.
Eval <- crossprod(yhat - das[["val"]]$y)/nval
c(Treino = Etra, Validação = Eval)
})
cv <- cbind(Flexibilidade, as.data.frame(t(cv)))
# Capacidade de predição contra o grau do polinômio usado.
xyplot(Treino + Validação ~ Flexibilidade,
data = cv,
type = "o",
ylab="REQM",
auto.key = list(corner = c(0.95, 0.95)))
A validação por holdout, embora seja simples de implementar e barata computacionalmente, não é muito utilizada na prática. Com uma partição apenas, não se tem uma avaliação precisa sobre a qualidade do algoritmo (não se pode ficar dependente de uma boa partição). Ademais, em algumas situações, não se tem um número suficiente de dados para dividir em dois grupos satisfatoriamente, gerando modelos mal ajustados e aproximações para o erro do teste superestimadas.
Menos dados \(\Rightarrow\) geralmente menos informação \(\Rightarrow\) maior variabilidade
Validação cruzada K-fold
Como dito anteriormente, separar os dados em somente duas partes disjuntas pode trazer resultados divergentes, dependendo da informação contida em cada conjunto (especialmente quando os dados são escassos). A abordagem de validação cruzada por k-fold minimiza esses problemas. O método consiste em dividir os dados em \(K\) partes iguais, ajusta-se o modelo utilizando \(K-1\) partes, e a parcela restante fica destinada à validação. Esse processo é repetido \(K\) vezes (em cada momento uma partição diferente será a validação), em seguida os resultados são combinados obtendo a médias dos erros obtidos.
Definição
Sejam \(K\) partes denotadas por \(C_1,C_2,\dots, C_K\), em que \(C_k\) representa o índice da \(k\)-ésima parte.
Considere ainda que temos \(n_k\) observações na partição \(k\) (se \(n\) é múltiplo de \(K\), então \(n_k = \frac{n}{K}\)).
Calcule: \[ CV_{(K)} = \sum_{k=1}^{K} \frac{n_k}{n}EQM_k, \] em que \(EQM_k = \frac{\sum_{i \in C_k} (y_i - \hat y_i)^2}{n_k}\), e \(\hat y_i\) é o valor ajustado da observação \(i\), obtido dos dados com a \(k\)-ésima parte removida.
No exemplo abaixo, temos 200 observações simuladas. Utilizamos o método 10-fold, e em cada painel apresentamos uma partição dos dados. A curva verdadeira está representada em azul e os pontos em vermelho são as observações de validação. Em cada \(fold\), o modelo foi ajustado com as observações em preto (obtendo o modelo estimado), e avaliadas com as vermelhas (calculando o \(EQM_k\)). Assim, a média ponderada dos 10 \(EQM's\) serve de estimativa do erro de teste (dados ainda não rotulados), além de critério para escolha do modelo.
#-----------------------------------------------------------------------
# Método k-fold.
# Gerando um conjunto maior de observações do mesmo modelo.
set.seed(1)
da <- data.frame(x = runif(200, 0, 100))
da$y <- eta(da$x) + epsilon(da$x)
# nrow(da)
# Número de observações e quantidade de folds para dividir os dados.
n <- nrow(da)
k <- 10
da$i <- ceiling((1:n)/(n/k))
# Parte os dados em k conjuntos disjuntos.
das <- split(da, f = da$i)
# str(das)
# Replica a fórmula para exibir os dados.
form <- replicate(k, y ~ x)
names(form) <- sprintf("fold %d", 1:k)
xyplot.list(form,
data = da,
type = "n",
as.table = TRUE) +
# Curva do modelo verdadeiro.
layer(panel.curve(eta, col = "dodgerblue3", lwd = 3)) +
# Observações do conjunto de treinamento.
layer(panel.points(x = x[da$i != packet.number()],
y = y[da$i != packet.number()])) +
# Observações do conjunto de avaliação.
layer(panel.points(x = x[da$i == packet.number()],
y = y[da$i == packet.number()],
pch = 19,
col = "red")) +
# Ajuste do modelo aos dados de treinamento.
layer(panel.smoother(x = x[da$i != packet.number()],
y = y[da$i != packet.number()],
form = y ~ poly(x, degree = 5),
method = lm)) +
# Linhas de referência.
layer(panel.grid(), under = TRUE)
Seguindo com o exemplo, apresentamos o resultado da validação segundo a \(REQM\). Na amostra de treinamento (gráfico à esquerda), as curvas são bastante similares, o que era de se esperar, uma vez que trata-se de 90% dos dados, e os 10 folds gerados não mudam muito entre si. Por outro lado, note que para os dados de validação (à direita), dependendo do fold, temos diferentes valores para o \(REQM\) (a linha em preto representa a média geral). Entretanto, independente da partição utilizada na validação, a decisão sobre a flexibilidade do modelo não se altera (neste caso, polinômio de grau 4).
cen <- expand.grid(fold = 1:k,
deg = 1:15)
kfol <- mapply(f = cen$fold,
d = cen$deg,
FUN = function(f, d) {
# Ajuste do modelo aos dados de treinamento.
j <- da$i != f
m0 <- lm(y ~ poly(x, degree = d),
data = subset(da, j))
# Erro de treinamento.
ErrTra <- deviance(m0)/sum(j)
# Erro de validação.
yhat <- predict(m0, newdata = das[[f]])
ErrVal <- crossprod(yhat - das[[f]]$y)/length(yhat)
return(c(Treino = ErrTra, Validação = ErrVal))
})
kfol <- cbind(cen, as.data.frame(t(kfol)))
# Os folds juntos para um mesmo tipo de erro.
xyplot(Treino + Validação ~ deg,
groups = fold,
data = kfol,
ylab="REQM",
xlab="Flexibilidade",
type = "o") +
layer(panel.xyplot(x = x,
y = y,
type = "a",
col = "black",
lwd = 2))
Validação cruzada Leave-One-Out
Validação cruzada Leave-One-Out (LOOCV) é muito parecido com validação cruzada holdout. Ambos separam os dados em duas partes, mas o LOOCV utiliza somente uma observação na fase de validação. I.e., treina-se o modelo com \(n-1\) dados, \(\left\{(x_1,y_1), \dots, (x_{k-1},y_{k-1}),(x_{k+1},y_{k+1}),\dots,(x_n,y_n)\right\}\) e o avalia utilizando a observação restante, \((x_k,y_k)\). Repetimos esse procedimento \(n\) vezes, excluindo em cada momento uma observação diferente.
A validação cruzada via LOOCV apresenta vantagens sobre holdout, pois possui menor vício (não superestima o erro do teste), além de não ser tão influenciada pela aleatorização das partições. Todavia, como veremos na Seção Bias-Variance Trade-Off, o vício não é a única fonte de preocupação, devemos considerar ainda a variância (veja essa variabilidade na figura abaixo).
#-----------------------------------------------------------------------
# k-fold ou leave-one-out (n-fold)
# Cenários.
cen <- expand.grid(fold = 1:n,
deg = 1:15)
# Obtendo os erros para cada cenário.
nfol <- mapply(f = cen$fold,
d = cen$deg,
FUN = function(f, d) {
# Ajuste do modelo aos dados de treinamento.
m0 <- lm(y ~ poly(x, degree = d),
data = da[-f, ])
# Erro de treinamento.
ErrTra <- deviance(m0)/(n - 1)
# Erro de validação.
yhat <- predict(m0, newdata = da[f, ])
ErrVal <- (yhat - da$y[f])^2
return(c(Treino = ErrTra, Validação = ErrVal))
})
nfol <- cbind(cen, as.data.frame(t(nfol)))
names(nfol)[4] <- "Validação"
xyplot(Treino + Validação ~ deg,
groups = fold,
data = nfol,
ylab="REQM",
xlab="Flexibilidade",
scales = list(y = "free"),
type = "o") +
layer(panel.xyplot(x = x,
y = y,
type = "a",
col = "black",
lwd = 2))
Neste caso, LOOCV é pior avaliado quando comparado com outros \(k-fold\), com \(k<n\). Isso ocorre devido a correlação (positiva) entre as estimativas. Quanto menos dobras (partições), menos overlap entre os dados de treinamento. Além disso, deve-se analisar os custos computacionais envolvidos em cada etapa, o LOOCV pode ser muito caro, dependendo do conjunto de dados. Na figura que segue, temos uma comparação entre os resultados via Leave-one-out (em preto) versus 10-fold (em vermelho). Muito embora o comportamento médio seja bastante semelhante, a variabilidade do LOOCV é superior, podendo trazer implicações na estimação do erro de predição.
xyplot(Validação ~ deg,
data = nfol,
ylab="REQM",
xlab="Flexibilidade",
type = c("p", "a")) +
as.layer(xyplot(Validação ~ (deg + 0.15),
data = kfol,
col = "red",
type = c("p", "a")))
Validação cruzada: certo e errado
Um fato importante que merece destaque é sobre separar os dados para validação depois de se realizar algum pré-processamento. Imagine que temos 800 variáveis, a após determinado filtro (por correlação das classes, por exemplo) chega-se em 200 variáveis, como na figura abaixo.
Não podemos realizar validação cruzada considerando apenas as 200 variáveis, pois o procedimento de reduzir o espaço das características gerou um aprendizado que não pode ser ignorado. O correto é separar os dados para validação antes de qualquer manipulação nos dados, como apresentado na figura que segue.
Qual valor de \(K\) escolher?
Diante de tantas possibilidades de validação, surge a dúvida sobre definir o valor de \(K\). Divisões simples (holdout) fazem um bom trabalho quando o número de hiperparâmetros é baixo e a quantidade de dados de treinamento e validação é elevada. Entretanto, o estimador do erro de teste é viciado. Já quando \(K=n\), tem-se um custo de processamento considerável, e apesar do estimador ser aproximadamente não viciado, ele apresenta alta variabilidade. Na prática, são recomendados \(K=5\) ou 10 (veja Kohavi, 1995), tendo assim um compromisso entre custo computacional, viés e variância.
Bias-Variance Trade-Off
Um bom desempenho do método em um conjunto de teste requer um baixo erro quadrático médio: \[ \mathrm{E} \left[y_0 - h(\mathbf{x}_0)\right]^2 = \mathrm{V}ar\left[h(\mathbf{x}_0)\right] + \left[\textrm{Vício}(h(\mathbf{x}_0)\right]^2 + \mathrm{V}ar(\varepsilon). \]
Em geral, quanto mais flexível o modelo, maior a sua variância. Essa componente se refere ao quanto \(h(\mathbf{x}_0)\) muda quando a estimamos utilizando diferentes dados de treino. Note que idealmente as estimativas de \(h(\mathbf{x}_0)\) não deveriam mudar muito entre os conjuntos.
Por outro lado, quanto mais simples o modelo, maior será seu vício. O viés representa o erro de aproximar um problema real (extremamente complicado) por uma função simples.
Para ilustrar essa questão, considere a situação abaixo. A linha laranja representa a curva real do modelo, e em cada painel temos curvas ajustadas de acordo com o grau de flexibilidade do polinômio (de 1 a 15). Note que no primeiro gráfico, temos um modelo ajustado bastante viesado, e a medida que elevamos o grau do polinômio, o vício diminui, mas a variabilidade aumenta substancialmente.
#-----------------------------------------------------------------------
# Decomposição em variância e vício.
# Valor considerado para avaliar a predição (x_0).
x0 <- 50
# Gerando dados artificiais.
da <- data.frame(x = runif(100, 0, 100))
da$y <- eta(da$x) + epsilon(da$x)
# Observação futura a ser prevista pelo modelo.
y <- eta(x0) + epsilon(x0)
# Ajuste do modelo.
m0 <- lm(y ~ x, data = da)
h <- predict(m0, newdata = list(x = x0))
#-----------------------------------------------------------------------
# Repetindo o processo várias vezes variando o grau do polinômio.
# Valor da covariável e estimativa de Y pelo modelo verdadeiro.
x0 <- 50
eta_x0 <- eta(x0)
# Graus do polinômio a serem avaliados e quantidade de repetições.
deg <- 1:15
rpt <- 100
# Aplica para cada grau.
resul <- lapply(deg,
FUN = function(d) {
t(replicate(rpt, {
da <- data.frame(x = sample(0:100, size = 30))
da$y <- eta(da$x) + epsilon(da$x)
m0 <- lm(y ~ poly(x, degree = d),
data = da)
y <- eta_x0 + epsilon(x0)
h <- predict(m0,
newdata = list(x = x0))
names(h) <- NULL
return(c(yobs = y,
hfit = h))
}))
})
# Calcula os termos da decomposição do EQM.
eqm_decom <- function(m) {
c(eqm = mean((m[, "yobs"] - m[, "hfit"])^2),
Vy = var(m[, "yobs"]) * (rpt - 1)/rpt,
Vh = var(m[, "hfit"]) * (rpt - 1)/rpt,
Bh = (eta_x0 - mean(m[, "hfit"]))^2)
}
resul <- t(sapply(resul, eqm_decom))
resul <- cbind(deg = deg, as.data.frame(resul))
#-----------------------------------------------------------------------
# Visualizando este conceito.
# Simulando dados.
da <- data.frame(x = sample(0:100, size = 30))
da$y <- eta(da$x) + epsilon(da$x)
# Cria lista com fórmulas.
dmax <- 15
form <- replicate(dmax, y ~ x)
names(form) <- 1:dmax
xyplot.list(form,
data = da,
type = "n",
as.table = TRUE,
panel = function(...) {
# Gráficos dos pontos originais (não exibidos).
panel.xyplot(...)
# Simulando a resposta e ajustando o modelo.
d <- panel.number()
xseq <- 0:100
resul <- replicate(30, {
xnew <- sample(0:100, size = 30)
ynew <- eta(xnew) + epsilon(xnew)
panel.points(x = xnew,
y = ynew,
pch = 20,
cex = 0.8,
alpha = 0.25)
m0 <- lm(ynew ~ poly(xnew, degree = d))
yfit <- predict(m0,
newdata = list(xnew = xseq))
panel.lines(x = xseq,
y = yfit,
col = "black",
alpha = 0.7)
return(yfit)
})
# Ajuste médio ponto a ponto.
yfitm <- rowMeans(resul)
panel.lines(x = xseq,
y = yfitm,
col = "green2",
lwd = 3)
# O verdadeiro modelo.
panel.curve(eta,
col = "orange",
lwd = 3)
panel.abline(v = x0, lty = 2, col = "gray")
})
Bootstrap
O Bootstrap é frequentemente utilizado para quantificar a incerteza associada a determinado estimador ou método de aprendizagem. Por exemplo, pode ser usado para estimar o erro padrão dos coeficientes do ajuste de uma regressão linear (ou seus intervalos de confiança). Na verdade, o poder dessa ferramenta reside no fato de ser aplicável a uma gama de métodos estatísticos, incluindo aqueles para os quais uma medida de variação é difícil de obter (e não é fornecido automaticamente pelo software).
Suponha que temos uma amostra de tamanho \(n\) de uma população. É esperado que essa amostra guarde as mesmas características da população. A ideia do bootstrap consiste em considerá-la como a população (“população bootstrap”), e realizar amostragem dessa “população”. Sendo assim, coletamos aleatoriamente (com reposição) \(B\) amostras de tamanho \(n\) (chamadas de amostras bootstrap). Para cada uma das \(B\) amostras, somos capazes de calcular qualquer estimativa de interesse, além do desvio padrão associado. A figura que segue resume esse processo.
Erro de predição e Bootstrap
Na validação cruzada, o \(k\)-ésimo \(fold\) de validação é distinto dos demais \(k-1\) \(folds\) usados no treinamento. Ou seja, não há overlap entre os dados de treino e validação. Esse fato é crucial para o sucesso da abordagem. Por outro lado, as amostras bootstrap apresentam um significativo overlap entre os dados (cerca de \(2/3\)), e a sua utilização para predizer o erro do teste pode trazer valores subestimados. Por esse motivo, consideramos a validação cruzada como uma abordagem mais simples e atrativa para estimar o erro de predição, e deixamos para utilizar bootstrap em outros casos que veremos durante o curso.