Multinomial logistic regression

Fit models for women’s labor force participation

library(car)
library(nnet)   # for `multinom()`
library(dplyr)
data(Womenlf, package = "carData")

make not.work the reference category

Womenlf <- within (Womenlf, {
  partic <- ordered(partic, levels=c('not.work', 'parttime', 'fulltime'))})

Fit model with main effects

wlf.multinom <- multinom(partic ~ hincome + children, 
                         data=Womenlf, Hess=TRUE)

## # weights:  12 (6 variable)
## initial  value 288.935032 
## iter  10 value 211.454772
## final  value 211.440963 
## converged

summary(wlf.multinom, Wald=TRUE)

## Call:
## multinom(formula = partic ~ hincome + children, data = Womenlf, 
##     Hess = TRUE)
## 
## Coefficients:
##          (Intercept)      hincome childrenpresent
## parttime   -1.432321  0.006893838      0.02145558
## fulltime    1.982842 -0.097232073     -2.55860537
## 
## Std. Errors:
##          (Intercept)    hincome childrenpresent
## parttime   0.5924627 0.02345484       0.4690352
## fulltime   0.4841789 0.02809599       0.3621999
## 
## Value/SE (Wald statistics):
##          (Intercept)    hincome childrenpresent
## parttime   -2.417573  0.2939197      0.04574407
## fulltime    4.095266 -3.4607098     -7.06407045
## 
## Residual Deviance: 422.8819 
## AIC: 434.8819

Anova(wlf.multinom)

## Analysis of Deviance Table (Type II tests)
## 
## Response: partic
##          LR Chisq Df Pr(>Chisq)    
## hincome    15.153  2  0.0005123 ***
## children   63.559  2  1.579e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# overall test?
#car::linearHypothesis(wlf.multinom, "parttime, fulltime")

Examine coefficients

coef(wlf.multinom)

##          (Intercept)      hincome childrenpresent
## parttime   -1.432321  0.006893838      0.02145558
## fulltime    1.982842 -0.097232073     -2.55860537

exp(coef(wlf.multinom))

##          (Intercept)   hincome childrenpresent
## parttime    0.238754 1.0069177      1.02168740
## fulltime    7.263353 0.9073454      0.07741263

Wald tests & p-values

stats <- summary(wlf.multinom, Wald=TRUE)
z <- stats$Wald.ratios
p <- 2 * (1 - pnorm(abs(z)))
zapsmall(p)

##          (Intercept)   hincome childrenpresent
## parttime   0.0156244 0.7688193       0.9635142
## fulltime   0.0000422 0.0005388       0.0000000

Plot fitted values in two panels

op <- par(mfrow=c(1,2))

predictors <- expand.grid(hincome=1:45, children=c('absent', 'present'))
p.fit <- predict(wlf.multinom, predictors, type='probs')
Hinc <- 1:max(predictors$hincome)
for ( kids in c("absent", "present") ) {
    data <- subset(data.frame(predictors, p.fit), children==kids)
    plot( range(Hinc), c(0,1), type="n",
        xlab="Husband's Income", ylab='Fitted Probability',
        main = paste("Children", kids))
    lines(Hinc, data[, 'not.work'], lwd=3, col="black", lty=1)
    lines(Hinc, data[, 'parttime'], lwd=3, col="blue",  lty=2)
    lines(Hinc, data[, 'fulltime'], lwd=3, col="red",   lty=3)
  if (kids=="absent") {
  legend(5, 0.97, lty=1:3, lwd=3, col=c("black", "blue", "red"),
    legend=c('not working', 'part-time', 'full-time'))
    }
}

par(op)

Effect plots

wlf.effects <- allEffects(wlf.multinom)
plot(wlf.effects, style='stacked')

plot(Effect(c("hincome", "children"), wlf.multinom), 
     style = "stacked", key.args=list(x=.75, y=.25),
     colors = c(grey(.85), "pink", "lightblue")
     )

plotting probabilities

get fitted probabilities

options(digits=3)
predictors <- expand.grid(hincome=1:50, children=c('absent', 'present'))
fit <- data.frame(predictors, 
                  predict(wlf.multinom, predictors, type='probs'))

Re-shape for plotting

library(tidyr)

tidyr::gather()

plotdat <- fit |>
  gather(key="Level", value="Probability", not.work:fulltime) 


plotdat <- fit |>
  pivot_longer(not.work:fulltime, 
               names_to = "Level",
               values_to = "Probability") 
head(plotdat)

## # A tibble: 6 x 4
##   hincome children Level    Probability
##     <int> <fct>    <chr>          <dbl>
## 1       1 absent   not.work      0.128 
## 2       1 absent   parttime      0.0307
## 3       1 absent   fulltime      0.842 
## 4       2 absent   not.work      0.138 
## 5       2 absent   parttime      0.0335
## 6       2 absent   fulltime      0.828

library(directlabels)
gg <-
ggplot(plotdat, aes(x = hincome, y = Probability, colour = Level)) + 
  geom_line(size=1.5) + theme_bw() + 
  facet_grid(~ children, labeller = label_both) +
  theme_bw(base_size = 14)
direct.label(gg, list("top.bumptwice", dl.trans(y = y + 0.2)))