Análisis de Datos Categóricos (SOC3070)

class: center, middle, inverse, title-slide

.title[
# Análisis de Datos Categóricos (SOC3070)
]
.subtitle[
## Efectos en Regresión Logística
]
.author[
### Mauricio Bucca Profesor Asistente, Sociología UC
]
.date[
### <a href="https://github.com/mebucca">github.com/mebucca</a>
]

---

class: inverse, center, middle

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

---
##Configuración

- 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}$`

--

donde,

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

---
##Configuración

- 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}$`

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

- 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}$`

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}}$$`
 
--

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

---
## Un ejemplo empírico

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

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

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)
```

```
## 
## Call:
## glm(formula = everaffair_d ~ factor(sex) + ym, family = binomial(link = "logit"), 
## data = affairsdata)
## 
## Coefficients:
## Estimate Std. Error z value Pr(>|z|) 
## (Intercept) -1.71395 0.20633 -8.307 < 2e-16 ***
## factor(sex)male 0.22157 0.19064 1.162 0.245139 
## ym 0.05843 0.01727 3.383 0.000718 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
## Null deviance: 675.38 on 600 degrees of freedom
## Residual deviance: 662.15 on 598 degrees of freedom
## AIC: 668.15
## 
## Number of Fisher Scoring iterations: 4
```
]

---
class: inverse, center, middle

## Efectos marginales sobre el logit

---
## Efectos marginales sobre el logit

.center[ 
![](class_10_files/figure-html/unnamed-chunk-3-1.png)
]

---
## 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}$`

.pull-left[

```
## Estimate Std. Error
## (Intercept) -1.71395429 0.20633183
## factor(sex)male 0.22157111 0.19064190
## ym 0.05842959 0.01727402
```
]
.pull-right[
![](class_10_files/figure-html/unnamed-chunk-5-1.png)
]

---
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}$$`

--

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$`,

`$$\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}}$`.

--

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

`$$\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}$`

--

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}}$`"
]

.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)]
```

```
##                    Estimate Std. Error
## (Intercept)     -1.71395429 0.20633183
## factor(sex)male  0.22157111 0.19064190
## ym               0.05842959 0.01727402
```

``` r
# exp(coeffs)
exp(summary(logit_affairs_sex_ym)$coefficients[,c(1)])
```

```
##     (Intercept) factor(sex)male              ym 
##        0.180152        1.248036        1.060170
```
    
]

--
 
.pull-right[
 ![](class_10_files/figure-html/unnamed-chunk-7-1.png)
]

---
class: inverse, center, middle

## Efectos marginales sobre la probabilidad de éxito

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

.center[ 
![](class_10_files/figure-html/unnamed-chunk-8-1.png)
]

---
## 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}$$`

--
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)
]

--

- .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}$$`

---
## 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}$$`

--

- .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}$$`
 
--

- .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}$$`

--

- .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)$$`

--

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.

--

- 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

--

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[

---
## 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})}}$$`

![](class_10_files/figure-html/unnamed-chunk-11-1.png)
]

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

![](class_10_files/figure-html/unnamed-chunk-12-1.png)
]
---
## 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.

--
.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)
```

```
## 
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 %
## sex mean(male) - mean(female) 0.0407 0.03505 1.16 0.245 2.0 -0.02799
## ym mean(dY/dX) 0.0107 0.00308 3.47 <0.001 10.9 0.00466
## 97.5 %
## 0.1094
## 0.0167
## 
## Type: response 
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
```

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

]

---
## 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") # probabilidades predichas

me_ym <- tibble(p = p_hat, beta_ym = beta_ym, me = beta_ym * p_hat * (1 - p_hat) 
)
```

.pull-left[

```
## # A tibble: 601 × 3
## p beta_ym me
## <dbl> <dbl> <dbl>
## 1 0.2874 0.05843 0.01197 
## 2 0.1854 0.05843 0.008824
## 3 0.3021 0.05843 0.01232 
## 4 0.3507 0.05843 0.01331 
## 5 0.1902 0.05843 0.009000
## 6 0.1643 0.05843 0.008024
## 7 0.1584 0.05843 0.007789
## 8 0.3507 0.05843 0.01331 
## 9 0.3021 0.05843 0.01232 
## 10 0.1971 0.05843 0.009245
## # ℹ 591 more rows
```
]

.pull-right[

``` r
ame = round(mean(me_ym$me),4)
print(paste0("AME Years of Marriage: ", ame))
```

```
## [1] "AME Years of Marriage: 0.0107"
```
]

---
## 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}$$`
 
--

``` 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)
```

```
## [1] 0.0107
```

---
## 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
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)
```

```
##   Promedio ym Promedio male 
##     8.1776955     0.4758735
```

``` r
slopes(logit_affairs_sex_ym2, newdata =   expand.grid(male=mean_male, ym=mean_ym))
```

```
## 
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## male 1 - 0 0.0410 0.03528 1.16 0.245 2.0 -0.02818 0.1101
## ym dY/dX 0.0108 0.00315 3.43 <0.001 10.7 0.00461 0.0169
## 
## Type: response 
## Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, male, ym, everaffair_d
```

---
## 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
slopes(logit_affairs_sex_ym, 
                        newdata =  expand.grid(sex = c("male", "female"),
                                  ym=c(4,7,15)
                                  )) %>% filter(term=="ym")
```

```
## 
## Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0.01007 0.00259 3.89 < 0.001 13.3 0.00499 0.0151
## 0.00882 0.00226 3.91 < 0.001 13.4 0.00440 0.0132
## 0.01104 0.00315 3.51 < 0.001 11.1 0.00487 0.0172
## 0.00981 0.00281 3.49 < 0.001 11.0 0.00430 0.0153
## 0.01331 0.00437 3.05 0.00232 8.8 0.00474 0.0219
## 0.01232 0.00420 2.93 0.00335 8.2 0.00409 0.0205
## 
## Term: ym
## Type: response 
## Comparison: dY/dX
## Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, sex, ym, everaffair_d
```

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

| Método | Definición | Ventaja | Limitación |
|--------|------------|---------|------------|
| **AME** | Promedio del efecto marginal sobre todos los casos | Resume heterogeneidad en 1 número | Puede ocultar variación |
| **MEM** | Efecto en covariables fijadas en su promedio | Fácil de interpretar | “Promedio” puede no ser realista |
| **MER** | Efectos en valores representativos | Ilustra heterogeneidad | Depende de la elección de escenarios |

---
class: inverse, center, middle

##Hasta la próxima clase. Gracias!

Mauricio Bucca 
https://mebucca.github.io/ 
github.com/mebucca