multinomial.r

``` r
## Multinomial Logit Example

## Libraries
library(tidyverse) # Graphing and other tools
library(palmerpenguins) # The data
library(nnet)  # Fitting the model
library(marginaleffects) # Extracting the effects (v0.8 or higher)

### Note:
## The marginaleffects package has many examples and a good 
## discussion of the many sorts of conditional effects you can 
## extract and draw: https://vincentarelbundock.github.io/marginaleffects/

## My book has a more general discussion of strategies for plotting 
## model effects stuff in general. It uses a slightly older package
## for calculating marginal effects, but the general points still apply.
## See e.g. https://socviz.co/modeling.html for stuff on models specifically.

## Data
penguins |>
  select(species, island, body_mass_g, sex)
#> # A tibble: 344 × 4
#>    species island    body_mass_g sex   
#>    <fct>   <fct>           <int> <fct> 
#>  1 Adelie  Torgersen        3750 male  
#>  2 Adelie  Torgersen        3800 female
#>  3 Adelie  Torgersen        3250 female
#>  4 Adelie  Torgersen          NA <NA>  
#>  5 Adelie  Torgersen        3450 female
#>  6 Adelie  Torgersen        3650 male  
#>  7 Adelie  Torgersen        3625 female
#>  8 Adelie  Torgersen        4675 male  
#>  9 Adelie  Torgersen        3475 <NA>  
#> 10 Adelie  Torgersen        4250 <NA>  
#> # … with 334 more rows

## We'll try to predict `species`
penguins |>
  count(species)
#> # A tibble: 3 × 2
#>   species       n
#>   <fct>     <int>
#> 1 Adelie      152
#> 2 Chinstrap    68
#> 3 Gentoo      124


out <-  multinom(species ~ body_mass_g + sex,
                 data = penguins, trace = FALSE)

summary(out)
#> Call:
#> multinom(formula = species ~ body_mass_g + sex, data = penguins, 
#>     trace = FALSE)
#> 
#> Coefficients:
#>           (Intercept)  body_mass_g     sexmale
#> Chinstrap   -1.685209 0.0002690784  -0.1598625
#> Gentoo     -75.467266 0.0187355090 -14.4242193
#> 
#> Std. Errors:
#>            (Intercept)  body_mass_g    sexmale
#> Chinstrap 0.0498123789 5.005832e-05 0.31011613
#> Gentoo    0.0008968176 1.542868e-04 0.00488104
#> 
#> Residual Deviance: 282.982 
#> AIC: 294.982


## Extract the average marginal effects
ame <- marginaleffects(out, type = "probs")

ame |>
  summary()
#>       Group        Term      Contrast     Effect Std. Error   z value
#> 1    Adelie body_mass_g         dY/dX -1.125e-04  1.233e-05   -9.1277
#> 2    Adelie         sex male - female  1.835e-01  3.965e-02    4.6285
#> 3 Chinstrap body_mass_g         dY/dX -4.926e-06  7.648e-06   -0.6441
#> 4 Chinstrap         sex male - female  6.938e-02  4.158e-02    1.6686
#> 5    Gentoo body_mass_g         dY/dX  1.175e-04  1.695e-05    6.9310
#> 6    Gentoo         sex male - female -2.529e-01  2.321e-03 -108.9954
#>     Pr(>|z|)      2.5 %     97.5 %
#> 1 < 2.22e-16 -1.367e-04 -8.838e-05
#> 2 3.6827e-06  1.058e-01  2.613e-01
#> 3   0.519522 -1.992e-05  1.006e-05
#> 4   0.095201 -1.212e-02  1.509e-01
#> 5 4.1794e-12  8.425e-05  1.507e-04
#> 6 < 2.22e-16 -2.575e-01 -2.484e-01
#> 
#> Model type:  multinom 
#> Prediction type:  probs


## We get one marginal effect for each term. So now we can plot by group.
## We'll use the `plot_cap()` function from the `marginaleffects` pacakge.
## It will plot adjusted predictions against the predictors. We could
## also do this manually via ggplot.
out |>
  plot_cap(condition = "body_mass_g",
           type = "probs") +
  facet_wrap(~ group, ncol = 1) +
  labs(title = "Conditional  Adjusted Predictions",
       x = "Body Mass in grams",
       y = "Probability")
```

![](https://i.imgur.com/OKFPlEd.png)

``` r

## Conditional Contrasts
plot_cco(
  out,
  effect = list(body_mass_g = c(3550, 4750)),
  condition = "group",
  type = "probs") +

  labs(title = "Conditional Contrasts Example",
       x = NULL,
       y = "Contrast in Body Mass on Species") +
  coord_flip()
```

![](https://i.imgur.com/yIierML.png)

``` r



## There are other ways to do it ... I guess the
## hard part is calculating the exact contrast/
## effect/comparison you want. Presenting it
## is mostly a matter of trying to help the
## reader/viewer grasp the logic of the comparison
## and understand the uncertainty involved in the
## estimate.
```

<sup>Created on 2022-10-08 with [reprex v2.0.2](https://reprex.tidyverse.org)</sup>