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

Adiciona material de validação cruzada.

parent 79dfd104
---
title: "Validação Cruzada"
author: Prof. Walmes M. Zeviani & Prof. Eduardo V. Ferreira
date: 2017-08-24
#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)
```
# Definindo o modelo
```{r, messages = FALSE}
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(lattice)
library(latticeExtra)
```
```{r}
#-----------------------------------------------------------------------
# Gerando dados.
curve(10 + 0.05 * x + 2 * sin(x/10), 0, 100)
# 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.
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 = "orange", lwd = 3)
```
# Validação cruzada *holdout*
```{r}
#-----------------------------------------------------------------------
# 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 = "orange", lwd = 3)) +
layer(panel.points(x = x,
y = y,
pch = 19,
col = "red"),
data = das[["val"]])
# Calculação a validação cruzada do handout.
deg <- 1:dmax
cv <- sapply(deg,
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(ErrTra = Etra, ErrVal = Eval)
})
cv <- cbind(deg, as.data.frame(t(cv)))
# Capacidade de predição contra o grau do polinômio usado.
xyplot(ErrTra + ErrVal ~ deg,
data = cv,
type = "o",
auto.key = list(corner = c(0.95, 0.95)))
```
# Validação cruzada com $k$-*folds*
```{r}
#-----------------------------------------------------------------------
# Método k-fold.
# Gerando um conjunto maior de observações do mesmo modelo.
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 = "orange", 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)
# Avaliando cenários varrendo os fold e variando o grau do polinômio.
cen <- expand.grid(fold = 1:k,
deg = 1:15)
# Avaliando cada cenário.
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(ErrTra = ErrTra, ErrVal = ErrVal))
})
kfol <- cbind(cen, as.data.frame(t(kfol)))
# Curva de erro em cada fold.
xyplot(ErrTra + ErrVal ~ deg | fold,
data = kfol,
as.table = TRUE,
auto.key = list(corner = c(0.95, 0.95)),
type = "o")
# Os folds juntos para um mesmo tipo de erro.
xyplot(ErrTra + ErrVal ~ deg,
groups = fold,
data = kfol,
type = "o") +
layer(panel.xyplot(x = x,
y = y,
type = "a",
col = "black",
lwd = 2))
```
# Validação cruzada com *leave-one-out*
```{r}
#-----------------------------------------------------------------------
# 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(ErrTra = ErrTra, ErrVal = ErrVal))
})
nfol <- cbind(cen, as.data.frame(t(nfol)))
names(nfol)[4] <- "ErrVal"
xyplot(ErrTra + ErrVal ~ deg,
groups = fold,
data = nfol,
scales = list(y = "free"),
type = "o") +
layer(panel.xyplot(x = x,
y = y,
type = "a",
col = "black",
lwd = 2))
# Leave-one-out vs k-fold.
xyplot(ErrVal ~ deg,
data = nfol,
type = c("p", "a")) +
as.layer(xyplot(ErrVal ~ (deg + 0.15),
data = kfol,
col = "red",
type = c("p", "a")))
```
Quando usar o leave-one-out?
# Decomposição do erro quadrático médio de predição
Vamos fazer a decomposição do erro quadrático médio de predição de uma
nova observação $Y$ cujo vetor de covariáveis é $x_0$. O verdadeiro
modelo para a média de $Y$ é representado pela função $\eta(x)$,
conforme mostra a igualdade abaixo,
$$
Y_{x_0} = \eta(x_0) + e,
$$
em que $e$ é um termo aleatório de média 0 e variância constante
$\sigma^2$.
Então o erro quadrático médio é função da variância do erro, da
variância das predições usando o modelo estimado denotado por
$\hat{h}(x)$ e do vício, que é função da discrepância entre $\eta(x_0)$
e a média de $\hat{h}(x_0)$.
$$
\begin{align*}
\mathbb{E}[(Y_{x_0} - \hat{h}(x_0))^2]
&= \mathbb{E}[Y^2 + \hat{h}^2 - 2Y\hat{h}]\\
&= \underbrace{\mathbb{E}[Y^2]}_{
\mathbb{V}[Y] = \mathbb{E}[Y^2] - \mathbb{E}^2[Y]}
+ \mathbb{E}[\hat{h}^2] - \mathbb{E}[2Y\hat{h}]\\
&= \mathbb{V}[Y] + \mathbb{E}^2[Y]
+ \mathbb{V}[\hat{h}] + \mathbb{E}^2[\hat{h}]
- 2\mathbb{E}[Y\hat{h}] \\
&= \mathbb{V}[Y] + \mathbb{V}[\hat{h}]
+ (h^2 - 2h\mathbb{E}[\hat{h}] + \mathbb{E}^2[\hat{h}])\\
&= \mathbb{V}[Y] + \mathbb{V}[\hat{h}] + (h - \mathbb{E}[\hat{h}])^2\\
&= \sigma^2 + \mathbb{V}[\hat{h}] + \mathbb{B}^2[\hat{h}]
\end{align*}
$$
```{r}
#-----------------------------------------------------------------------
# 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))
# Quadrado do erro de predição.
(y - h)^2
#-----------------------------------------------------------------------
# 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
# deg <- c(1, 7)
# rpt <- 5000
# 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))
resul
xyplot(eqm + Vy + Vh + Bh ~ deg,
data = resul,
auto.key = list(corner = c(0.95, 0.95))) +
glayer(panel.smoother(..., se = FALSE))
```
```{r, eval = FALSE, include = FALSE}
#-----------------------------------------------------------------------
# Verificando a decomposição do EQM para estimação de parâmetros.
theta <- c(A = 10, B = 2)
fx <- function(x) {
theta["A"] * x/(theta["B"] + x) + rnorm(length(x), 0, 0.1)
}
x <- seq(2, 30, by = 2)
plot(y ~ x)
s <- 1000
r <- replicate(s, {
y <- fx(x)
n0 <- nls(y ~ A * x/(B + x), start = theta)
coef(n0)
})
e <- rowMeans(sweep(r,
MARGIN = 1,
STATS = rbind(10, 2),
FUN = function(x, y) {
(x - y)^2
}))
b <- (rowMeans(r) - theta)^2
v <- apply(r, MARGIN = 1, FUN = var) * (s - 1)/s
# v <- apply(r, MARGIN = 1, FUN = var)
e - (b + v)
```
```{r}
#-----------------------------------------------------------------------
# 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")
})
```
# *Leave-one-out*, DFfits e DFbetas
Para um material mais completo, visite:
<http://leg.ufpr.br/~walmes/cursoR/fiocruz/fiocruz04diagn.html>.
```{r}
# Dados.
plot(dist ~ speed, data = cars)
# Ajuste do modelo.
m0 <- lm(dist ~ speed + I(speed^2), data = cars)
# Índice das observações.
i <- 1:nrow(cars)
# Artefatos do modelo ajustado.
X <- model.matrix(m0)
iXlX <- solve(t(X) %*% X)
siXlX <- sqrt(diag(iXlX))
Hi <- sqrt(diag(X %*% iXlX %*% t(X)))
fit0 <- fitted(m0)
beta0 <- coef(m0)
loo <- lapply(i,
FUN = function(ii) {
m1 <- update(m0, data = cars[-ii, ])
fit <- predict(m1,
newdata = cars[ii, ])
s <- sqrt(deviance(m1)/df.residual(m1))
se_beta <- s * siXlX
se_fit <- s * Hi[ii]
list(dfbetas = (beta0 - coef(m1))/se_beta,
dffits = (fit0[ii] - fit)/se_fit)
})
head(dfbetas(m0))
head(t(sapply(loo, "[[", "dfbetas")))
head(dffits(m0))
head(sapply(loo, "[[", "dffits"))
```
html_document:
highlight: haddock
fig_caption: yes
fig_align: center
# code_folding: hide
toc: true
toc_dep: 3
# toc_float:
# collapsed: false
# smooth_scroll: false
number_sections: true
css: ../config/vignette-style.css
includes:
in_header: ../config/MathJax.html
before_body: ../config/vignette-header.html
after_body: ../config/vignette-footer.html
mathjax: default
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