Econometria III

.title[
# Econometria III
]
.subtitle[
## Modelos de Escolha Qualitativa
]
.author[
### Rafael Bressan
]
.date[
### Esag </br> 08/11/2023
]

---

---

# Onde estamos?

**Hoje**

1. Modelos de escolha qualitativa!

2. Um aplicativo legal! 😎

---

# Modelos de Escolha Qualitativa

---

---

# Modelos de Escolha Qualitativa

Até agora, nossos modelos ficaram assim:

`$$\begin{align}
y &= b_0 + b_1 x + e \\
e &\sim D\left(0,\sigma^2\right)
\end{align}$$`

* A suposição de distribuição em `$e$`:

* Em princípio implica que `$y \in \mathbb{R}$`.

* Resultados de testes, renda, taxas de criminalidade, etc. são todos resultados contínuos. ✅
]

Mas alguns resultados são claramente binários (ou seja, `VERDADEIRO` ou `FALSO`):

* Ou você trabalha ou não,

* Ou você tem filhos ou não,

* Ou você comprou um produto ou não,

* Você jogou uma moeda e saiu cara ou coroa.
]

---

# Resultados Binários

* Resultados restritos a `FALSO` vs `VERDADEIRO`, ou `0` vs `1`.

* Teríamos `$y \in \{0,1\}$`.

* Nessas situações, estamos principalmente interessados em estimar a **probabilidade de resposta** ou a **probabilidade de sucesso**:

`$$p(x) = \Pr(y=1 | x)$$`
* como `$p(x)$` muda quando mudamos `$x$`?

>Se aumentarmos `$x$` em uma unidade, como a probabilidade de `$y=1$` mudaria?
     
---

# Lembrando o Experimento de Bernoulli

.pull-left[
Lembre-se da [Distribuição de Bernoulli?](https://en.wikipedia.org/wiki/Bernoulli_distribution): Chamamos uma variável aleatória `$y \in \{0,1\}$` tal que

`$$\begin{align}
\Pr(y = 1) &= p \\
\Pr(y = 0) &= 1-p \\
p &\in[0,1]
\end{align}$$`

uma variável aleatória de *Bernoulli*.
]

`$$\begin{align}
\Pr(y = 1 | X = x) &= p(x) \\
\Pr(y = 0 | X = x) &= 1-p(x) \\
p(x) &\in[0,1]
\end{align}$$`

* Particularmente: *valor esperado* (ou seja, a média) de `$Y$` dado `$x$`

$$
E[y | x] = 1 \times p(x) + 0 \times (1-p(x)) = p(x)
$$

* Muitas vezes modelamos **expectativas condicionais** 😉
]
---

# O Modelo de Probabilidade Linear (MPL)

* A opção mais simples. Modele a probabilidade de resposta como

$$
\Pr(y = 1 | x) = p(x) = \beta_0 + \beta_1 x_1 + \dots + \beta_K x_K
$$
* Interpretação: *uma mudança de 1 unidade em `$x_1$`, digamos, resulta em uma mudança de `$\beta_1$` em `$p(x)$`.*

## Exemplo: Mroz (1987)

* Participação feminina no mercado de trabalho

* Como o status de `inlf` (*na força de trabalho*) depende da renda familiar da mulher solteira, sua educação, idade e número de filhos pequenos?
---

# Mroz 1987

```r
data(mroz, package = "wooldridge")
plot(factor(inlf) ~ age,
  data = mroz,
  ylevels = 2:1,
  ylab = "in labor force?"
)
```

---

# Rodando o MPL

```r
LPM <- lm(inlf ~ nwifeinc + educ + exper + I(exper^2) + age + I(age^2) + kidslt6,
  data = mroz
)
broom::tidy(LPM)
```

```
## # A tibble: 8 × 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  0.322     0.486       0.662  5.08e- 1
## 2 nwifeinc    -0.00343   0.00145    -2.36   1.86e- 2
## 3 educ         0.0375    0.00735     5.10   4.33e- 7
## 4 exper        0.0383    0.00577     6.63   6.44e-11
## 5 I(exper^2)  -0.000565  0.000189   -2.98   2.96e- 3
## 6 age         -0.00112   0.0225     -0.0497 9.60e- 1
## 7 I(age^2)    -0.000182  0.000258   -0.706  4.80e- 1
## 8 kidslt6     -0.260     0.0341     -7.64   6.72e-14
```
]

* Apenas `inlf` assume somente dois valores, 0 ou 1.

* Resultados: se a renda da mulher solteira aumenta em 10 (ou seja, 10.000 USD), `$p(x)$` cai em 0,034 (isso é um efeito pequeno!),

* uma criança pequena adicional reduziria a probabilidade de trabalho em 0,26 (isso é grande).

* Até agora, tudo simples.️]

---

# MPL: Prevendo probabilidades negativas?!

.pull-left[
<img src="02-probit_pt_files/figure-html/unnamed-chunk-3-1.svg" style="display: block; margin: auto;" />
]

.pull-right[
<br>
<br>
* As previsões do MPL de `$p(x)$` não são garantidas no intervalo unitário `$[0,1]$`.

* Lembre-se: `$e \sim D\left(0,\sigma^2\right)$`

* aqui, algumas probabilidades menores que zero!

* Particularmente irritante se você quiser *previsões*: O que é probabilidade -0,3? 🤔]

---

# MPL no modelo saturado: sem problemas!

```r
mroz %<>%
  # classify age into 3 and huswage into 2 classes
  mutate(
    age_fct = cut(age, breaks = 3, labels = FALSE),
    huswage_fct = cut(huswage, breaks = 2, labels = FALSE)
  ) %>%
  mutate(classes = paste0("age_", age_fct, "_hus_", huswage_fct))

LPM_saturated <- mroz %>%
  lm(inlf ~ classes, data = .)
broom::tidy(LPM_saturated)
```

```
## # A tibble: 6 × 5
##   term               estimate std.error statistic  p.value
##   <chr>                 <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)          0.611     0.0277    22.0   2.98e-83
## 2 classesage_1_hus_2  -0.611     0.350     -1.75  8.11e- 2
## 3 classesage_2_hus_1  -0.0257    0.0404    -0.635 5.25e- 1
## 4 classesage_2_hus_2  -0.277     0.203     -1.37  1.72e- 1
## 5 classesage_3_hus_1  -0.149     0.0494    -3.01  2.72e- 3
## 6 classesage_3_hus_2  -0.111     0.350     -0.317 7.51e- 1
```

]

* Cada classe: `$p(x)$` *dentro daquela célula*.

]

---

# MPL no modelo saturado: sem problemas!

.left-wide[
<img src="02-probit_pt_files/figure-html/saturated-1.svg" style="display: block; margin: auto;" />
]

* Cada segmento de linha: `$p(x)$` *dentro daquela célula*.

* Por exemplo. as mulheres da faixa etária mais jovem e de menor renda do marido (classe `age_1_hus_1`) têm a maior probabilidade de trabalhar (0.611).
]

---

# Modelos de Resposta Binária Não-Lineares

Nesta classe de modelos mudamos a forma como modelamos a probabilidade de resposta `$p(x)$`. Em vez da estrutura linear simples de cima, escrevemos

$$
\Pr(y = 1 | x) = p(x) = G \left(\beta_0 + \beta_1 x_1 + \dots + \beta_K x_K \right)
$$
--

* ***quase*** idêntico ao MPL!

* exceto que o *índice linear* `$\beta_0 + \beta_1 x_1 + \dots + \beta_K x_K$` agora está dentro da ***função de ligação*** `$G(\cdot)$` (i.e. _link function_).

* Propriedade principal de `$G$`: transforma qualquer `$z\in \mathbb{R}$` em um número no intervalo `$(0,1)$`.

* Isso resolve nosso problema de previsões fora do intervalo `$\left[0, 1\right]$` para probabilidades.
---

# Modelos de Resposta Binária Não-Lineares

## `$G$`: *probit* e *logit*

.left-wide[
<img src="02-probit_pt_files/figure-html/cdfs-1.svg" style="display: block; margin: auto;" />
]

<br>
Para **probit** e **logit**:

1. qualquer valor `$x$` resulta em um valor `$p(x)$` entre 0 e 1.

1. estritamente crescentes.

1. Logit tem *caudas mais longas*.

]

---

# Modelos de Resposta Binária Não-Lineares

## `$G$`: *probit* e *logit*

.left-wide[
<img src="02-probit_pt_files/figure-html/cdfs2-1.svg" style="display: block; margin: auto;" />
]

<br>
**Probit**:

* `$G(z)=\Phi(z)$` Distribuição Normal

**Logit**:

* `$G(z)=\Lambda(z)=\frac{1}{1+\exp(-z)}$` Distribuição Logística
]

---
# Modelos de Resposta Binária Não-Lineares

## Rodando probit e logit no `R`: a função `glm`

* Usamos a função `glm` para rodar um **modelo linear generalizado**

* Isso *generaliza* nosso modelo linear padrão. Temos que especificar uma `família` e um `link`:

```r
probit <- glm(inlf ~ age,
  data = mroz,
  family = binomial(link = "probit")
)

logit <- glm(inlf ~ age,
  data = mroz,
  family = binomial(link = "logit")
)
```

---

# Interpretação

```r
modelsummary::modelsummary(list("probit" = probit, "logit" = logit))
```

<table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:center;"> probit </th>
   <th style="text-align:center;"> logit </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:center;"> 0.707 </td>
   <td style="text-align:center;"> 1.136 </td>
  </tr>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:center;"> (0.248) </td>
   <td style="text-align:center;"> (0.398) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> age </td>
   <td style="text-align:center;"> -0.013 </td>
   <td style="text-align:center;"> -0.020 </td>
  </tr>
  <tr>
   <td style="text-align:left;box-shadow: 0px 1.5px">  </td>
   <td style="text-align:center;box-shadow: 0px 1.5px"> (0.006) </td>
   <td style="text-align:center;box-shadow: 0px 1.5px"> (0.009) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Num.Obs. </td>
   <td style="text-align:center;"> 753 </td>
   <td style="text-align:center;"> 753 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AIC </td>
   <td style="text-align:center;"> 1028.9 </td>
   <td style="text-align:center;"> 1028.9 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> BIC </td>
   <td style="text-align:center;"> 1038.1 </td>
   <td style="text-align:center;"> 1038.1 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Log.Lik. </td>
   <td style="text-align:center;"> -512.442 </td>
   <td style="text-align:center;"> -512.431 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> F </td>
   <td style="text-align:center;"> 4.828 </td>
   <td style="text-align:center;"> 4.858 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> RMSE </td>
   <td style="text-align:center;"> 0.49 </td>
   <td style="text-align:center;"> 0.49 </td>
  </tr>
</tbody>
</table>
]

* logit: -0.02 para logit,

* impacto da idade na probabilidade de trabalhar é **negativo**

* No entanto, **quão** negativo? Não podemos dizer!]

---

# Interpretação

$$
\Pr(y = 1 | \text{age})= G \left(x \beta\right) = G \left(\beta_0 + \beta_1 \text{age} \right) 
$$
e o *efeito marginal* de `idade` na probabilidade de resposta é

`$$\frac{\partial{\Pr(y = 1 | \text{age})}}{ \partial{\text{age}}} = g \left(\beta_0 + \beta_1 \text{age} \right) \beta_1$$`

* `$g$` é definida como `$g(z) = \frac{dG}{dz}(z)$` - a primeira derivada de `$G$` (sendo `$G$` uma distribuição, `$g$` é a função densidade).

* dado que `$G$` é não-linear, isso significa que `$g$` não será constante. Você pode experimentar isso usando este [aplicativo aqui](https://floswald.shinyapps.io/marginal_effects_of_logit_probit/):

---

# Interpretação

***Não há um único efeito marginal*** nesses modelos, pois isso depende de *onde avaliamos* a expressão anterior. Na prática, existem duas abordagens comuns:

1. reporte o **efeito parcial na média** (PEA): `$$PEA(X_j)=g(\bar{x} \beta) \beta_j$$`

1. relate o **efeito parcial médio** (APE): `$$APE(X_j)=\frac{1}{n} \sum_{i=1}^N g(x_i \beta) \beta_j$$`

Felizmente, existem pacotes disponíveis que nos ajudam a calcular esses efeitos marginais com bastante facilidade. Um deles se chama [`marginaleffects`](https://vincentarelbundock.github.io/marginaleffects/index.html), e o usamos da seguinte forma:

---

# Interpretação

```r
f <- "inlf ~ age + kidslt6 + nwifeinc" # setup a formula
glms <- list()
glms$probit <- glm(
  formula = f,
  data = mroz,
  family = binomial(link = "probit")
)
glms$logit <- glm(
  formula = f,
  data = mroz,
  family = binomial(link = "logit")
)
# now the marginal effects versions
mfx_mean <- lapply(glms, marginaleffects, newdata = "mean")
mfx_avg <- lapply(glms, marginaleffects)
```

---

# Interpretação

<table style='NAborder-bottom: 0; width: auto !important; margin-left: auto; margin-right: auto; font-family: "Arial Narrow", "Source Sans Pro", sans-serif; margin-left: auto; margin-right: auto;' class="table lightable-classic">
 <thead>
<tr>
<th style="empty-cells: hide;" colspan="1"></th>
<th style="padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #111111; margin-bottom: -1px; ">PEA</div></th>
<th style="padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #111111; margin-bottom: -1px; ">APE</div></th>
</tr>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:center;"> probit </th>
   <th style="text-align:center;"> logit </th>
   <th style="text-align:center;"> probit  </th>
   <th style="text-align:center;"> logit  </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> age </td>
   <td style="text-align:center;"> -0.013*** </td>
   <td style="text-align:center;"> -0.013*** </td>
   <td style="text-align:center;"> -0.013*** </td>
   <td style="text-align:center;"> -0.013*** </td>
  </tr>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:center;"> (0.002) </td>
   <td style="text-align:center;"> (0.002) </td>
   <td style="text-align:center;"> (0.002) </td>
   <td style="text-align:center;"> (0.002) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> kidslt6 </td>
   <td style="text-align:center;"> -0.300*** </td>
   <td style="text-align:center;"> -0.303*** </td>
   <td style="text-align:center;"> -0.290*** </td>
   <td style="text-align:center;"> -0.292*** </td>
  </tr>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:center;"> (0.039) </td>
   <td style="text-align:center;"> (0.040) </td>
   <td style="text-align:center;"> (0.036) </td>
   <td style="text-align:center;"> (0.036) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> nwifeinc </td>
   <td style="text-align:center;"> -0.004** </td>
   <td style="text-align:center;"> -0.004** </td>
   <td style="text-align:center;"> -0.004** </td>
   <td style="text-align:center;"> -0.004** </td>
  </tr>
  <tr>
   <td style="text-align:left;box-shadow: 0px 1.5px">  </td>
   <td style="text-align:center;box-shadow: 0px 1.5px"> (0.002) </td>
   <td style="text-align:center;box-shadow: 0px 1.5px"> (0.002) </td>
   <td style="text-align:center;box-shadow: 0px 1.5px"> (0.002) </td>
   <td style="text-align:center;box-shadow: 0px 1.5px"> (0.002) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Num.Obs. </td>
   <td style="text-align:center;"> 753 </td>
   <td style="text-align:center;"> 753 </td>
   <td style="text-align:center;"> 753 </td>
   <td style="text-align:center;"> 753 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Log.Lik. </td>
   <td style="text-align:center;"> -478.395 </td>
   <td style="text-align:center;"> -478.377 </td>
   <td style="text-align:center;"> -478.395 </td>
   <td style="text-align:center;"> -478.377 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> F </td>
   <td style="text-align:center;"> 21.784 </td>
   <td style="text-align:center;"> 20.280 </td>
   <td style="text-align:center;"> 21.784 </td>
   <td style="text-align:center;"> 20.280 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> RMSE </td>
   <td style="text-align:center;"> 0.47 </td>
   <td style="text-align:center;"> 0.47 </td>
   <td style="text-align:center;"> 0.47 </td>
   <td style="text-align:center;"> 0.47 </td>
  </tr>
</tbody>
<tfoot><tr><td style="padding: 0; " colspan="100%">
<sup></sup> + p &lt; 0.1, * p &lt; 0.05, ** p &lt; 0.01, *** p &lt; 0.001</td></tr></tfoot>
</table>

---
layout: false
class: separator, middle

# Qualidade do Ajuste em Modelos Binários

---
layout: true

---
# Qualidade do Ajuste em Modelos Binários

* Não existe `$R^2$` universalmente aceito para modelos binários.

* Podemos pensar em um *pseudo* `$R^2$` que compara nosso modelo com outro sem regressores:

```r
glms$probit0 <- update(glms$probit, formula = . ~ 1)
1 - as.vector(logLik(glms$probit) / logLik(glms$probit0))
```

```
## [1] 0.07084972
```

* Mas isso não é super informativo (ao contrário do `$R^2$`). As alterações no valor da log-verossimilhança são altamente não lineares.

---
# Qualidade do Ajuste em Modelos Binários

* Vamos verificar **precisão** - qual é a proporção prevista corretamente! `round(fitted(x))` atribui `1` se a prob `$> 0.5$` previsto.

```r
prop_tbl <- prop.table(table(
  true = mroz$inlf,
  pred = round(fitted(glms$probit))
)) |>
  as.data.table() |>
  dcast(true ~ pred, value.var = "N")
kbl(prop_tbl,
  digits = 4,
  caption = "Matriz de Confusão"
) |>
  add_header_above(c(
    " " = 1,
    "Pred" = 2
  ))
```

<table>
<caption>Matriz de Confusão</caption>
 <thead>
<tr>
<th style="empty-cells: hide;border-bottom:hidden;" colspan="1"></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Pred</div></th>
</tr>
  <tr>
   <th style="text-align:left;"> true </th>
   <th style="text-align:right;"> 0 </th>
   <th style="text-align:right;"> 1 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 0 </td>
   <td style="text-align:right;"> 0.1700 </td>
   <td style="text-align:right;"> 0.2616 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1 </td>
   <td style="text-align:right;"> 0.1222 </td>
   <td style="text-align:right;"> 0.4462 </td>
  </tr>
</tbody>
</table>

---

# Curvas ROC (Receiver Operating Characteristics)

## [Nomenclatura](https://en.wikipedia.org/wiki/Receiver_operating_characteristic)

**Taxa de Falsos Positivos** (FPR): `$FP/N=$` 0.6061538

**Acurácia**: `$(TP + TN)/ (P+N)=$` 0.6162019
]

.panelset[
.panel[.panel-name[TPR]
<table>
 <thead>
<tr>
<th style="empty-cells: hide;border-bottom:hidden;" colspan="1"></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Pred</div></th>
</tr>
  <tr>
   <th style="text-align:left;"> true </th>
   <th style="text-align:right;"> 0 </th>
   <th style="text-align:right;"> 1 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 0 </td>
   <td style="text-align:right;color: black !important;background-color: white !important;"> 0.1700 </td>
   <td style="text-align:right;color: black !important;background-color: white !important;"> 0.2616 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1 </td>
   <td style="text-align:right;color: white !important;background-color: darkgreen !important;"> 0.1222 </td>
   <td style="text-align:right;color: white !important;background-color: darkgreen !important;"> 0.4462 </td>
  </tr>
</tbody>
</table>

<br>
.center[Relacionado a **Poder**]
]

.panel[.panel-name[FPR]
<table>
 <thead>
<tr>
<th style="empty-cells: hide;border-bottom:hidden;" colspan="1"></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Pred</div></th>
</tr>
  <tr>
   <th style="text-align:left;"> true </th>
   <th style="text-align:right;"> 0 </th>
   <th style="text-align:right;"> 1 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 0 </td>
   <td style="text-align:right;color: white !important;background-color: darkgreen !important;"> 0.1700 </td>
   <td style="text-align:right;color: white !important;background-color: darkgreen !important;"> 0.2616 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1 </td>
   <td style="text-align:right;color: black !important;background-color: white !important;"> 0.1222 </td>
   <td style="text-align:right;color: black !important;background-color: white !important;"> 0.4462 </td>
  </tr>
</tbody>
</table>

<br>
.center[Relacionado a **Erro Tipo I**]
]

.panel[.panel-name[ACC]
<table>
 <thead>
<tr>
<th style="empty-cells: hide;border-bottom:hidden;" colspan="1"></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Pred</div></th>
</tr>
  <tr>
   <th style="text-align:left;"> true </th>
   <th style="text-align:right;"> 0 </th>
   <th style="text-align:right;"> 1 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 0 </td>
   <td style="text-align:right;color: white !important;background-color: darkgreen !important;"> 0.1700 </td>
   <td style="text-align:right;color: black !important;background-color: white !important;"> 0.2616 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1 </td>
   <td style="text-align:right;color: black !important;background-color: white !important;"> 0.1222 </td>
   <td style="text-align:right;color: white !important;background-color: darkgreen !important;"> 0.4462 </td>
  </tr>
</tbody>
</table>
]

]
]

---

# Curvas ROC

* O corte de 0,5 é arbitrário. E se todas as probabilidades previstas forem `$> 0,5$`, mas nos dados houver cerca de 50% de zeros?

* Vamos escolher um *corte arbitrário* `$c \in (0,1)$` e verificar a precisão de cada valor. Isso dá uma visão melhor.

---
# Curvas ROC

* Podemos confrontar a **taxa de verdadeiros positivos** (TPR) com a **taxa de falsos positivos** (FPR), para cada valor de corte `$c\in[0, 1]$`.

1. TPR: número de mulheres corretamente previstas para trabalhar dividido pelo número de mulheres que trabalham.
     
     2. FPR: número de mulheres incorretamente previstas para trabalhar dividido pelo número de mulheres não trabalhadoras.
    
--

* Plotar ***FPR vs TPR*** para cada `$c$` define a curva **ROC**.

* Um bom modelo tem uma curva ROC no canto superior esquerdo: FPR = 0, TPR = 1.

---

# Curvas ROC

```r
pred_prob <- prediction(fitted(glms$probit), mroz$inlf)
par(mfrow = c(1, 2), mar = lowtop)
acc <- ROCR::performance(pred_prob, "acc")
max_acc <- which.max(acc@y.values[[1]])
plot(acc)
points(acc@x.values[[1]][max_acc], acc@y.values[[1]][max_acc], col = "darkgreen")
roc <- ROCR::performance(pred_prob, "tpr", "fpr")
plot(roc)
points(roc@x.values[[1]][max_acc], roc@y.values[[1]][max_acc], col = "darkgreen")
abline(0, 1, lty = 2, col = "red")
```

<img src="02-probit_pt_files/figure-html/unnamed-chunk-10-1.svg" style="display: block; margin: auto;" />
]

* ROC sempre acima da linha de 45 graus. Melhor do que atribuição aleatória (jogar uma moeda)!

]

---
# Estimação de Máxima Verossimilhança

* Ao invés de usar média e variância condicionais, usaremos a distribuição condicional completa

* Amostra *iid* `$\lbrace y_i, x_i \rbrace_{i=1}^N$`. Queremos estimar a distribuição `$f(y|x)$`

* **Hipótese:** esta distribuição é conhecida, a não ser por um número finito de parâmetros fixos.
    + Impomos um modelo paramétrico para a densidade condicional
    
---

# Estimação de Máxima Verossimilhança

## Exemplo Probit

* Não observamos `$y_i^*$`, apenas `$y_i$`
]

.pull-right[
`$$\begin{equation*}
y_i = \begin{cases}
1 \text{ se } y_i^* > 0\\ 
0 \text{ se } y_i^* \leq 0
\end{cases}
\end{equation*}$$`
]

---

# Estimação de Máxima Verossimilhança

## Exemplo Probit

* Dadas as especificações acima:

.pull-left[
`$$\begin{align*}
P(y_i=1|x_i)&=P(y_i* > 0|x_i)\\
&=P(\varepsilon_i > -x_i\theta |x_i)\\
&=1-\Phi(-x_i\theta)\\
&=\Phi(x_i\theta)
\end{align*}$$`
]
--
.pull-right[
`$$\begin{align*}
P(y_i=0|x_i)&=P(y_i* \leq 0|x_i)\\
&=P(\varepsilon_i \leq -x_i\theta |x_i)\\
&=\Phi(-x_i\theta)\\
&=1-\Phi(x_i\theta)
\end{align*}$$`
]

---
# Estimação de Máxima Verossimilhança

## Exemplo Probit

* A função de probabilidade condicional pode ser escrita como:

`$$p(y_i|x_i)=\Phi(x_i\theta)^{y_i}\cdot [1-\Phi(x_i\theta)]^{1-y_i}$$`
--

* Esta é a probabilidade de ocorrência de um ponto `$(x_i, y_i)$`. Para o conjunto de dados observados temos a ***função de verossimilhança***

`$$\ell(\theta)=\Pi_{i=1}^N \Phi(x_i\theta)^{y_i}\cdot [1-\Phi(x_i\theta)]^{1-y_i}$$`

* Máxima Verossimilhança consiste em ***maximizar a função de verossimilhança*** em relação aos parâmetros, `$\theta$`.

* Podemos maximizar a ***log-verossimilhança*** por praticidade. `$\mathcal{L}(\theta)=\log(\ell(\theta))=\sum_{i=1}^N y_i\log(\Phi(x_i\theta))+(1-y_i)\log(1-\Phi(x_i\theta))$`

---
# Estimação de Máxima Verossimilhança

## Exemplo Probit

`$$\hat\theta = \arg \max_{\theta} \mathcal{L}(\theta)$$`

* O vetor de ***score*** é dado pela derivada da log-verossimilhança em relação a cada um dos parâmetros

`$$s(\theta)=\nabla_\theta \mathcal{L}(\theta)=\left[\frac{\partial\mathcal{L}(\theta)}{\partial\theta_1}, \ldots, \frac{\partial\mathcal{L}(\theta)}{\partial\theta_k}\right]^{\prime}$$`
* A Hessiana é a matriz de segundas derivadas

`$$H(\theta)=\nabla^2_\theta \mathcal{L}(\theta)$$`

---
# Estimação de Máxima Verossimilhança

## Exemplo Probit

Exercício: Calcule o _score_ do modelo probit.

---
# Estimação de Máxima Verossimilhança

## Exemplo Probit - Solução

`$\mathcal{L}(\theta)=\sum_{i=1}^N y_i\log(\Phi(x_i\theta))+(1-y_i)\log(1-\Phi(x_i\theta))$`

`$$\begin{align}
s(\theta)&=\frac{d\mathcal{L}(\theta)}{d\theta}=\sum_{i=1}^N y_i x_i\frac{\phi(x_i \theta)}{\Phi(x_i \theta)} - (1-y_i)x_i\frac{\phi(x_i \theta)}{1-\Phi(x_i \theta)}\\
&=\sum_{i=1}^N \frac{[1-\Phi(x_i \theta)]y_i x_i\phi(x_i \theta) - \Phi(x_i \theta)(1-y_i)x_i\phi(x_i \theta)}{\Phi(x_i \theta) [1-\Phi(x_i \theta)]}\\
&=\sum_{i=1}^N \frac{x_i\phi(x_i \theta)}{\Phi(x_i \theta) [1-\Phi(x_i \theta)]}[y_i-\Phi(x_i \theta)]
\end{align}$$`

---
# Estimação de Máxima Verossimilhança

## Exemplo Exponencial

- Suponha que tenhamos um conjunto simples de dados `$\{x_i\}_{i=1}^N$` `$iid$` e queremos ajustar uma distribuição exponencial

- `$x_i\sim Exp(\lambda)$`. `$f(x; \lambda)=\frac{1}{\lambda}e^{-x_i/\lambda}$`. `$x_i\geq 0$`.

- Qual é a função de log-verossimilhança?
- Encontre o estimador de máxima verossimilhança `$\hat\lambda_{mle}$`

---

# Estimação de Máxima Verossimilhança

## Exemplo Exponencial

- Função de verossimilhança: `$L(\lambda)=\Pi_{i=1}^N f(x_i; \lambda)=\Pi_{i=1}^N \frac{1}{\lambda}e^{-x_i/\lambda}$`

- Log-verossimilhança: `$\mathcal{L}(\lambda)=\sum_{i=1}^N \log \frac{1}{\lambda}e^{-x_i/\lambda}$`

- Ao final encontramos ...

.pull-right[
.center[
<img src="https://media.giphy.com/media/Jlirxn597uSeW7kPAQ/giphy.gif", height="250" />
]
]

---
# Método dos Momentos Generalizado

## Método dos Momentos (MM)

* Podemos utilizar momentos populacionais para identificar parâmetros
    + `$E[u_i]=0$` é um exemplo
--
* De fato, ***tanto MQO quanto VI*** podem ser entendidos como estimadores MM!

* Quais são os momentos utilizados no MQO?

* `$E[u]=0$` e `$E[xu]=0$`, que são as equações normais do MQO

* O MM resolve o análogo amostral destas equações

---

# Método dos Momentos Generalizado

## Método dos Momentos (MM)

* Vamos considerar um _setup_ mais geral:
* `$\theta_0\in\mathbb{R}^p$` é um vetor dos parâmetros da população
* `$m(\theta_0)\equiv E[g(Z_i, \theta_0)]=0$`, são as equações de momento
* `$g(Z_i, \theta)$` é uma função vetorial `$r\times 1$`
* `$Z_i\equiv (X_i, Y_i)$` é uma observação do conjunto de dados 
* `$\theta_0$` é a única solução para a equação de momento
* MM vai resolver o análogo amostral:

`$$\hat{m}(\hat\theta)=\frac{1}{N}\sum_{i=1}^N g(Z_i, \hat\theta)=0$$`
]

.right-thin[

![](https://media.giphy.com/media/YmnRIWXjSkct9ZZtNZ/giphy.gif)
]
---

# Método dos Momentos Generalizado

## Clarificando

.left-wide[
* Considere o MQO. Podemos escrever duas equação de momento que identificam os parâmetros do modelo

`$$\hat{m}(\hat\beta)=\begin{cases}
\frac{1}{N}\sum_{i=1}^N(y_i-\hat\beta_0-\hat\beta_1 x_i)=0\\
\frac{1}{N}\sum_{i=1}^Nx_i(y_i-\hat\beta_0-\hat\beta_1 x_i)=0
\end{cases}$$`

* O MM escreve momentos populacionais em termos dos parâmetros e resolve o sistema de equações, simples!
]
--

---
# Método dos Momentos Generalizado

## GMM

* Podemos generalizar o método dos momentos. Daí o nome ***generalized method of moments*** - **GMM**

* Minimizar a norma `${\displaystyle \|{\hat {m}}(\theta )\|_{W}^{2}={\hat {m}}(\theta )^{\mathsf {T}}\,W{\hat {m}}(\theta)}$`

* Onde `$W$` é uma matriz positiva definida

* `${\displaystyle {\hat {\theta }}=\operatorname {arg} \min _{\theta \in \Theta }{\bigg (}{\frac {1}{N}}\sum _{i=1}^{N}g(Z_{i},\theta ){\bigg )}^{\mathsf {T}}{W}{\bigg (}{\frac {1}{N}}\sum_{i=1}^{N}g(Z_{i},\theta ){\bigg )}}$`

--
* Especialmente útil quando temos mais equações de momento do que parâmetros

---
# Método dos Momentos Generalizado

## Aplicações de GMM

* MQO, VI e MLE podem ser obtidos a partir de um GMM

* Aplicações diversas em
    + _Cross-Section_: modelos não-lineares com variáveis endógenas 
    + Séries Temporais: correlação serial e heterocedasticidade
    + Painel: extensões do modelo de efeitos fixos

* Veja mais em [Wooldridge (2001) JEP](https://pubs.aeaweb.org/doi/pdfplus/10.1257/jep.15.4.87)

---

# Leitura Recomendada

* WOOLDRIDGE, Jeffrey M. Introdução à econometria: uma abordagem moderna. São Paulo: Cengage Learning, 2016. Tradução da 4ª edição norte-americana por José Antonio Ferreira. Capítulo 17 Modelo com Variáveis Dependentes Limitadas e Correções da Seleção Amostral.

* GUJARATI, Damodar N.; PORTER, Dawn C. Econometria básica. Porto Alegre: Amgh Editora, 2011. - 5. ed. Capítulo 15 Modelos de regressão de resposta qualitativa

* WOOLDRIDGE, Jeffrey M. Econometric Analysis of Cross Section and Panel Data. MIT press, 2010. Second Edition. Chapter 15 Binary Response Models

---

class: title-slide-final, middle
background-image: url(../../img/logo/UdescEsag.jpeg)
background-size: 350px
background-position: 9% 19%

# ATÉ A PRÓXIMA AULA!

.footnote[
[1]: Este slides foram baseados nas aulas de econometria da [SciencesPo Department of Economics](https://github.com/ScPoEcon/ScPoEconometrics-Slides)
]

|                                                                                                            |                                   |
| :--------------------------------------------------------------------------------------------------------- | :-------------------------------- |
| <a href="https://github.com/rfbressan/econometria3_slides">.ScPored[<i class="fa fa-link fa-fw"></i>] | Slides |
| <a href="http://github.com/rfbressan">.ScPored[<i class="fa fa-github fa-fw"></i>]                          | @rfbressan                      |
| <a href="https://raw.githack.com/rfbressan/econometria3_slides/master/lectures/02-probit/lista_I_pt.html">.ScPored[<i class="fa fa-list fa-fw"></i>] | Lista de Exercícios I |