Procurando \(f\)

O objetivo do aprendizado estatístico é estimar \(f(x)\), a forma funcional do componenete sistemático. Procuramos o método estatístico que nos dê um \(\hat{f}(x)\) que mais se aproxime de \(f(x)\). Nas modelagens que vimos até agora assumimos uma forma funcional determinada para \(f(x)\) e, a partir desta premissa, estimamos os parâmetros desta função por MQO, por ML ou por algum outra técnica. Modelos deste tipo são chamados de paramétricos e podem se aplicar tanto a dados quantitativos, quando são chamados de regressão, quanto a dados categóricos, quando são chamados de modelos de classificação.

Outra estratégia de modelagem não assume uma forma funcional para o componente sistemático do modelo. Procuramos nos aproximar o melhor possível dos dados utilizando esses mesmos dados como guia, evitando formas funcionais muito complexas.

Modeos paramétricos tema avantagem de serem mais interpretáveis, o que permite seu uso para infierência (quais são os preditores suficientes para nosso modelo, qual a relação entre a variável resposta e os preditores, et.). A desvantagem é que são menos flexíveis, isto é, se a forma funcional for muito diferente da forma assumida o erro será grande.

Modelos não paramétricos são mais flexíveis e são mais precisos, o que os tornam muito úteis para predição. A desvantagem deles é que são pouco interpretáveis. Há sempre um trade-off entre interpretabilidade e flexibilidade.

Como escolher o melhor método?

Ao ajustar um modelo estamos interessados em quão bem os valores preditos correspondem aos valores observados. Calculamos a qualidade do ajuste por meio do Erro Médio Quadrado, no caso de VD quantitativa ou regressão, e da Taxa de Erro para VDs qualitativas ou classificação:

\[ \text{EMQ} = \frac{1}{n} \sum_{i=1}^{n}(y_{i} - \hat{f}(x_{i}))^2\]

e

\[ \text{TE} = \frac{1}{n} \sum_{i=1}^{n} I(y_{i} \neq \hat{y_{i}})\]

Como vimos acima podemos tornar nossos modelos flexíveis o bastante para se ajustar quase que perfeitamente aos dados. O problema é que, quando trabalhamos com predição, não nos interessa tanto que nosso modelo se ajuste às observações nas quais ele foi treinado, mas sim a observações ainda não vistas. Não nos interessa, por exemplo, saber a probabilidade de um golpe de estado que já aconteceu (a probabilidade é 1), mas sim conseguir prever se haverá um novo golpe ou não. Também não nos interessa apenas descrever os determinantes do desempenho passado de alunos, mas sim qual será o desempenho conforme variações nestes determinantes.

Um modelo bem flexível se ajusta muito bem aos nossos dados, nos dando EQM de treino. Mas se usarmos esse modelo em novos dados, seu ajuste não será muito bom, isto é o EQM nestes novos dados será grande.

Se tivermos dados o suficiente podemos usar nossa estimativa \(\hat f\) para calacular o MSE ou o TE de observações teste e escolher o método que nos dá o menor erro no teste.

Podemos decompor o EQM de teste em:

\[ \mathbb{E}(y_o - \hat{f}(x_o))^2 = \text{Var}(\hat{f_o}) + \text{Viés}(\hat{f}(x_o))^2 + \text{Var}(\epsilon) \]

Modelos mais flexíveis tem maior variância na forma funcional e menor viés. O contrário ocorre com modelos menos flexíveis. Isso faz com tenhamos um trade-off entre viés e variância.

A partir de certa quantidade de flexibilidade, o MSE do teste começa a subir apesar do MSE do treino diminuir continuamente. Quando o MSE treino é pequeno, mas o MSE teste é grande temos o que se chama de overfitting. A partir de certo ponto nosso modelo começa a capturar padrões nos dados que ocorrem ao acaso e não tem nada a ver com o processo ‘verdadeiro’ que gerou aqueles dados. Quando ajustamos esse modelo ao teste aqueles padrões aleatórios não estão presentes e o ajuste não é tão bom. Esse padrão de “U” no MSE de teste ocorre sempre em qualquer base de dados e em qualquer modelo.

Portanto devemos procurar ajustar nosso modelo a novos dados para evitar overffiting e garantir validade externa a eles. No entanto, muitas vezes não temos dados o suficiente para separar uma base de teste. Quando isso ocorre podemos usar um conjunto de técnicas que permitem extrair várias bases teste de um mesmo conjunto de dados. Uma destas técnicas é o chamado bootstrep, quando utilizando amostras de nossos dados sorteadas com reposição. Isso equvale a treinarmos nosso modelo em vários dados de treino. Outro conjunto de técnicas recebem o nome de cross validation ou validação cruzada.

Na validação simples fazemos como no exemplo que acabamos de ver no R: sorteamos aleatoriamente algumas observações de nossos dados e usamos para treino, o resto será a base teste. Neste proceso acabamos com aenas uma base de teste. Na validação cruzada tentamos maximizar o uso dos dados para validação. No limite temos a Leave-One-Out Cross-Validation (LOOCV) onde separamos uma observação para teste e usamos todas as outras para treinar o modelo, repetindo o procedimento para todas as observações. Esse processo exige muita computação. Um métod mais eficiente é o K-Fold Cross-Validation.

No k-fold CV dividimos, aleatoriamente, nossa base em k grupos (folds). Um grupo servirá como teste e os outros para treino, repetindo para todos os grupos. Para cada teste obtemos um MSE.Tiramos a média de todos esses k’s MSE e obtemos o erro quadrado médio do k-fold CV.

No R:

library(ISLR)
data("Auto")
library(boot)
data(Auto)
set.seed(17)
cv.erro.10 <- rep(0,10)
for(i in 1:10){
  glm.fit <- glm(mpg~poly(horsepower, i), data = Auto)
  cv.erro.10[i] <- cv.glm(Auto, glm.fit, K = 10)$delta[1]
}
cv.erro.10
##  [1] 24.27207 19.26909 19.34805 19.29496 19.03198 18.89781 19.12061 19.14666
##  [9] 18.87013 20.95520

Modelos não paramétricos de classificação

Os modelos que vimos até agora foram estimados utilizando métodos paramétricos. Nós estabelecemos uma forma funcional para nosso modelo que tenta se aproximar da verdadeira forma da distribuição \(f(Y \mid X)\) e então estimamos parâmetros, isto é, valores que sintetizam a relação entre \(Y\) e \(X\). Se a forma funcional que estabelecemos para nosso modelo for uma distribuição normal, por exemplo, estimamos a média e variância e temos assim uma noção de como \(Y\) varia com \(X\), do mesmo modo se nossa distribuição for da família da binomial estimando \(\pi\) obtemos uma decrição de \(f(Y\mid X)\).

Um possível problema do método paramétrico é que o modelo que usamos nem sempre reflete corretamente a verdadeira forma de \(f(Y \mid X)\). Um outro método para conhecermos a distribuição de \(f(Y\mid X)\) não parte de uma forma funcional, mas tenta apreender esta forma diretamente dos dados. Estes métodos tentam se aproximar o máximo possível de todos os pontos dos dados evitando, ao mesmo tempo, variar muito.

Um exemplo de método não paramétrico é, como vimos, o uso de lowess para estimar o voto no plebiscito no Chile. Existem vários outros métodos não paramétricos ou semi-paramétricos para VD’s tanto quantitativas quanto qualitativas. Um deles é o K-Nearest Neighbors (KNN). Por este método atribuimos a observação \(i\) à categoria A ou B com uma probabilidade dada pela quantidade de vizinhos que pertencem a uma ou outra categoria.

\[ \text{Pr}(Y = j\mid X = x) = \frac{1}{K} \sum_{i \in \mathcal{N(x)}} I(y_i = j)\]

Onde \(j\) é uma das categorias e \(\mathcal{N}\) é a vizinhança de \(i\). O KNN atribui \(x\) à categoria \(j\) com a maior probabilidade de ocorrência entre os vizinhos.

Observem que o KNN é uma abordagem completamente não paramétrica e não nos diz quais variáveis são importantes, não obtemos nada parecido com uma tabela de coeficientes como vimos quando usamos outros modelos, isto é, um modelo pouco interpretável. Esse tipo de modelo é usado quando estamos interessados mais em predição do que em inferência. Neste último caso estamos interessados na relação entre a variévale dependente e as variáveis explicativas, isto é, de que forma um determinado preditor, como a adesão ao status quo no caso do plebiscito, impacta na VD, o voto no sim ou não. Observamos os coeficientes para testar nossas teorias sobre essas relações, nossos modelos teóricos. Na predição estamos interessados em estimar futuros eventos com base nos dados e algumas suposições sobre como aqueles dados surgiram. No KNN, por exemplo, com base no conhecimento que adquirirmos ao observar dados da relação entre \(Y\) e \(X\), conseguims predizer, com algum erro, as categorias a que pertencerão um novo \(\tilde X\). Vamos falar mais de predição na próxima aula.

Outro modelo não paramétrico muito utilizado são as Árvores de Decisão.

Métodos como MQO dão resultados excelentes produzindo estimadores não enviesados e com a menor variância. O método da máxima verossimilhança também garante estiamdores com mínima variância e não viesados e, dependendo de algumas condições de regularidade, garante que os estimadores são assimptoticamente eficientes, isto é consistentes (\(\hat \theta_{\text{plim}} = \theta\)) e normalmente distribuídos (\(\hat \theta \sim N[\theta, \textbf{I}(\theta)^{-1}]\)). Mas essas propriedades dependem de que a forma funcional que assumimos esteja correta.

Recentemente e graças ao aumento no poder de processamento dos computadores, métodos mais flexíveis tem sido propostos para estimar modelos. Vimos que métodos não-paramétricos não assumem uma forma funcional para a relação entre a variável resposta e as variáveis explicativas. Estes métodos buscam essa forma funcional nos dados. Quando a relação entre nossa VD e nossas variáveis explicativas é complexa esses métodos dão resultados superiores, mensurados pos alguma função de erro.

Vimos também que reduzir o erro em um conjunto de dados não garante que o modelo se ajuste bem a futuras observações. Na verdade modelos muito flexíveis podem acabar se sobre-ajustando aos dados gerando erros maiores quando ajustado a novas observações. Isso se deve a um compromisso entre viés e variância. Ao flexibilizar nossos modelos para diminuir o viés de nossa predição aumentamos sua variância.

Ao treinar nossos modelos em diferentes bases conseguimos atingir um equilíbrio entre viés e variância. Ajustamos modelos a uma parte dos dados, aprendemos a forma funcional \(\hat f(Y\mid X)\), testamos esse ajuste em outra base de dados, que chamamos de teste, verificamos a qualidade do ajuste e selecionamos o modelo com o melhor ajuste. Vimos que cross-validation e bootstrap permitem aproveitar ao máximo nossos dados para treinar e testar várias vezes nossos modelos.

Alémde um modelo que otimize a relação entre viés e variância, o ideal seria termos um métdo que tenha boa flexibilidade, para poder se ajustar a qualquer relação entre \(Y\) e \(X\) e, ao mesmo tempo, ser facilmente interpretável. Uma regressão linear é de fácil interpretação, mas pouco flexível. Um KNN é muito flexível, mas pouco interpretável. Um bom compromisso entre esses dois objetivos nos é dado pela técnica das Árvores de Decisão.

Vamos supor que queremos estimar a probabilidade de alguém votar em Bolsonaro dadas certas caracteríticas pessoais (sexo, renda, escolaridade, religião etc.). Podemos proceder testando cada uma das variáveis explicativas e vendo qual minimiza o erro de classificação. Comecemos pela renda. Dividimos a renda em dois grupos, os que ganham mais que x e os que ganham menos que x. Agora contamos quantos sujeitos em cada grupo votou no Bolsonaro. Se naquele grupo a maioria votou em Bolsonaro, classifico todo grupo como bolosnarista e calculo qual seria meu erro com relação à distribuição real. Podemos fazer isso para várias divisões e verificar quais geram o menor erro de classificação. Fazemos a mesma coisa com as outras variáveis. Escolhemos a variável \(k\) e o ponto da divisão \(c\) que nos dá o menor erro e dividimos nosso espaço amostral neste ponto. Agora temos dois “conjuntos” de dados, uma parte cujos valor de \(k\) é menor que \(c\) e outra onde é maior. Repito esse procedimento nesses dois subconjuntos e vou, dessa maneira, criando vários subgrupos conforme variáveis e pontos de corte que minimizam meu erro de classificação.

load("df3.RData")
library(tree)
## Registered S3 method overwritten by 'tree':
##   method     from
##   print.tree cli
bozo.tree <- tree(bozo ~ idade + sexo + ideo + evan + esc + rend + regiao, df3)

summary(bozo.tree)
## 
## Classification tree:
## tree(formula = bozo ~ idade + sexo + ideo + evan + esc + rend + 
##     regiao, data = df3)
## Variables actually used in tree construction:
## [1] "ideo" "sexo" "rend" "esc" 
## Number of terminal nodes:  5 
## Residual mean deviance:  1.104 = 1098 / 995 
## Misclassification error rate: 0.308 = 308 / 1000
plot(bozo.tree)
text(bozo.tree, cex = 0.7, label = "yprob")

plot(bozo.tree)
text(bozo.tree, cex = 0.7)

plot(bozo.tree)
text(bozo.tree, pretty = 0, cex = 0.7)

bozo.tree
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 1000 1312.0 0 ( 0.63500 0.36500 )  
##    2) ideo: esq 230  116.2 0 ( 0.93043 0.06957 ) *
##    3) ideo: spart,dir 770 1061.0 0 ( 0.54675 0.45325 )  
##      6) sexo: Feminino 427  552.4 0 ( 0.65105 0.34895 )  
##       12) rend: Até 2 S.M.,Não sabe 200  220.4 0 ( 0.76000 0.24000 )  
##         24) esc: Analfabeto/ primario incompleto 27    0.0 0 ( 1.00000 0.00000 ) *
##         25) esc: Primario completo/ Ginasial incompleto,Ginasial completo,Colegial incompleto,Colegial completo,Superior incompleto,Superior completo,Pós graduação 173  204.3 0 ( 0.72254 0.27746 ) *
##       13) rend: De 2 a 3 S.M.,De 3 a 5 S.M.,De 5 5 a 10 S.M.,De 10 a 20 S.M. 227  311.9 0 ( 0.55507 0.44493 ) *
##      7) sexo: Masculino 343  466.0 1 ( 0.41691 0.58309 ) *

As várias divisões nos dados podem ser ilustradas por meio de uma árvore de decisão, daí o nome do método. Na parte superior temos a variável com maior impacto e o ponto de sua divisão que gerou o menor erro. No caso do voto em Bolsonaro a variável de maior impacto é ideologia. O ponto de divisão é “esquerda” e vemos que a probabilidade de alguém votar em bolonaro se identificando com a esquerda é de 7%. Vemos que aqueles que não se identificam com partidos ou se identificam com a direita tem maior probabilidade de votar em Bolsonaro, mas aqui há outras divisões. A primeira divisão se deve ao sexo. Mulheres tem menor probabilidade de votar em Bolsonaro (42%), mas elas se dividem. A renda é a variável com mais impacto nesta terceira divisão. Aqueles que sabem sua renda e possuem renda maior que 2 SM tem 55% de probabilidae de votar em bolsonaro. Os que têm renda inferior a 2 SM ou que não sabem sua renda se dividem entre aqueles com baixa escolaridade, com 0% de probabilidade de votar em Bolsonaro e aqueles com maior escolaridade, com 44% de probabilidade de votar nele.

As Árvores de Decisão são bem fáceis de interpretar e nos dão de forma direta a relevância de cada variável na predição de uma variável resposta. A extensão dele para o caso de variáveis resposta quantitativas é direto. A diferença é que no lugar de minimizar o erro de classificação minimizamos a soma dos resíduos ao quadrado (RSS). Podemos calcular o erro de nossa árvore ajustando os dados em uma base de treino e prevendo resultados em uma base de teste

set.seed(2)
treino <- sample(1:nrow(df3), 600)
df3.teste <- df3[-treino,]
bozo.teste <- df3.teste$bozo
bozo.tree2 <- tree(bozo ~ idade + sexo + ideo + evan + esc + rend + regiao, df3, subset = treino)
bozo.pred <- predict(bozo.tree2, df3.teste, type = "class")
table(bozo.pred, bozo.teste)
##          bozo.teste
## bozo.pred   0   1
##         0 201  88
##         1  47  64

Nosso acerto de predição é de 73%. Vamos examinar a tabela de erro, chamada de tabela de confusão. As linhas correspondem às categorias preditas e as colunas aos valores ‘verdadeiros’. A diagonal principal nos dá as previsões corretas e sua soma sobre o total nos dá a exatidão do modelo (accuracy). Os verdadeiros positivos sobre o total de previsões de positivos nos dá a precisão (precision) do modelo. Quando classificamos as observações podemos cometer dois tipos de erro. Podemos classificar um eleitor de Bolsonaro como tendo votado em outros (erro tipo I) ou podemos classificar quem votou em outro candidato como tendo votado em Bolsonaro (erro tipo II). Deixamos de indetificar 83 bolsonarista e classificamos errado 36 não bolosnaristas. Portanto enquanto o erro geral não é tão alto, dos 139 bolsonaristas da amostra identificamos apenas 40%. Essa razão entre os verdadeiros positivos e todos os realmente positivos nos dá a sensibilidade (sensitivity) do teste, ou o poder do teste. Por outro lado dos 261 não bolsonaristas da amostra classificamos errado apenas 14% esse valor nos dá a especificidade (specificity) do teste. Podemos aumentar a sensibilidade estabelecendo um limite menor para classificarmos um bolosnarista enquanto tal.

bozo.pred2 <- predict(bozo.tree2, df3.teste, type = "vector")
table(ifelse(bozo.pred2[,2] > 0.3,1,0), bozo.teste)
##    bozo.teste
##       0   1
##   0 143  34
##   1 105 118
#plot(roc(bozo.teste,bozo.pred2[,2]))

Por mais intuitivo que seja e por mais que possamos obter bons ajustes o método das Árvores não é tão preciso quanto o da Máxima Verossimilhança ou o MQO. Uma única árvore é muito sensível a escolha da primeira variável e ponto de corte, podendo dar resultados enviesados e pode acabar sobre-ajustando os dados. Uma maneira de suprir as deficiências das Árvores de Decisão é ajustar várias árvores aos dados usando bootstrap e depois agregar os resultados, tirando a média no caso de VD quanti ou escolhendo a classificação de maior frequência no caso de VD quali. A esse procedimento se dá o nome de bootstrap aggregation ou Bagging. Como não temos mais apenas uma árvore não temos a representação gráfica das Árvores de Decisão nem a mesma facilidade de interpretar os resultados. No entanto, é possível verificar quais variáveis são as mais importantes observando todos as divisões em todas as árvores e o quanto elas reduziram o erro, ponderado pelo número de árvores em que foi utilizada.

Embora o bagging resolva o problema de ajustar uma única árvore, ele não ajuda a resolver o problema da sensibilidade às condições iniciais. Uma maneira de garantir que várias árvores diferentes sejam treinadas é selecionar diferentes subconjustos de parâmetros e ajustar várias árvores com esses parâmetros. Isso evita que um dos parâmetros acabe dominando todas as árvores. Por termos várias árvores diferentes esse método se chama Random Forrest. Quando o RF usa todos os parâmetros ele se iguala ao bagging. Em geral ele usa \(m = \sqrt{p}\) parâmetros.

library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:ggplot2':
## 
##     margin
set.seed(1)
bozo.rf <- randomForest(bozo ~ idade + sexo + ideo + evan + esc + rend + regiao + ideo, df3, subset = treino, importance = TRUE)
bozorf.pred <- predict(bozo.rf, df3.teste)
table(bozorf.pred, bozo.teste)
##            bozo.teste
## bozorf.pred   0   1
##           0 202  89
##           1  46  63
importance(bozo.rf)
##                 0          1 MeanDecreaseAccuracy MeanDecreaseGini
## idade  -0.3874372  4.4049873             2.627281        28.827368
## sexo    9.6247756 12.2619883            15.269305        15.831632
## ideo   25.5018499 37.5519361            40.832620        38.849229
## evan    1.9536450  0.4687578             1.793893         9.817236
## esc     2.7999293  2.4346582             3.902736        36.828790
## rend    7.3016211 14.1051668            14.464518        34.037944
## regiao -2.5125511  8.6561861             3.682514        23.452778
varImpPlot(bozo.rf)

Outra maneira de usar várias árvores se dá pela técnina do Boosting. Neste método inicializamos os valores preditos e os resíduos em \(\hat f(x) = 0\) e \(r_{i} = 0\). Em cada iteração \(b\) ajustamos uma árvore com \(d\) divisões dos dados treino \((X,r)\) obtendo \(\hat f^b(x)\) e atualizamos \(\hat f(x) \leftarrow \hat f(x) + \lambda \hat f^b(x)\) e os resíduos \(r_{i} \leftarrow r_{i} - \lambda \hat f^b(x)\), isto é usamos parte \(\lambda\) do valor predito para a atualização. Depois de fazer isso em todas as iteraçãoes obtemos o valor predito \(\hat f(x) = \sum_{b = 1}^{B} \lambda \hat f^b(x)\). Em outras palavras ajustamos pequenas árvores aos resíduos reduzindo esses resíduos aos poucos até obter uma árvore com o menor resíduo. Ao ajustar poucas árvores pequenas (pequeno \(d\)) e ao “descontar” o impacto do ajuste (\(\lambda\)) o algoritmo aprende divagar “visitando” árvores que de outro modo não seriam consideradas.

library(gbm)
## Loaded gbm 2.1.8
set.seed(1)
bozo2 <- as.logical(as.integer(df3$bozo)-1)
df3$bozo2 <- bozo2
bozo.boost <- gbm(bozo2 ~ idade + sexo + ideo + evan + esc + rend + regiao, data = df3[treino,],
                  distribution = "bernoulli",
                  n.trees = 100,
                  interaction.depth = 1,
                  shrinkage = 0.1)
summary(bozo.boost)

##           var   rel.inf
## ideo     ideo 36.887806
## rend     rend 19.297919
## esc       esc 17.958458
## sexo     sexo 14.254870
## regiao regiao  6.436398
## idade   idade  3.002641
## evan     evan  2.161909
plot(bozo.boost, i = "ideo")

bozoboost.pred <- predict(bozo.boost, df3.teste, type="response")
## Using 100 trees...
table(ifelse(bozoboost.pred > 0.5,1,0), bozo.teste)
##    bozo.teste
##       0   1
##   0 201  81
##   1  47  71

Como base de comparação vamos calcular o erro de predição ajustando um logit aos dados

bozo.glm <- glm(bozo ~ idade + sexo + ideo + evan + esc + rend + regiao, data = df3[treino,], family = binomial(link = "logit"))
bozoglm.pred <- predict(bozo.glm, df3.teste, type="response")
table(ifelse(bozoglm.pred > 0.5,1,0), bozo.teste)
##    bozo.teste
##       0   1
##   0 197  75
##   1  51  77