Machine Learning

leg.ufpr.br/~walmes/ensino/ML

1 Definindo o modelo

Para demonstrar o uso funcionamento da validação cruzada, será usado um conjunto de dados simulado apenas pela conveniência de serem conhecidos parâmetros.

As observações serão provenientes do modelo \[ y = \beta_{0} + \beta_{1} x + \beta_{3} \sin(x/\beta_{4}) + \epsilon, \]

em que \(\beta = [\beta_0, \beta_1, \beta_2, \beta_3]\) são os parâmetros de regressão e \(\epsilon\) um erro normal de média 0 e variância \(\sigma^2\). Os reais valores usados para esses parâmetros na simulação estão expostos dentro do código.

#-----------------------------------------------------------------------
# Pacotes.

rm(list = objects())
library(lattice)
library(latticeExtra)
#-----------------------------------------------------------------------
# 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 ruído 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)
legend("bottomright",
       bty = "n",
       legend = c("Sinal verdadeiro", "Valores observados"),
       lty = c(1, 0),
       lwd = c(3, 1),
       col = c("orange", "black"),
       pch = c(NA, 1))

2 Validação cruzada holdout

Na validação cruzada holdout é deixado de fora, para avaliação do desempenho de predição, uma porporção dos dados. Com a porção complementar, um modelo é ajustado aos dados, ocasião na qual são estimadas as quantidade livres do modelo, ou seja, os parâmetros.

#-----------------------------------------------------------------------
# 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)
## List of 2
##  $ tra:'data.frame': 80 obs. of  2 variables:
##   ..$ x: num [1:80] 75.9 80.3 71.6 24.5 63 ...
##   ..$ y: num [1:80] 15.4 17.8 15.1 12 14.5 ...
##  $ val:'data.frame': 20 obs. of  2 variables:
##   ..$ x: num [1:20] 40.5 89.9 88.8 36.6 95.7 ...
##   ..$ y: num [1:20] 10.7 15.4 14.6 10.7 16 ...
# Criando uma lista de fórmulas para repetir os dados nos paineis.
dmax <- 15
form <- replicate(dmax, y ~ x)
names(form) <- sprintf("pol. grau: %d", 1:dmax)

# Diagramas de dispersão dos dados um ajuste variando o grau.
xyplot.list(form,
            data = das[["tra"]],
            cex = 0.6,
            layout = c(3, 5),
            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,
                       cex = 0.6,
                       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 com os dados de treino e de teste.
                 yhat_inter <- predict(m0, newdata = das[["tra"]]) # Dados de treino.
                 yhat_exter <- predict(m0, newdata = das[["val"]]) # Dados de teste.
                 # Erros de predição
                 Etra <- crossprod(das[["tra"]]$y - yhat_inter)/ntra
                 Eval <- crossprod(das[["val"]]$y - yhat_exter)/ntra
                 # Retorno.
                 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)))

3 Validação cruzada com \(k\)-folds

Na validação cruzada do tipo \(k\)-folds, os dados são dividos nos registros em 5 conjuntos de tamanho aproximadamente igual (igual quando \(n\) for múltiplo de \(k\)).

Ao todo são \(k\) situações no qual se faz exatamente o que foi feito no procedimento holdout. Todavia, com essa divisão, cada registro ficará como observação no conjunto de teste pelo menos uma vez. Isso faz com que a avaliação do modelo tenha maior cobertura em termos de condições.

#-----------------------------------------------------------------------
# Método k-fold.

# Gerando um conjunto maior de observações do mesmo modelo.
da <- data.frame(x = runif(207, 0, 100))
da$y <- eta(da$x) + epsilon(da$x)
nrow(da)
## [1] 207
# 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))

# Número de registros por fold.
table(da$i)
## 
##  1  2  3  4  5  6  7  8  9 10 
## 20 21 21 20 21 21 20 21 21 21
# Parte os dados em k conjuntos disjuntos (gera uma lista).
das <- split(da, f = da$i)
str(das)
## List of 10
##  $ 1 :'data.frame':  20 obs. of  3 variables:
##   ..$ x: num [1:20] 34.7 84.2 61.4 55.9 68 ...
##   ..$ y: num [1:20] 11.1 14.9 11.7 11.2 12.2 ...
##   ..$ i: num [1:20] 1 1 1 1 1 1 1 1 1 1 ...
##  $ 2 :'data.frame':  21 obs. of  3 variables:
##   ..$ x: num [1:21] 12.5 32.5 100 70.5 52.2 ...
##   ..$ y: num [1:21] 13.6 11.1 12.4 14.1 10 ...
##   ..$ i: num [1:21] 2 2 2 2 2 2 2 2 2 2 ...
##  $ 3 :'data.frame':  21 obs. of  3 variables:
##   ..$ x: num [1:21] 25.4 26.2 63.7 51.8 29.5 ...
##   ..$ y: num [1:21] 11.56 13 14.27 10.57 9.39 ...
##   ..$ i: num [1:21] 3 3 3 3 3 3 3 3 3 3 ...
##  $ 4 :'data.frame':  20 obs. of  3 variables:
##   ..$ x: num [1:20] 24.731 88.563 0.686 0.377 91.049 ...
##   ..$ y: num [1:20] 13.84 14.69 9.56 9.31 14.06 ...
##   ..$ i: num [1:20] 4 4 4 4 4 4 4 4 4 4 ...
##  $ 5 :'data.frame':  21 obs. of  3 variables:
##   ..$ x: num [1:21] 85.4 60.6 69.7 87.2 87 ...
##   ..$ y: num [1:21] 15.1 12.5 15.2 15.8 15.3 ...
##   ..$ i: num [1:21] 5 5 5 5 5 5 5 5 5 5 ...
##  $ 6 :'data.frame':  21 obs. of  3 variables:
##   ..$ x: num [1:21] 94.04 28.34 51.76 99.39 4.23 ...
##   ..$ y: num [1:21] 14.68 9.66 10.02 13.39 11.15 ...
##   ..$ i: num [1:21] 6 6 6 6 6 6 6 6 6 6 ...
##  $ 7 :'data.frame':  20 obs. of  3 variables:
##   ..$ x: num [1:20] 75.9 74.2 70.9 30.4 76 ...
##   ..$ y: num [1:20] 16.1 14.5 14.8 11.5 16.8 ...
##   ..$ i: num [1:20] 7 7 7 7 7 7 7 7 7 7 ...
##  $ 8 :'data.frame':  21 obs. of  3 variables:
##   ..$ x: num [1:21] 76.3 58.1 41.5 62 34 ...
##   ..$ y: num [1:21] 16.4 11.7 10.2 13.9 11.9 ...
##   ..$ i: num [1:21] 8 8 8 8 8 8 8 8 8 8 ...
##  $ 9 :'data.frame':  21 obs. of  3 variables:
##   ..$ x: num [1:21] 76.62 3.76 29.88 86.4 56.36 ...
##   ..$ y: num [1:21] 16.6 10.4 12.5 15.5 12.2 ...
##   ..$ i: num [1:21] 9 9 9 9 9 9 9 9 9 9 ...
##  $ 10:'data.frame':  21 obs. of  3 variables:
##   ..$ x: num [1:21] 65.8 18.9 45.4 59.1 51.1 ...
##   ..$ y: num [1:21] 13.8 14.4 12.2 13 11.3 ...
##   ..$ i: num [1:21] 10 10 10 10 10 10 10 10 10 10 ...
# 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",
            xlab = "Grau do polinômio",
            ylab = "Erro",
            cex = 0.6,
            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,
                       cex = 0.6,
                       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))
                   # Predição com os dados de treino e de teste.
                   yhat_inter <- predict(m0, newdata = subset(da, j)) # Dados de treino.
                   yhat_exter <- predict(m0, newdata = das[[f]])      # Dados de teste.
                   # Erros de predição
                   Etra <- crossprod(subset(da, j)$y - yhat_inter)/length(yhat_inter)
                   Eval <- crossprod(das[[f]]$y - yhat_exter)/length(yhat_exter)
                   # Retorno.
                   c(ErrTra = Etra, ErrVal = Eval)
               })

kfol <- cbind(cen, as.data.frame(t(kfol)))

# Curva de erro em cada fold.
xyplot(ErrTra + ErrVal ~ deg | factor(fold),
       data = kfol,
       xlab = "Grau do polinômio",
       ylab = "Erro",
       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,
       xlab = "Grau do polinômio",
       ylab = "Erro",
       groups = fold,
       data = kfol,
       type = "o") +
    layer(panel.xyplot(x = x,
                       y = y,
                       type = "a",
                       col = "black",
                       lwd = 2))

4 Validação cruzada com leave-one-out

#-----------------------------------------------------------------------
# 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,
       xlab = "Grau do polinômio",
       ylab = "Erro",
       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?

5 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}] + (\eta^2 - 2h\mathbb{E}[\hat{h}] + \mathbb{E}^2[\hat{h}])\\ &= \mathbb{V}[Y] + \mathbb{V}[\hat{h}] + (\eta - \mathbb{E}[\hat{h}])^2\\ &= \sigma^2 + \mathbb{V}[\hat{h}] + \mathbb{B}^2[\hat{h}] \end{align*} \]

#-----------------------------------------------------------------------
# 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
##        1 
## 3.144327
#-----------------------------------------------------------------------
# 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
##    deg      eqm        Vy         Vh           Bh
## 1    1 6.103615 1.0511001 0.08176052 5.320882e+00
## 2    2 3.716172 0.8337139 0.16464649 2.493738e+00
## 3    3 4.351430 1.0059225 0.17915714 2.651481e+00
## 4    4 1.763831 1.2664840 0.14194413 1.779472e-01
## 5    5 1.456391 1.1950296 0.15341861 1.441929e-01
## 6    6 1.046582 0.8814642 0.17189460 1.277271e-02
## 7    7 1.325178 1.0732927 0.25014617 7.429692e-04
## 8    8 1.403162 1.0330854 0.27993645 6.131462e-03
## 9    9 1.096374 1.0049492 0.31254871 2.976572e-07
## 10  10 1.667162 1.1013899 0.35280593 5.732081e-05
## 11  11 1.591532 1.1700133 0.32756392 2.486194e-03
## 12  12 1.343440 0.9537399 0.32944838 1.151621e-03
## 13  13 1.279023 0.7870739 0.41499978 1.454529e-02
## 14  14 1.580736 1.0510679 0.63860877 7.219584e-03
## 15  15 1.645274 1.0178021 0.59534610 3.064336e-03
xyplot(eqm + Vy + Vh + Bh ~ deg,
       data = resul,
       type = c("p", "o"),
       xlab = "Grau do polinômio",
       ylab = "Componentes do erro quadrático médio",
       auto.key = list(corner = c(0.95, 0.95),
                       text = c(expression(EQM(hat(y))),
                                expression(Var(y)),
                                expression(Var(hat(y))),
                                expression(Bias(hat(y))))))#  +

    # glayer(panel.smoother(..., se = FALSE))
#-----------------------------------------------------------------------
# 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) <- sprintf("Pol. grau: %d", 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")
            })

6 Leave-one-out, DFfits e DFbetas

Para um material mais completo, visite: http://leg.ufpr.br/~walmes/cursoR/fiocruz/fiocruz04diagn.html.

# 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))
##    (Intercept)        speed   I(speed^2)
## 1 -0.274891030  0.246714807 -0.218599848
## 2  0.109198833 -0.098005996  0.086837495
## 3 -0.188270827  0.147739655 -0.117691200
## 4  0.158683853 -0.124522200  0.099195894
## 5 -0.002372816  0.001681285 -0.001215252
## 6 -0.080020109  0.046278022 -0.025675268
head(t(sapply(loo, "[[", "dfbetas")))
##       (Intercept)        speed   I(speed^2)
## [1,] -0.274891030  0.246714807 -0.218599848
## [2,]  0.109198833 -0.098005996  0.086837495
## [3,] -0.188270827  0.147739655 -0.117691200
## [4,]  0.158683853 -0.124522200  0.099195894
## [5,] -0.002372816  0.001681285 -0.001215252
## [6,] -0.080020109  0.046278022 -0.025675268
head(dffits(m0))
##            1            2            3            4            5 
## -0.281893050  0.111980344 -0.223308547  0.188215355 -0.003221695 
##            6 
## -0.137841851
head(sapply(loo, "[[", "dffits"))
##            1            2            3            4            5 
## -0.281893050  0.111980344 -0.223308547  0.188215355 -0.003221695 
##            6 
## -0.137841851
Machine Learning Prof. Eduardo V. Ferreira & Prof. Walmes M. Zeviani
25px