class: inverse, center, middle

```{r xaringan-themer, include=FALSE, warning=FALSE}
library(tidyverse)
library(xaringanthemer)
style_duo_accent(primary_color ="#001F3F", secondary_color = "#225522",
                 background_color = "#f8f7f3",
                 header_font_google = google_font("Archivo"),
                 text_font_google   = google_font("Inconsolata"), 
                 link_color= "#800000"

)
```

#Regresión Logística
##Interpretación

---
##Configuración

--

<br>

- Tenemos $y_{1}, \dots y_{n}$ son $n$ variables independientes con distribución $\text{Bernoulli}(p_{i})$


--

- $\mathbb{E}(y_{i} \mid x_{i1}, \dots,x_{ik}) = \mathbb{P}(y_{i}=1 \mid x_{i1}, \dots,x_{ik}) = p_{i}$ 

<br>
--

donde, 


$$\ln \frac{p_{i}}{1 - p_{i}} = \beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}$$

---
##Configuración

<br>

- Tenemos $y_{1}, \dots y_{n}$ son $n$ variables independientes con distribución $\text{Bernoulli}(p_{i})$


- $\mathbb{E}(y_{i} \mid x_{i1}, \dots,x_{ik}) = \mathbb{P}(y_{i}=1 \mid x_{i1}, \dots,x_{ik}) = p_{i}$ 

<br>

donde, 

$$\underbrace{\ln \frac{p_{i}}{1 - p_{i}}}_{\text{Link logit}(p_{i})} = \beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}$$
---
##Configuración


<br>

- Tenemos $y_{1}, \dots y_{n}$ son $n$ variables independientes con distribución $\text{Bernoulli}(p_{i})$



- $\mathbb{E}(y_{i} \mid x_{i1}, \dots,x_{ik}) = \mathbb{P}(y_{i}=1 \mid x_{i1}, \dots,x_{ik}) = p_{i}$ 

<br>

donde, 

$$\ln \frac{p_{i}}{1 - p_{i}} = \underbrace{\beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}}_{\text{Predictor lineal  } \eta_{i}}$$
<br>
--

- Interpretación: ¿cuál es el significado de $\beta_{0},\beta_{1}, \dots, \beta_{k}$?

---
## Un ejemplo empírico

```{r,  include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
# load data on extra-marital affairs from package "Ecdat"
library("Ecdat")
library("viridis")
library("tidyverse")
library("modelr")
library("cowplot")
library("margins")
library("marginaleffects")


theme_set(theme_cowplot())

data(Fair)
affairsdata <- Fair %>% as_tibble()

# create a binary variable indicating wether persons has ever had an affair
affairsdata <- affairsdata %>% 
  mutate(everaffair = case_when(nbaffairs == 0 ~ "Never", nbaffairs > 0 ~ "At least once") ) %>%
  # map into 0/1 code
  mutate(everaffair_d = case_when(nbaffairs == 0 ~ 0, nbaffairs > 0 ~ 1))
```

.pull-left[
Continuando con los datos de infidelidad, ajustaremos el siguiente modelo:

<br>
$$\underbrace{\ln \frac{p_{i}}{ 1 - p_{i}}}_{\text{logit}(p_{i})}    = \beta_{0} + \beta_{1}\text{male}_{i} + \beta_{2}\text{ym}_{i}$$
<br>

donde:

- $p_{i} =\mathbb{P}(\text{everaffair}_{i}=1)$

- $\text{logit}(p_{i})$ es el .bold[log odds] de tener un affair

- $p_{i}$ y $\text{logit}(p_{i})$ son una función género (male) y años de matrimonio (ym).

]

--
.pull-right[
```{r}
logit_affairs_sex_ym <- 
  glm(everaffair_d ~ factor(sex) + ym, 
      family=binomial(link="logit"), 
      data=affairsdata)
summary(logit_affairs_sex_ym)
```
]

---
class: inverse, center, middle

## Efectos marginales sobre el logit 

---
## Efectos marginales sobre el logit 

<br>
<br>

.center[ 
```{r, echo=FALSE, fig.width=6, fig.height=5, warning=FALSE, message=FALSE}
fakedata <- tibble(x=seq(-6,6,by=0.01), logit=1 + 0.7*x, odds=exp(logit), probabilities=1/(1+exp(-logit)))

fakedata %>%  ggplot(aes(x=x, y=logit)) +  
  geom_line(size=2, alpha = 1) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
  guides(fill="none", color="none") +
  theme(axis.text.y = element_text(size = 22), axis.text.x = element_text(size = 22),
  axis.title.y = element_text(size = 24), axis.title.x = element_text(size = 24), 
  legend.text = element_text(size = 18), legend.position="bottom") +
  labs(x="x", y="logit(x)") 
```
]


---
## Efectos marginales sobre el logit 

Dado el siguiente modelo de regresión logística: 

$$\text{logit}(p_{i}) = \ln \frac{p_{i}}{1 - p_{i}} = \beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}$$
--

- El intercepto $\beta_{0}$ corresponde al logit de la probabilidad de éxito cuando $x_{1} = \dots = x_{k} = 0$

--

- El efecto marginal de $x_{k}$ sobre el logit de la probabilidad de éxito está dado por:


.pull-left[
.content-box-blue[
$$\frac{\partial\text{logit}(p_{i})}{\partial x_{k}} = \beta_{k}$$
]
]
.pull-right[
.content-box-yellow[
"Un cambio infinitesimal en $x_{k}$ $\big(\partial x_{k}\big)$ se traduce en un cambio en $\partial x_{k} \beta_{k}$ unidades en el logit de p"
] 
]

--

¿Por qué?
--
 Asumamos que sólo la variable $x_{k}$ cambia en una unidad, de $x_{k}=c$ a $x_{k}=c+1$. Entonces:

--

\begin{align}
  \frac{\Delta\text{logit}(p_{i}) }{\Delta x_{k}} &= (\beta_{0} + \beta_{1}x_{1} +  \dots + \beta_{k}(c+1)) - (\beta_{0} + \beta_{1}x_{1}  \dots + \beta_{k}c)  \\ 
  &=  \beta_{k}(c+1) - \beta_{k}c = \beta_{k}
\end{align}


---
## Efectos marginales sobre el logit 

En nuestro ejemplo: $\ln \frac{p_{i}}{1 -p_{i}} = \beta_{0} + \beta_{1}*\text{male}_{i} + \beta_{2}*\text{ym}_{i}$

<br>

.pull-left[
```{r, echo=FALSE}
summary(logit_affairs_sex_ym)$coefficients[,c(1,2)]
```
]
.pull-right[
```{r, echo=FALSE, fig.width=6, fig.height=5, warning=FALSE, message=FALSE}

# plot the result

grid <- affairsdata  %>% data_grid(sex,ym=seq(0,90),.model=logit_affairs_sex_ym)
predictions <- cbind(grid,logit_hat = predict(logit_affairs_sex_ym, newdata = grid)) 

predictions %>% ggplot(aes(x=ym, y=logit_hat, group=sex, colour=sex)) +  
  geom_line(size=2, alpha = 1) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
guides(fill="none", color="none") +
  theme(axis.text.y = element_text(size = 22), axis.text.x = element_text(size = 22),
  axis.title.y = element_text(size = 24), axis.title.x = element_text(size = 24), 
  legend.text = element_text(size = 18), legend.position="bottom") +
  labs(x="years of marriage", y="logit(p)") +
  annotate('text', x = 37, y = -1, label = "beta[ym]==0.0584", parse = TRUE, size=8) 
```
]

---
class: inverse, center, middle

## Efectos multiplicativos sobre las odds de éxito 

---
## Efectos multiplicativos sobre las odds 

Dado el siguiente modelo de regresión logística: 

$$\text{logit}(p_{i}) = \ln \frac{p_{i}}{1 - p_{i}} = \beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}$$

<br>
--

exponenciando a ambos lados obtenemos las .bold[odds de éxito]:

$$\frac{p_{i}}{1 - p_{i}} =  e^{\beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}}$$

--

equivalentemente

.content-box-blue[
$$\frac{p_{i}}{1 - p_{i}} =  e^{\beta_{0}} \cdot e^{\beta_{1} x_{i1}}  \dots e^{\beta_{k} x_{ik}}$$
]

---
## Efectos multiplicativos sobre las odds 

$$\frac{p_{i}}{1 - p_{i}} =  e^{\beta_{0}} \cdot e^{\beta_{1} x_{i1}}  \dots e^{\beta_{k} x_{ik}}$$
--

Examinando esta ecuación observamos que cuando $x_1 = \dots = x_k = 0$, 


<br>


$$\frac{p_{i}}{1 - p_{i}}=  e^{\beta_{0}}$$
Es decir, cuando $x_1 = \dots = x_k = 0$ las odds de éxito están dadas por $e^{\beta_{0}}$.


<br>
--

Generalizando,  observamos que cuando $x_1 = c_1  \dots = x_k = c_k$,  las odds de éxito están dadas por

<br>

$$\frac{p_{i}}{1 - p_{i}} =  e^{\beta_{0}} \cdot e^{\beta_{1} c_{1}}  \dots e^{\beta_{k} c_{k}}$$


---
## Efectos multiplicativos sobre las odds: odds ratios

Considera la situación en que $i$ y $j$ son dos observaciones con $x_{k}=c$ y $x_{k}=c+1$, respectivamente. El resto de las covariables toman valores idénticos. 
--
 Las odds de éxito para cada caso son:


- $p_{i}/(1 - p_{i}) = e^{\beta_{0}} \cdot e^{\beta_{1} x_{i1}}  \dots (e^{\beta_{k}})^{c}$

- $p_{j}/(1 - p_{j}) = e^{\beta_{0}} \cdot e^{\beta_{1} x_{i1}}  \dots (e^{\beta_{k}})^{c+1}$


<br>
--

El ratio de las odd de éxito entre $j$ e $i$ está dado por:

\begin{align}
\frac{p_{j}/(1 - p_{j})}{p_{i}/(1 - p_{i})} &= \frac{e^{\beta_{0}} \cdot e^{\beta_{1} x_{i1}}  \dots (e^{\beta_{k}})^{c+1}}{e^{\beta_{0}} \cdot e^{\beta_{1} x_{i1}}  \dots (e^{\beta_{k}})^{c}} \\ \\
&= e^{\beta_{k}}
\end{align}

En otras palabras, manteniendo otros factores constantes, $e^{\beta_{k}}$ representa la odds ratio entre el caso con $x_{k}$ aumentado en una unidad, y el caso con $x_{k}$ en un nivel basal dado. 


---
## Efectos multiplicativos sobre las odds 


.content-box-yellow[
"Un cambio en $\Delta$ unidades de $x_{k}$ multiplica las las odds de éxito por $e^{\Delta \beta_{k}}$"
] 

<br>
.bold[Propiedades]:

--

- $e^{\beta_{k}}$ está restringido al rango $[0,\infty+)$. Es una constante que "comprime" o amplifica las odds de éxito.

--

- Si $\beta_{k} < 0  \to  (0 < e^{\beta_{k}} < 1)$. Es decir, un aumento en $x_{k}$ está asociado con una reducción (multiplicativa) de las odds de éxito.

--

- Si $\beta_{k} = 0  \to  (e^{\beta_{k}} =1)$. Es decir, un cambio en $x_{k}$ está asociado a un cambio nulo en las odds de éxito.

--

- Si $\beta_{k} > 0  \to  (e^{\beta_{k}} > 1)$. Es decir, un aumento en $x_{k}$ está asociado a aumento (multiplicativo) en de las odds de éxito.

---
## Efectos multiplicativos sobre las odds 

En nuestro ejemplo: $\ln \frac{p_{i}}{1 -p_{i}} = \beta_{0} + \beta_{1}*\text{male}_{i} + \beta_{2}*\text{ym}_{i}$, por tanto:
  
  
.pull-left[
  $\frac{p_{i}}{1 -p_{i}} = e^{\beta_{0}} \cdot  e^{\beta_{1}*\text{male}_{i}} \cdot e^{\beta_{2}*\text{ym}_{i}}$
      
```{r}
# coeffs
summary(logit_affairs_sex_ym)$coefficients[,c(1,2)]
    
# exp(coeffs)
exp(summary(logit_affairs_sex_ym)$coefficients[,c(1)])
```
    
]

--
  
.pull-right[
  ```{r, echo=FALSE, fig.height=6, warnings=FALSE}
  
  # plot the result
  
  grid <- affairsdata  %>%data_grid(sex,ym=seq(0,90),.model=logit_affairs_sex_ym)
  predictions <- cbind(grid,logit_hat = predict(logit_affairs_sex_ym, newdata = grid)) %>%
    mutate(odds = exp(logit_hat)) 
  
  predictions %>% ggplot(aes(x=ym, y=odds, group=sex, colour=sex)) +  
    geom_line(size=2, alpha = 1) +
    scale_color_viridis_d() +  scale_fill_viridis_d() +
      guides(fill="none", color="none") +
    theme(axis.text.y = element_text(size = 22), axis.text.x = element_text(size = 22),
          axis.title.y = element_text(size = 24), axis.title.x = element_text(size = 24), 
          legend.text = element_text(size = 18), legend.position="bottom") +
    labs(x="years of marriage", y="Odds of Affair") +
    annotate('text', x = 60, y = 1, label = "e^beta[ym]==1.06", parse = TRUE, size=8) +
    annotate('text', x = 58, y = 20, label = "e^beta[male]==1.21", parse = TRUE, size=8) 
  ```
]


---
class: inverse, center, middle

## Efectos marginales sobre la probabilidad de éxito 


---
## Efectos marginales sobre la probabilidad de éxito 

<br>
<br>

.center[ 
```{r, echo=FALSE, fig.width=6, fig.height=5, warning=FALSE, message=FALSE}
fakedata <- tibble(x=seq(-6,6,by=0.01), logit=1 + 0.7*x, odds=exp(logit), probabilities=1/(1+exp(-logit)))

fakedata %>%  ggplot(aes(x=x, y=probabilities)) +  
  geom_line(size=2, alpha = 1) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
  guides(fill="none", color="none") +
  theme(axis.text.y = element_text(size = 22), axis.text.x = element_text(size = 22),
  axis.title.y = element_text(size = 24), axis.title.x = element_text(size = 24), 
  legend.text = element_text(size = 18), legend.position="bottom") +
  labs(x="x", y="p(x)") 
```
]


---
## Efectos marginales sobre la probabilidad de éxito


Dado el siguiente modelo de regresión logística: 

$$\ln \frac{p_{i}}{1 - p_{i}} = \beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}$$

<br>
--
Queremos saber el .bold[efecto marginal] de los predictores sobre la .bold[probabilidad] de éxito. Formalmente
--
$$\frac{\partial p_{i}}{\partial x_{k}}$$
--
.bold[Paso a paso...]

.img-right[
![mostaza](mostaza-merlo.jpg)
]

<br>
--

- .bold[Paso 1]: reexpresar la ecuación de regresión como

<br>

$$\ln(p_{i}) - \ln(1-p_{i}) = \beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}$$


---
## Efectos marginales sobre la probabilidad de éxito


- .bold[Paso 1]: reexpresar la ecuación de regresión como


$$\ln(p_{i}) - \ln(1-p_{i}) = \beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}$$

<br>
--

- .bold[Paso 2]: usando diferenciación implícita y tomando derivada con respecto a $x_{k}$ en ambos lados, obtenemos:

$$\frac{\partial p_{i} }{\partial x_{k}} \frac{1}{p_{i}} + \frac{\partial p_{i} }{\partial x_{k}} \frac{1}{(1 - p_{i})} = \beta_{k}$$
<br>
--

- .bold[Paso 3]: re-agrupando:

$$\frac{\partial p_{i} }{\partial x_{k}} \frac{1}{p_{i}(1 - p_{i})}= \beta_{k}$$

---
## Efectos marginales sobre la probabilidad de éxito


- .bold[Paso 3]: re-agrupando:

$$\frac{\partial p_{i} }{\partial x_{k}} \frac{1}{p_{i}(1 - p_{i})}= \beta_{k}$$

<br>
--

- .bold[Paso 4]: re-ordebando los términos obtenemos el efecto marginal que buscamos:

.content-box-blue[
$$\frac{\partial p_{i}}{\partial x_{k}}= \beta_{k} \cdot p_{i}(1 - p_{i})$$
]

--

recordar que ...
  $$p_{i} = \frac{1}{1 + e^{-(\beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik})}}$$

---
## Efectos marginales sobre la probabilidad de éxito

- .bold[Paso 5]: reemplazando,

$$\frac{\partial p_{i}}{\partial x_{k}}= \beta_{k} \cdot \Bigg(\frac{1}{1 + e^{-(\beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}})}\Bigg) \Bigg(1 - \frac{1}{1 + e^{-(\beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik})}}\Bigg)$$

<br>
--

Corolario:

--

- Dado que $p_{i}(1-p_{i})$ es estrictamente positivo,
 - $\beta_{k}>0$ indica que variable $x_{k}$ tienen un efecto positivo sobre la probabilidad de éxito
 - $\beta_{k}<0$ indica que $x_{k}$ tiene un efecto negativo.

<br>
--

- El .bold[efecto marginal] de $x_{k}$ no depende sólo de $\beta_{k}$ sino varía dependiendo del valor de $x_{k}$ y del valor de todas las otras covariables.


---
## Efectos marginales sobre la probabilidad de éxito

$$\frac{\partial p_{i}}{\partial x_{k}}= \beta_{k} \cdot \Bigg(\frac{1}{1 + e^{-(\beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik}})}\Bigg) \Bigg(1 - \frac{1}{1 + e^{-(\beta_{0} + \beta_{1} x_{i1} + \dots + \beta_{k} x_{ik})}}\Bigg)$$

<br>
--

Corolario, 

--

- A medida que la $p_{i}$ se aproxima a 0, $\frac{\partial p_{i}}{\partial x_{k}}$ tiende a 0.

- A medida que la $p_{i}$ se aproxima a 1, $\frac{\partial p_{i}}{\partial x_{k}}$ tiende a 0.

- A medida que la $p_{i}$ se acerca a 0.5, $\frac{\partial p_{i}}{\partial x_{k}}$ tiende a $\beta_{k}/4$.


---
## Efectos marginales sobre la probabilidad de éxito

En nuestro ejemplo: $\ln \frac{p_{i}}{1 -p_{i}} = \beta_{0} + \beta_{1}*\text{male}_{i} + \beta_{2}*\text{ym}_{i}$, por tanto:
$$p_{i} = \frac{1}{1 + e^{-(\beta_{0} + \beta_{1}*\text{male}_{i} + \beta_{2}*\text{ym}_{i})} }$$

.pull-left[
```{r, echo=FALSE}
summary(logit_affairs_sex_ym)$coefficients[,c(1,2)]
```
]
.pull-right[
```{r, echo=FALSE, fig.width=6, fig.height=5}

# plot the result

grid <- affairsdata  %>% data_grid(sex,ym=seq(0,90),.model=logit_affairs_sex_ym)
predictions <- cbind(grid,logit_hat = predict(logit_affairs_sex_ym, newdata = grid)) %>%
                mutate(p_hat = 1/(1 + exp(-logit_hat)))

predictions %>% ggplot(aes(x=ym, y=p_hat, group=sex, colour=sex)) +  
  geom_line(size=2, alpha = 1) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
guides(fill="none", color="none") +
  theme(axis.text.y = element_text(size = 22), axis.text.x = element_text(size = 22),
  axis.title.y = element_text(size = 24), axis.title.x = element_text(size = 24), 
  legend.text = element_text(size = 18), legend.position="bottom") +
  labs(x="years of marriage", y="p") +
  annotate('text', x = 60, y = 0.3, label = "slope==0.0584*p*(1-p)", parse = TRUE, size=8) 
```
]

---
## Efectos marginales sobre la probabilidad de éxito

.pull-left[
$$p_{i} = \frac{1}{1 + e^{-(\beta_{0} + \beta_{1}*\text{male}_{i} + \beta_{2}*\text{ym}_{i})}}$$

```{r, echo=FALSE, fig.width=6, fig.height=5}

# plot the result

grid <- affairsdata  %>% data_grid(sex,ym=seq(0,90),.model=logit_affairs_sex_ym)
predictions <- cbind(grid,logit_hat = predict(logit_affairs_sex_ym, newdata = grid)) %>%
                mutate(p_hat = 1/(1 + exp(-logit_hat)))

predictions %>% ggplot(aes(x=ym, y=p_hat, group=sex, colour=sex)) +  
  geom_line(size=2, alpha = 1) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
guides(fill="none", color="none") +
  theme(axis.text.y = element_text(size = 22), axis.text.x = element_text(size = 22),
  axis.title.y = element_text(size = 24), axis.title.x = element_text(size = 24), 
  legend.text = element_text(size = 18), legend.position="bottom") +
  labs(x="years of marriage", y="p") 
```
]

--

.pull-right[
$$\frac{\partial p_{i}}{\partial \text{ym} }= \beta_{\text{ym}} \cdot p_{i}(1 - p_{i})$$

```{r, echo=FALSE, fig.width=6, fig.height=5}

beta_ym = logit_affairs_sex_ym$coefficients[3]

grid <- affairsdata  %>% data_grid(sex,ym=seq(0,90),.model=logit_affairs_sex_ym)
predictions <- cbind(grid,logit_hat = predict(logit_affairs_sex_ym, newdata = grid)) %>%
                mutate(p_hat = 1/(1 + exp(-logit_hat))) %>%
                mutate(me = beta_ym*(p_hat)*(1-p_hat))

predictions %>% ggplot(aes(x=ym, y=me, group=sex, colour=sex)) +  
  geom_line(size=2, alpha = 1) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
guides(fill="none", color="none") +
  theme(axis.text.y = element_text(size = 22), axis.text.x = element_text(size = 22),
  axis.title.y = element_text(size = 24), axis.title.x = element_text(size = 24), 
  legend.text = element_text(size = 18), legend.position="bottom") +
  labs(x="years of marriage", y="efecto 'marginal' years of marriage") 
```
]
---
## Efectos marginales sobre la probabilidad de éxito

- Efectos marginales son _esencialmente_ heterogeneos. No hay un efecto sino muchos. 

--

- Heterogeneidad crece con la complejidad del modelo: número de predictores, interacciones, etc. 

--

- En la práctica, muchas veces queremos UN número que resuma el efecto marginal. 

<br>
--
.pull-left[
![For god sakes just give me the damn number](https://i.makeagif.com/media/8-29-2018/ior4IF.gif)
]

--

Cantidades de interes:
.pull-right[

* Average Marginal Effects (AME)

* Marginal Effects at the Mean (MEM)

* Marginal Effects at Representative Values (MER)

]

---
## Efectos marginales sobre la probabilidad de éxito: AME

.bold[AME]: Efecto marginal promedio (en muestra)

--

En `R` pueden usar el nuevo paquete `marginaleffects` (o comando `margins` en `Stata` )

```{r}
avg_slopes(logit_affairs_sex_ym)
```

--

.pull-left[
.bold[¿De donde vienen estos números?]

]

.img-right[
![boric](boric.png)
]

---
## Efectos marginales sobre la probabilidad de éxito: AME

 .bold[(a) Solución análítica], usando 
$$\frac{\partial p_{i}}{\partial \text{ym} }= \beta_{\text{ym}} \cdot p_{i}(1 - p_{i})$$
--

En `R`:

```{r}
beta_ym = logit_affairs_sex_ym$coefficients[3] # coeficiente ym 
p_hat <- predict(logit_affairs_sex_ym, type="response") # prob predichas cada obs.
me_ym = tibble(me_ym=beta_ym*p_hat*(1-p_hat)) # efector marginal para cada obs.
```

--

.pull-left[
```{r, echo=FALSE, fig.width=4.5, fig.height=3.7, message=FALSE}
me_ym
```
]

.pull-right[
```{r}
ame = round(mean(me_ym$me_ym),4)
print(paste0("AME Years of Marriage: ", ame))
```
]

---
## Efectos marginales sobre la probabilidad de éxito: AME


 .bold[(b) Aproximación numérica], usando la definición la derivada:
--

$$\frac{\partial p_{i}}{\partial \text{ym}} \approx  \frac{p_{i}(x_{1}, \dots ,\text{ym} = c + \delta) - p_{i}(x_{1}, \dots ,\text{ym} = c )}{\delta}$$
<br>
--

```{r}
p_hat <- predict(logit_affairs_sex_ym, type="response") # prob predichas cada obs.

delta = 0.01 # cantidad cambio
affairsdata_delta <- affairsdata %>% mutate(ym = ym + delta) # cambio marginal en ym
# probabilidad predicha con ym aumentado marginalmente
p_hat_delta <- predict(logit_affairs_sex_ym, newdata=affairsdata_delta ,type="response") 

me_ym_approx <- (p_hat_delta - p_hat)/delta # approx. numerica marginal effect.

# average marginal effect
round(mean(me_ym_approx),4)
```

---
## Efectos marginales sobre la probabilidad de éxito: MEM

.bold[MEM]: Efecto marginal al valor promedio (en muestra)

--

.pull-left[
- ¿Qué promedio?

  - ¿Promedio de "years of marriage" para cada género por separado?
  
  - ¿Promedio global de "years of marriage" por cada género por separado?
  
  - ¿Promedio global de "years of marriage" y la media género (como var. continua)?
]

--

.pull-right[

![ah](http://25.media.tumblr.com/tumblr_lx16y1gdGm1qls3mfo1_500.gif)


]


--

Depende de que lo que queremos ...

---
## Efectos marginales sobre la probabilidad de éxito: MEM

Implentamos la tercera opción usando `marginaleffects` en `R`: efecto marginal para años de matrimonio y género fijos a su promedio globa.

--

```{r, warning=FALSE, message=FALSE}
affairsdata <- affairsdata %>% mutate(male = 1*(sex=="male"))
logit_affairs_sex_ym2 <- glm(everaffair_d ~ male + ym, family=binomial(link="logit"), 
                             data=affairsdata)

mean_male = mean(affairsdata$male);  mean_ym = mean(affairsdata$ym)
c("Promedio ym"=mean_ym, "Promedio male"=mean_male)
```

--

```{r, warning=FALSE, message=FALSE}
marginaleffects(logit_affairs_sex_ym2, newdata =   expand.grid(male=mean_male, ym=mean_ym))
```

---
## Efectos marginales sobre la probabilidad de éxito: MEM

.bold[Solución analística], usando 
$$\frac{\partial p_{i}}{\partial \text{ym} }= \beta_{\text{ym}} \cdot p_{i}(1 - p_{i})$$

--

```{r}
beta_ym = logit_affairs_sex_ym$coefficients[3] # coeficiente ym 
grid <- affairsdata %>% data_grid(male=mean_male, ym=mean_ym, .model=logit_affairs_sex_ym2); grid
p_hat <- predict(logit_affairs_sex_ym2, newdata=grid, type="response") # prob predichas 

mem_ym = beta_ym*p_hat*(1-p_hat); round(mem_ym,4) # efecto marginal 

```

---
## Efectos marginales sobre la probabilidad de éxito: MEM

.bold[Aproximación numérica], usando: 

$$\frac{\partial p_{i}}{\partial \text{ym}} \approx  \frac{p_{i}(x_{1}, \dots ,\text{ym} = c + \delta) - p_{i}(x_{1}, \dots ,\text{ym} = c )}{\delta}$$

<br>
--

```{r}
delta = 0.01 # cantidad cambio

grid <- affairsdata %>% data_grid(male=mean_male, ym=c(mean_ym, mean_ym+ delta), .model=logit_affairs_sex_ym2)
p_hat <- predict(logit_affairs_sex_ym2, newdata=grid, type="response") # prob predichas 

mem_ym_approx = (p_hat[2] - p_hat[1])/delta
round(mem_ym_approx,4) # efecto marginal 
```

---
## Efectos marginales sobre la probabilidad de éxito: MER

Una versión más general de MEM es calcular los efectos marginales para "valores representativos". 



--

En `R`


```{r, warning=FALSE, message=FALSE}
marginaleffects(logit_affairs_sex_ym, 
                        newdata =  expand.grid(sex = c("male", "female"),
                                  ym=c(4,7,15)
                                  )) %>% filter(term=="ym")
```




---
class: inverse, center, middle


##Hasta la próxima clase. Gracias!


<br>
Mauricio Bucca <br>
https://mebucca.github.io/ <br>
github.com/mebucca




