class: inverse, center, middle

```{r xaringan-themer, include=FALSE, warning=FALSE}
library(tidyverse)
library(xaringanthemer)
style_duo_accent(primary_color ="#556B2F", secondary_color = "#317589",
                 background_color = "#f8f7f3",
                 header_font_google = google_font("Archivo"),
                 text_font_google   = google_font("Inconsolata"), 
                 link_color= "#C19A6B"

)
```

# Regresión Logística Multinomial


---
## Estructura de un modelo de regresión logística multinomial 

$$\newcommand{\vect}[1]{\boldsymbol{#1}}$$
Un modelo de regresión logística multinomial generaliza la regresión logística (binomial) a situaciones en que la variable dependiente es una .bold[variable discreta con dos o más valores no ordenados] (ejemplo: voto entre tres candidatos, elección de barrio, etc).

<br>
--

.bold[Configuración]

- Tenemos $n$ observaciones (individuos) independientes: $i = 1, \dots, n$

--

- Para cada observación observamos datos $y_{i}, \dots , y_{n}$ que actúan como variable dependiente, donde $y_{i} \in \{j:1,2, \cdots, J\}$
  
  - Las $J$ categorías de $y_{i}$ no siguen necesariamente un orden.

--

- Asumimos que estos datos son realizaciones de $n$ variables aleatorias que siguen una distribución Multinomial con probabilidades desconocidas: $Y_{i} \sim \text{Multinomial}(\vect{p_{i}}: p_{1}, \dots, p_{J})$

--

- Dichas probabilidades, $\vect{p_{i}}$, varían de individuo en individuo en función de ciertas covariables.


---
## Distribución Multinomial

<br>
.bold[Distribución Multinomial:]

- Dado un experimento con resultados posibles $\{1,2, \dots, J\}$, con respectiva probabilidad de "éxito" $\{p_1,p_2, \dots, p_J\}$ 

--

- Si repetimos el experimento $n$ veces: ¿Cuál es la probabilidad de obtener la siguiente cantidad de éxitos en cada categoría: $\{n_1,n_2, \dots, n_J\}$, donde $n_1 + n_2 + \dots + n_J =n$?
  
<br>
--

.bold[Ejemplo]: Al tirar un dado justo 12 veces, ¿cuál es la probabilidad de obtener cada lado 2 veces?

```{r}
dmultinom(x=c(2,2,2,2,2,2), size=12 ,p=c(1/6,1/6,1/6,1/6,1/6,1/6))
```

---
## Distribución Multinomial (n=1)


.bold[Distribución Multinoulli o Categórica:]

Si la cantidad de intentos es igual a 1 (n=1), la probabilidad de éxito en una de las categorías (y fracaso en todas las otras) viene definida por la siguiente función de probabilidad:
 
--

$$\quad \mathbb{P}(Y = j ) = p_{1}^{1[y=1]}p_{2}^{1[y=2]} \cdots p_{J}^{1[y=J]}    \quad \text{ donde } \quad y_{j} \in \{0,1\} \quad \text{y} \quad \sum^{J}_{j=1}1[y=j]=1$$

--

.bold[Ejemplo]: Al tirar un dado justo 1 vez, ¿cuál es la probabilidad de obtener el número 4?

```{r}
dmultinom(x=c(0,0,0,1,0,0), size=1 ,p=c(1/6,1/6,1/6,1/6,1/6,1/6))
```

--

- En una regresión logística multinomial cada observación en la variable dependiente ( $y_{i}$ ) viene de una distribución Multinoulli.

--

- Observamos los resultados ( $y_{i}=j$ ) y queremos estimar las probabilidades que los generan, ( $p_{ij}$ ).

--
 
- Específicamente, queremos estimar ( $p_{ij}$ ) como un función de covariables , usando número de parámetros $k$, tal que $k<n$.



---
## Regresión logística binomial


Una regresión logística binomial $y_{i} \in {0,1}$

.content-box-yellow[
$$\ln \frac{\mathbb{P}(y_{i}=1)}{\mathbb{P}(y_{i}=0)} =   \ln \frac{p_{i}}{1 - p_{i}} = \beta_{0} + \beta_{1}x_{1i} + \dots + \beta_{k}x_{ki}$$] 


<br>
--

Equivalentemente, podemos re-escribir el modelo como

<br>

$$\ln \frac{\mathbb{P}(y_{i}=1)}{\mathbb{P}(y_{i}=0)} =   \ln \frac{p_{1i}}{p_{0i}} =  \beta_{0} + \beta_{1}x_{1i} + \dots + \beta_{k}x_{ki}$$
<br>
--

Dado que $p_{0} + p_{1} = 1$, sólo necesitamos estimar una equación.


---
## Regresión logística multinomial, 3 categorías 


Si la variable dependiente toma tres valores, $y \in \{1,2,3\}$, y las probabilidades de obtener estos resultados son $p_{1},p_{2},p_{3}$, entonces podemos estimar estas tres probabilidades con dos modelos de regresión logística:

<br>
--

Usando $p_{3}$ como categoría de referencia:


$$(1) \quad \ln \frac{\mathbb{P}(y_{i}=1)}{\mathbb{P}(y_{i}=3)} =   \ln \frac{p_{1i}}{p_{3i}} =  \beta_{10} + \beta_{11}x_{1i} + \dots + \beta_{1k}x_{ki}$$
--
<br> y

$$(2) \quad \ln \frac{\mathbb{P}(y_{i}=2)}{\mathbb{P}(y_{i}=3)} =   \ln \frac{p_{2i}}{p_{3i}} =  \beta_{20} + \beta_{21}x_{1i} + \dots + \beta_{2k}x_{ki}$$
<br>
--

Describimos $\{(p_{1},p_{2},p_{3})_{i}, \quad  \dots \quad ,(p_{1},p_{2},p_{3})_{n}\}$ con $2(k + 1)$ parámetros.

---
## Regresión logística multinomial, J-categorías 

Generalizando, si la variable dependiente toma, $J$ valores -- $y \in \{1,2,\dots,J\}$ -- , y las probabilidades de obtener estos resultados son $p_{1},p_{2},\dots,p_{J}$, entonces podemos estimar estas $J$ probabilidades con  $J-1$ modelos de regresión logística:

--

Usando $p_{J}$ como categoría de referencia:

$$(1) \quad \ln \frac{\mathbb{P}(y_{i}=1)}{\mathbb{P}(y_{i}=J)} =   \ln \frac{p_{1i}}{p_{Ji}} =  \beta_{10} + \beta_{11}x_{1i} + \dots + \beta_{1k}x_{ki}$$
--

$$(2) \quad \ln \frac{\mathbb{P}(y_{i}=2)}{\mathbb{P}(y_{i}=J)} =   \ln \frac{p_{2i}}{p_{Ji}} =  \beta_{20} + \beta_{21}x_{1i} + \dots + \beta_{2k}x_{ki}$$
--

$$\vdots$$
--

$$(J-1) \quad \ln \frac{\mathbb{P}(y_{i}=J-1)}{\mathbb{P}(y_{i}=J)} =   \ln \frac{p_{(J-1)i}}{p_{Ji}} =  \beta_{(J-1)0} + \beta_{(J-1)1}x_{1i} + \dots + \beta_{(J-1)k}x_{ki}$$

---
## Regresión logística multinomial, J-categorías 

Una versión más compacta de lo expresado anteriormente lleva a la formulación estándar de un modelo de regresión multinomial:

<br>
--

.content-box-yellow[
$$\underbrace{\ln \frac{\mathbb{P}(y_{i}=j)}{\mathbb{P}(y_{i}=J)}}_{\text{log odds = log relative probability}}=   \ln \frac{p_{ji}}{p_{Ji}} = \overbrace{ \beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}^{\text{predictor lineal }j}$$
]

<br>
--

Describimos $\{(p_{1},\dots,p_{J})_{i}, \quad  \dots \quad ,(p_{1},\dots,p_{J})_{n}\}$ con $(J-1)(k+1)$ parámetros.

---
## Regresión logística multinomial, J-categorías 

De la formulación anterior es posible derivar el contraste en todos los restantes pares de categorías.
--

 Consideremos las categorías $a$ y $b$:

--

$$\ln \frac{p_{ai}}{p_{bi}} =  \ln \frac{p_{ai}/p_{Ji}}{p_{bi}/p_{Ji}} = \ln \frac{p_{ai}}{p_{Ji}} - \ln \frac{p_{bi}}{p_{Ji}}$$
--
Por tanto,


$$\ln \frac{p_{ai}}{p_{bi}} =  (\beta_{a0} + \beta_{a1}x_{1i} + \dots + \beta_{ak}x_{ki}) - (\beta_{b0} + \beta_{b1}x_{1i} + \dots + \beta_{bk}x_{ki})$$
--
En resumen:

.content-box-yellow[
$$\ln \frac{p_{ai}}{p_{bi}} =  (\beta_{a0} -\beta_{b0}) + (\beta_{a1} - \beta_{b1}) x_{1i} + \dots + (\beta_{ak} - \beta_{bk}) x_{ki}$$
]

donde cada $(\beta_{ak} -\beta_{bk})$ es un coeficiente en si mismo.
---
## Regresión logística multinomial expresada como probabilidades

Dado

$$\ln \frac{p_{ji}}{p_{Ji}} =  \beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}$$
<br>
--

-  Exponenciando en ambos lados obtenemos la probabilidad "relativa" entre las categorías $j$ y $J$:

$$\frac{p_{ji}}{p_{Ji}} =  e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}$$
<br>
--

- Luego, la probabilidad de obtener la categoría $j$ puede ser expresada como sigue:

$$p_{ji}  =  p_{Ji} \cdot e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}$$

---
## Regresión logística multinomial expresada como probabilidades

La expresión $p_{ji}  =  p_{Ji} \cdot e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}$ significa que:


<br>
--


$$p_{1i}  =  p_{Ji} \cdot e^{\beta_{10} + \beta_{11}x_{1i} + \dots + \beta_{1k}x_{ki}}$$
<br>
--


$$p_{2i}  =  p_{Ji} \cdot e^{\beta_{20} + \beta_{21}x_{1i} + \dots + \beta_{2k}x_{ki}}$$
<br>
--

$$\vdots$$
<br>
--
$$p_{(J-1)i}  =  p_{Ji} \cdot e^{\beta_{(J-1)0} + \beta_{(J-1)1}x_{1i} + \dots + \beta_{(J-1)k}x_{ki}}$$

<br>
--


Luego, falta sólo determinar $p_{Ji}$.

---
## Regresión logística multinomial expresada como probabilidades


Para determinar $p_{Ji}$ usamos los siguientes hechos:

--

.pull-left[
$$\text{(1)} \quad p_{ji}  =  p_{Ji} \cdot e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}$$
]

.pull-right[
$$\text{(2)} \quad   p_{1i} + p_{2i} + \dots + p_{Ji} = 1$$
]

<br>
--

Combinando (1) y (2) obtenemos que:

$$\sum^{J-1}_{j=1} p_{Ji} \cdot e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}} + p_{Ji} = 1, \quad \text{luego ... }$$
--


$$p_{Ji} \bigg(\sum^{J-1}_{j=1}  e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}} + 1 \bigg) = 1, \quad \text{por tanto ... }$$

--

.content-box-yellow[
$$p_{Ji} =\frac{1}{1 + \sum^{J-1}_{j=1}  e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}}$$
]

---
## Regresión logística multinomial expresada como probabilidades

Ahora sabemos que:

--

- $p_{Ji} =\frac{1}{1 + \sum^{J-1}_{j=1}  e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}}$ 

- $p_{ji}  =  p_{Ji} \cdot e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}$

<br>
--
Combinando ambos resultados obtenemos:


.content-box-yellow[
$$p_{ji} =\frac{e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}}{1 + \sum^{J-1}_{j=1}  e^{\beta_{j0} + \beta_{j1}x_{1i} + \dots + \beta_{jk}x_{ki}}}$$
]

---
class: inverse, center, middle

## Estimación

---
## Estimación


<br>
--

- Parámetros son estimados via MLE

--

- Alternativamente, minimizando función $softmax$  (ej, paquete `nnet` en `R`)

--

- Las $J-1$ equaciones del modelo de regresión logística son estimadas simultáneamente, imponiendo la restricción: $\sum_{j}p_{i} = 1$  ("constrained optimization")

--

- Es posible estimar $J-1$ regresiones logísticas separadamente pero no garantiza que $\sum_{j}p_{i} = 1$.



---
class: inverse, center, middle

## Regresión logística multinomial en la práctica

---
## Regresión logística multinomial en la práctica

Para ejemplificar el uso de regresión logística multinomial trabajaremos con los datos de intención de voto en el plebiscito de 1988. 

```{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("rsample")
library("arm")
library("DescTools")
library("caret")
library("carData")
library("nnet")
library("marginaleffects")


theme_set(theme_cowplot())

data(Chile)
plebs_1988 <- Chile %>% as_tibble()

# display the data as a tibble
plebs_1988 %>% arrange(statusquo)
```

<br>
`vote`: (A) Abstención; (N) NO; (U) Indecisa; (Y) SI

`statusquo`: apoyo al status-quo (+)

---
## Regresión logística multinomial en la práctica

```{r, message=F, warning=F, message=FALSE}
mlogit_vote_sq <- multinom(vote ~ statusquo, trace=F, data=plebs_1988); 
summary(mlogit_vote_sq)
```


---
## Regresión logística multinomial en la práctica

```{r, message=F, warning=F, message=FALSE}
mlogit_vote_sq <- multinom(vote ~ statusquo, trace=F, data=plebs_1988)
```

.pull-left[

```{r, message=F, warning=F, message=FALSE}
plebs_1988$intercept <- 1
X <- plebs_1988[, c("intercept", "statusquo")]; X
```
]

.pull-right[
```{r, message=F, warning=F, message=FALSE}
betas <- coef(mlogit_vote_sq); betas
```
]

---
## Regresión logística multinomial en la práctica


$$
\hat{Y} = X\beta  = 
\begin{bmatrix}
1 & x_{11} & x_{12} & \dots & x_{1p} \\
1 & x_{21} & x_{22} & \dots & x_{2p} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n1} & x_{n2} & \dots & x_{np}
\end{bmatrix}
\begin{bmatrix}
\beta_{01} & \beta_{02} & \dots & \beta_{0(K-1)} \\
\beta_{11} & \beta_{12} & \dots & \beta_{1(K-1)} \\
\vdots & \vdots & \ddots & \vdots \\
\beta_{p1} & \beta_{p2} & \dots & \beta_{p(K-1)}
\end{bmatrix}
$$


```{r, message=F, warning=F}
logit_predictions <- as.matrix(X) %*% t(betas)

#Log-odds observation 1
preds_1 <-logit_predictions[1,]; preds_1

#Predicted probabilities observation 1
probs_1 <- exp(preds_1)/(1 + sum(exp(preds_1)))
probs_1
```

---
## Regresión logística multinomial en la práctica

```{r, message=F, warning=F, message=FALSE}
mlogit_vote_sq <- multinom(vote ~ statusquo, trace=F, data=plebs_1988); 
predict(mlogit_vote_sq, type="prob")
```

---
## Regresión logística multinomial en la práctica


.pull-left[
```{r, echo=FALSE,warning=FALSE, message=FALSE}
# plot the result
grid <- plebs_1988  %>% data_grid(statusquo=seq_range(statusquo,30),.model=mlogit_vote_sq)

predictions <- cbind(grid,predict(mlogit_vote_sq, newdata=grid, type="prob")) %>%
              pivot_longer(-statusquo, names_to="outcome",values_to="p" ) %>%
              mutate(p_outcome = paste0("p_",outcome)) %>% dplyr::select(-outcome) %>%
              pivot_wider(names_from = p_outcome, values_from=p) %>%
              mutate(across(p_N:p_Y, ~  log(.x/p_A), .names = "logit_{.col}" )) %>%
              pivot_longer(-statusquo, names_to="quant", values_to="est")

predictions %>% filter(str_detect(quant, "^l")) %>% separate(quant,sep=8, into=c("quant","vote")) %>%
  ggplot(aes(x=statusquo, y=est, colour=vote, group=vote)) +
  geom_path(alpha=0.5, size=1.5) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
  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="top") +
  labs(x="Support status-quo", y="Logit(Vote=j)", colour="vote")
```
]

--

.pull-right[
```{r, echo=FALSE,  warning=FALSE, message=FALSE}
predictions %>% filter(str_detect(quant, "^p")) %>%  separate(quant,sep=2, into=c("quant","vote")) %>%
  ggplot(aes(x=statusquo, y=est, colour=vote, group=vote)) +
  geom_path(alpha=0.5, size=1.5) +
  scale_color_viridis_d() +  scale_fill_viridis_d() +
  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="top") +
  labs(x="Support status-quo", y="P(Vote=j)", colour="vote")
```
]


---
class: inverse, center, middle


##Hasta la próxima clase. Gracias!


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




