Simulando nossa incerteza sobre um modelo

Como vimos na primeira aula nossos modelos possuem dois componentes, o estocástico e o sistemático:

\[ Y_{i} \sim f(y_{i}|\theta_{i},\alpha) \\ \theta_{i} = g(X, \beta) \]

O componente estocástico descreve nossa incerteza sobre como se distribui nossa variável dependente. Podemos ter uma ideia de qual a probabilidade de um projeto de lei ser vetado, ou de um filho ter uma ocupação de maior status que o pai, mas, se não for zero ou um, não sabemos exatamente qual é esta probabilidade. Parte dessa incerteza é capturada pelo componente sistemático, mas sempre há um elemento irredutível na nossa incerteza quanto ao componente estocástico. Por outro lado os parâmetros de nosso componente sistemático também são estimados com algum grau de incerteza que, utilizando as técnicas vistas nas aulas anteriores, procuramos reduzir ao máximo.

Por conta da incerteza irredutível e da incerteza quanto à estimativa dos parâmetros não devemos apresentar apenas as estimativas de ponto, ou seja os coeficientes ou valores preditos obtidos, pois isso passa a ideia de uma falsa precisão. Devemos sempre apresentar as inferências ou predições derivadas de nossos modelos com seus intervalos de confiança. Em geral apresentamos os coeficientes com um intervalo de +/- 2 DP e a mesma coisa com os valores preditos. Mais difícil é apresentar esses intervalos de maneira combinada. Essa dificuldade ainda é maior se estamos trabalhando com modelos com interações ou com modelos mais complexos.

Podemos usar simulações para entender melhor nossos modelos. Simulações permitem incorporar nossas incertezas independente da complexidade do modelo. Como exemplo vamos usar o modelo do log da renda em função da altura e do sexo com interação entre as duas variáveis. (dados: Gelman e Hill survey de adultos nos EUA em 1994)

load("heights3.Rdata")
#attach(heights.clean)
log.earn <- log(heights.clean$earn)
earn.logmodel.3 <- lm (log.earn ~ height + male + height:male, data = heights.clean)
display (earn.logmodel.3, digits = 3)
## lm(formula = log.earn ~ height + male + height:male, data = heights.clean)
##             coef.est coef.se
## (Intercept)  8.388    0.844 
## height       0.017    0.013 
## male        -0.079    1.258 
## height:male  0.007    0.019 
## ---
## n = 1192, k = 4
## residual sd = 0.881, R-Squared = 0.09

Predição

A função predict() nos dá estimativas de ponto e os valores mínimo e máximo do intervalo de confiança de 95% para a variável dependente conforme ‘novos’ valores de nossas variáveis explicativas

x.new <- data.frame (height=68, male=1)
pred.interval <- predict (earn.logmodel.3, x.new, interval="prediction", 
  level=.95)

exp(pred.interval)
##        fit      lwr      upr
## 1 21435.57 3794.597 121088.9

Podemos construir um intervalo de predição em torno deste resultado representando nossa incerteza sobre o componente estocástico, ou seja sobre a distribuição predita de nossa variável dependente. Ao assumir que nossa VD segue uma distribuição gaussiana, ou normal, é simples construir um intervalo de predição. Para isso usamos a função rnorm(). Ela nos dá n valores extraídos de uma distribuição normal com média e desvio padrão definidos por nós. Neste caso pedimos 1000 valores com média 9.95 [por que?] e desvio padrão 0.88 [onde encontramos isso?] que extraímos do resultado do ajuste de nosso modelo

pred <- exp (rnorm (1000, 9.95, .88))
pred.original.scale <- rnorm (1000, 9.95, .88)

# Histograms (Figure 7.2)

par (mfrow=c(1,2))
hist (pred.original.scale, xlab="log(earnings)", main="")
hist (pred, xlab="earnings", main="")

mean(pred)
## [1] 31998.34
quantile(pred, 0.5)
##      50% 
## 21726.17
quantile(pred, c(0.25,0.75))
##      25%      75% 
## 11774.55 39136.84
quantile(pred, c(0.025,0.975))
##       2.5%      97.5% 
##   3722.509 122417.376

Simulações nos permitem fazer todo tipo de pergunta sobre nosso modelo, por exemplo, qual a diferença de renda entre homens e mulheres com a mesma altura?

pred.man <- exp (rnorm (1000, 8.4 + 0.017*68 - 0.079*1 + .007*68*1, .88))
pred.woman <- exp (rnorm (1000, 8.4 + 0.017*68 - 0.079*0 + .007*68*0, .88))
pred.diff <- pred.man - pred.woman
pred.ratio <- pred.man/pred.woman
mean(pred.diff)
## [1] 10528.8
mean(pred.ratio)
## [1] 3.094522

Incerteza em torno dos parâmetros

Outra fonte de incerteza ao estimar um modelo vem da estimação dos parâmetros. Podemos usar simulações para representar a incerteza nos coeficientes da regressão Aqui usamos a função sim() do pacote arm que calcula n simulações dos coeficientes estimados no modelo, isto é, nos dá a distribuição posterior do coeficiente:

n.sims <- 1000
sim.1 <- sim(earn.logmodel.3, n.sims)
height.coef <- sim.1@coef[,2]
mean (height.coef) # média do coeficiente 
## [1] 0.01691308
sd(height.coef) # o erro padrão do coeficiente
## [1] 0.01359873
quantile (height.coef, c(.025, .975)) # intervalo de 95%
##         2.5%        97.5% 
## -0.009449556  0.043983539
par (mfrow=c(1,1))
hist(height.coef, main = "Distribuição do coeficiente da altura", cex.main = 0.7) # visualizando a incerteza em torno do coeficiente

O que a função sim() faz é simular os coeficientes a partir de uma distribuição normal multivariada onde a média é o valor dos coeficientes e a variância é uma matriz de covariância dos coeficentes multiplicada pela variância do modelo.

# betas <- mvrnorm(1000, fit.1$coefficients, sd(fit.1$residuals)*vcov(fit.1))
# mean(betas[,2])
# sd (betas[,2])
# quantile (betas[,2], c(.025, .975))

Uma outra vantagem das simulações é ajudar é entender as interações

height.for.men.coef <- sim.1@coef[,2] + sim.1@coef[,4]
quantile (height.for.men.coef, c(.025, 0.5, .975))
##         2.5%          50%        97.5% 
## -0.002098528  0.023712548  0.051233205
hist(height.for.men.coef, main = "Distribuição do coeficiente da interação entre sexo e altura", cex.main = 0.7)

Predição incorporado incerteza quanto à distribuição da VD e quanto aos parâmetros

Agora que sabemos simular nossas incertezas quanto aos componentes estocástico e quanto ao componente sistemático de nosso modeo utilizando as funções rnorm() - ou qualquer outra distribuição -e sim(), podemos combinar as duas para fazer inferências ou predições para novos dados. Por exemplo podemos usar simulações para o caso dos poços em Bangladesh

wells <- read.table ("wells.dat")
fit.1 <- glm (switch ~ dist, data = wells, family=binomial(link="logit"))
display (fit.1, digits =4)
## glm(formula = switch ~ dist, family = binomial(link = "logit"), 
##     data = wells)
##             coef.est coef.se
## (Intercept)  0.6060   0.0603
## dist        -0.0062   0.0010
## ---
##   n = 3020, k = 2
##   residual deviance = 4076.2, null deviance = 4118.1 (difference = 41.9)

Primeiro verificamos a incerteza quanto aos coeficientes

n.sims <- 1000
sim.wells <- sim(fit.1, n.sims) # Fazemos 1000 simulações dos parâmetros (intercepto e efeito da distância)
head(sim.wells@coef)
##      (Intercept)         dist
## [1,]   0.6553858 -0.005996424
## [2,]   0.6386174 -0.007097600
## [3,]   0.6876117 -0.008107960
## [4,]   0.6006428 -0.006169664
## [5,]   0.5447238 -0.004429732
## [6,]   0.5657806 -0.005832619
mean(sim.wells@coef[,1])
## [1] 0.606741
mean(sim.wells@coef[,2])
## [1] -0.006247246

Vamos visualizar como se distribuem os pares de coeficiente

plot(sim.wells@coef[,1], sim.wells@coef[,2],
     xlab = expression(beta[0]),
     ylab = expression(beta[1]))

Agora visualizamos a incerteza em torno da estimativa do efeito da distância

plot(wells$dist, jitter(wells$switch,0.1), pch = 16, cex = 0.6,
     xlab = "Distância", ylab = "Prob. de Mudar de Poço", cex.lab = 0.7)
for(s in 1:n.sims){
  curve(invlogit(sim.wells@coef[s,1] + sim.wells@coef[s,2]*x),
        col = "gray",
        add = TRUE)
}
curve(invlogit(coef(fit.1)[1] + coef(fit.1)[2]*x), add = TRUE)

Suponha que temos 200 novas casas uniformemente espalhadas pelo território e queremos prever a chance deles mudarem para poços mais seguros se fizermos a mesma campanha

x.tilde <- cbind(rep(1,200),runif(200, 0, 300))
n.tilde <- nrow(x.tilde)
y.tilde <- array(NA, c(n.sims, n.tilde))
for(s in 1:n.sims){
  p.tilde <- invlogit(x.tilde %*% sim.wells@coef[s,])
  y.tilde[s,] <- rbinom(n.tilde, 1, p.tilde) # aqui usamos a distribuição binomial com tamanho 1, o que equivale a uma bernoulli, para nossa VD
}

y.tilde[1:5,1:5]
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    1    0    0    0
## [2,]    0    0    0    0    0
## [3,]    0    0    0    1    0
## [4,]    1    1    0    0    0
## [5,]    1    1    1    0    0
mean(rowMeans(y.tilde)) 
## [1] 0.430735
hist(rowMeans(y.tilde), main = "Distribuição das probabilidades preditas", cex.main = 0.7)

Agora observando como essa probabilidade muda com a distância

plot(x.tilde[,2], jitter(y.tilde[1,], 0.1), xlab = "Distância", ylab = "Prob. de mudar de poço", cex.lab = 0.7, pch = 16, cex = 0.6)
for(s in 1:n.sims){
  curve(invlogit(sim.wells@coef[s,1] + sim.wells@coef[s,2]*x),
        col = "gray",
        add = TRUE)
}
curve(invlogit(coef(fit.1)[1] + coef(fit.1)[2]*x), add = TRUE)

Usando simulação para verificar ajuste do modelo aos dados

Podemos usar nossas observações como se fossem os “novos” dados para verificar o ajuste de nosso modelo aos dados. Se a distrobuição dos valores preditos for muito diferente da distribuição de nossas observações o modelo não é muito bom.

n <- length(wells$switch)
X <- cbind(rep(1,n), wells$dist)
#head(X)
#head(model.matrix(fit.1))
y.tilde <- array(NA, c(n.sims, n))
for(s in 1:n.sims){
  p.tilde <- invlogit(X %*% sim.wells@coef[s,])
  y.tilde[s,] <- rbinom(n, 1, p.tilde)
}

mean(wells$switch)
## [1] 0.5751656
mean(rowMeans(y.tilde))
## [1] 0.5750934
hist(rowMeans(y.tilde), main = "Distribuição das probabilidades preditas", cex.main = 0.7)
abline(v = mean(wells$switch), lwd = 2, col = 2)

par(mfrow = c(1,2))
barplot(prop.table(table(wells$switch)), main = "Frequência das probabilidades observadas", cex.main = 0.7)
barplot(prop.table(table(y.tilde)), main = "Frequência das probabilidades preditas", cex.main = 0.7)

Bootstrap

Mais acima usamos a função sim() para simular a incerteza quanto aos nossos parâmetros. Essa incerteza foi calculada utilizando o erro padrão que obtivemos do ajuste do modelo em toda a base de dados. Outra maneira de quantificar nossa incerteza quanto aos parâmetros é calculá-los utilizando amostras de nossos dados sorteadas com reposição. Isso equvale a treinarmos nosso modelo em vários dados de treino.

library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:arm':
## 
##     logit
boot.fn <- function(dados, indices) return(coef(glm(switch ~ dist, data = wells, family=binomial(link="logit"),subset = indices)))
ene <- nrow(wells)
boot.fn(wells,1:ene)
##  (Intercept)         dist 
##  0.605959360 -0.006218819
coef(fit.1)
##  (Intercept)         dist 
##  0.605959360 -0.006218819
set.seed(1)
boot.fn(wells,sample(ene,ene, replace = T))
##  (Intercept)         dist 
##  0.640813760 -0.005999299
wells.boot <- boot(wells, boot.fn,1000)

# Vamos visualizar como se distribuem os pares de coeficiente

plot(wells.boot$t[,1], wells.boot$t[,2],
     xlab = expression(beta[0]),
     ylab = expression(beta[1]))

# Agora visualizamos a incerteza em torno da estimativa do efeito da distância

plot(wells$dist, jitter(wells$switch,0.1), xlab = "Distância", ylab = "Prob. de mudar de poço", cex.lab = 0.7, pch = 16, cex = 0.6)
n.sims <- wells.boot$R
for(s in 1:n.sims){
  curve(invlogit(wells.boot$t[s,1] + wells.boot$t[s,2]*x),
        col = "gray",
        add = TRUE)
}
curve(invlogit(coef(fit.1)[1] + coef(fit.1)[2]*x), add = TRUE)

Simulando modelos mais complexos

Simulações são o meio mais fácil para interpretarmos modelos mais complexos. Uma situação em que precisamos de um modelo mais complexo ocorre quando, por exemplo, temos uma mistura de dados qualitativos e quantitativos. Se, por exemplo, estamos interssados em fatores associados ao nível de coelesterol no sangue, mas nossa base só traz essa informações para os sujeitos que tinham colesterol acima de dado valor, podemos só trabalhar com os dados quantitativos, tirando da base os que tinham colesterol baixo, ou podemos de algum modo usar essa informação. Em economia modelos deste tipo são conhecidos como tobit.

Vamos deixar isso mais claro com um exemplo prático retirado de Gelman e Hill (6.7 e 7). Na base heights.clear que usamos acima só temos as observações com renda maior que zero. Deste modo não sabemos nada sobre a realação entre sexo, altura e renda para aqueles que não tem renda. Podemos usar a base com todas as observações (heights2) e usar dois estágios para poder incorporar a informação se o sujeito tem ou não renda (uma variável dicotômica).

Primeiro ajustamos um modelo logit para calculara probabilidade de um sujeito ter renda:

h2 <- heights2
earn.pos <- ifelse(h2$earn > 0, 1, 0)
earn.fit <- glm(earn.pos ~ height + male, family = binomial(link = "logit"), data = h2)
display(earn.fit)
## glm(formula = earn.pos ~ height + male, family = binomial(link = "logit"), 
##     data = h2)
##             coef.est coef.se
## (Intercept) -3.76     2.07  
## height       0.08     0.03  
## male         1.70     0.32  
## ---
##   n = 1379, k = 3
##   residual deviance = 989.9, null deviance = 1094.7 (difference = 104.8)

Depois ajustamos um modelo log-linear na renda

log.earn <- log(h2$earn)
earn.fit2 <- lm(log.earn ~ height + male, data = h2, subset = earn > 0)
display(earn.fit2)
## lm(formula = log.earn ~ height + male, data = h2, subset = earn > 
##     0)
##             coef.est coef.se
## (Intercept) 8.15     0.60   
## height      0.02     0.01   
## male        0.42     0.07   
## ---
## n = 1192, k = 3
## residual sd = 0.88, R-Squared = 0.09

Esses modelos nos dizem, por exemplo, que uma mulher de 1.60 metros (63 polegadas) tem uma probabilidade de \(\text{invlogit}(-3.85 + 0.08*66 + 1.7*0) = 0.80\) ou 80% de ter alguma renda. Se ela possui renda, o valor predito desta renda é de \(\text{exp}(8.15 + 0.02*66 + 0.42*0) = 12964\). Um maneira de combinar estas duas informações é usando simulações:

x.novo <- c(1, 68, 1) # (constante, altura, homem)

n.sims <- 1000

p.earn.pos <- invlogit(coef(earn.fit) %*% x.novo)
earn.pos.sim <- rbinom(n.sims,1,p.earn.pos)
earn.sim <- ifelse(earn.pos.sim == 0, 0, exp(rnorm(n.sims, coef(earn.fit2) %*% x.novo, sigma.hat(earn.fit2))))
quantile(earn.sim, c(0.025,0.25,0.5,0.75,0.975))
##       2.5%        25%        50%        75%      97.5% 
##   2532.111  11234.087  21526.859  38380.247 119047.923

Podemos incorporar nossa incerteza sobre os parâmetros

sim.1a <- sim(earn.fit, n.sims)
sim.1b <- sim(earn.fit2, n.sims)

p.earn.pos <- invlogit(sim.1a@coef  %*% x.novo)
earn.pos.sim <- rbinom(n.sims,1,p.earn.pos)
earn.sim <- ifelse(earn.pos.sim == 0, 0, exp(rnorm(n.sims, sim.1b@coef  %*% x.novo, sim.1b@sigma)))
quantile(earn.sim, c(0.025,0.25,0.5,0.75,0.975))
##      2.5%       25%       50%       75%     97.5% 
##      0.00  11261.76  19570.86  34185.19 108444.82

Se quisermos visualizar a relação entre renda e diferentes valores de altura e sexo podemos jogar tudo isso em uma função que nos retorna a média da renda (earn) para cada combinação:

rendaMedia <- function(altura, homem, sim.a, sim.b){
  x.novo <- c(1, altura, homem)
  p.earn.pos <- invlogit(sim.a@coef  %*% x.novo)
  earn.pos.sim <- rbinom(n.sims,1,p.earn.pos)
  earn.sim <- ifelse(earn.pos.sim == 0, 0, exp(rnorm(n.sims, sim.b@coef  %*% x.novo, sim.b@sigma)))
  return(mean(earn.sim))
}

Para aplicar essa função a uma sequência de vlores de altura e sexo usamos a função sapply()

altura <- seq(60, 75, 1)
renda.mulher <- sapply(altura, rendaMedia, homem = 0, sim.1a, sim.1b)
renda.homem <- sapply(altura, rendaMedia, homem = 1, sim.1a, sim.1b)

plot(altura, renda.mulher, type = "l", ylim = c(10000,40000), ylab = "Renda", cex.lab = 0.8)
lines(altura, renda.homem, col = "blue")
text(72,35000, "Homens", cex = 0.6)
text(72,15000, "Homens", cex = 0.6)

Juntando tudo: 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.

Na aula passada vimos que, 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.

O ideal seria termos um métdo que ao mesmo tempo 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)
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]))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

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)
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 aos 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 \(\lamda\) 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)
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