Load packages

library(car)
library(nnet)   # for `multinom()`
library(dplyr)
library(effects)
library(ggplot2)

Fit models for women’s labor force participation

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

Summaries

summary(wlf.multinom, Wald=TRUE)
## Call:
## multinom(formula = partic ~ hincome + children, data = Womenlf, 
##     Hess = TRUE)
## 
## Coefficients:
##          (Intercept)  hincome childrenpresent
## parttime       -1.43  0.00689          0.0215
## fulltime        1.98 -0.09723         -2.5586
## 
## Std. Errors:
##          (Intercept) hincome childrenpresent
## parttime       0.592  0.0235           0.469
## fulltime       0.484  0.0281           0.362
## 
## Value/SE (Wald statistics):
##          (Intercept) hincome childrenpresent
## parttime       -2.42   0.294          0.0457
## fulltime        4.10  -3.461         -7.0641
## 
## Residual Deviance: 423 
## AIC: 435

broom::tidy() puts this in a more convenient format

broom::tidy(wlf.multinom)
## # A tibble: 6 x 6
##   y.level  term            estimate std.error statistic  p.value
##   <chr>    <chr>              <dbl>     <dbl>     <dbl>    <dbl>
## 1 parttime (Intercept)     -1.43       0.592    -2.42   1.56e- 2
## 2 parttime hincome          0.00689    0.0235    0.294  7.69e- 1
## 3 parttime childrenpresent  0.0215     0.469     0.0457 9.64e- 1
## 4 fulltime (Intercept)      1.98       0.484     4.10   4.22e- 5
## 5 fulltime hincome         -0.0972     0.0281   -3.46   5.39e- 4
## 6 fulltime childrenpresent -2.56       0.362    -7.06   1.62e-12

Anova tests

Anova(wlf.multinom)
## Analysis of Deviance Table (Type II tests)
## 
## Response: partic
##          LR Chisq Df Pr(>Chisq)    
## hincome      15.2  2    0.00051 ***
## children     63.6  2    1.6e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Overall test B = 0

Test the hypothesis that all coefficients = 0. For linearHypothesis(), need to construct the combinations of the dichotomies and the predictors.

(H <- outer(c("parttime:", "fulltime:"),
           c("childrenpresent", "hincome"),
           paste0))
##      [,1]                       [,2]              
## [1,] "parttime:childrenpresent" "parttime:hincome"
## [2,] "fulltime:childrenpresent" "fulltime:hincome"
car::linearHypothesis(wlf.multinom, as.vector(H)) 
## Linear hypothesis test
## 
## Hypothesis:
## parttime:childrenpresent = 0
## fulltime:childrenpresent = 0
## parttime:hincome = 0
## fulltime:hincome = 0
## 
## Model 1: restricted model
## Model 2: partic ~ hincome + children
## 
##   Df Chisq Pr(>Chisq)    
## 1                        
## 2  4  58.4    6.2e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Examine coefficients

coef(wlf.multinom)
##          (Intercept)  hincome childrenpresent
## parttime       -1.43  0.00689          0.0215
## fulltime        1.98 -0.09723         -2.5586
exp(coef(wlf.multinom))
##          (Intercept) hincome childrenpresent
## parttime       0.239   1.007          1.0217
## fulltime       7.263   0.907          0.0774

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.016   0.769           0.964
## fulltime       0.000   0.001           0.000

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

IycgLS0tDQojJyB0aXRsZTogIk11bHRpbm9taWFsIGxvZ2lzdGljIHJlZ3Jlc3Npb24iDQojJyBhdXRob3I6ICJNaWNoYWVsIEZyaWVuZGx5Ig0KIycgZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiDQojJyBvdXRwdXQ6DQojJyAgIGh0bWxfZG9jdW1lbnQ6DQojJyAgICAgdGhlbWU6IHJlYWRhYmxlDQojJyAgICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KIycgLS0tDQoNCiMrIGVjaG89RkFMU0UNCmtuaXRyOjpvcHRzX2NodW5rJHNldCgNCiAgd2FybmluZyA9IEZBTFNFLCAgICMgYXZvaWQgd2FybmluZ3MgYW5kIG1lc3NhZ2VzIGluIHRoZSBvdXRwdXQNCiAgbWVzc2FnZSA9IEZBTFNFDQopDQoNCiMnICMjIExvYWQgcGFja2FnZXMNCmxpYnJhcnkoY2FyKQ0KbGlicmFyeShubmV0KSAgICMgZm9yIGBtdWx0aW5vbSgpYA0KbGlicmFyeShkcGx5cikNCmxpYnJhcnkoZWZmZWN0cykNCmxpYnJhcnkoZ2dwbG90MikNCg0KIycgIyMgRml0IG1vZGVscyBmb3Igd29tZW4ncyBsYWJvciBmb3JjZSBwYXJ0aWNpcGF0aW9uDQpkYXRhKFdvbWVubGYsIHBhY2thZ2UgPSAiY2FyRGF0YSIpDQoNCiMnICMjIyMgbWFrZSBgbm90LndvcmtgIHRoZSByZWZlcmVuY2UgY2F0ZWdvcnkNCg0KV29tZW5sZiA8LSB3aXRoaW4gKFdvbWVubGYsIHsNCiAgcGFydGljIDwtIG9yZGVyZWQocGFydGljLCBsZXZlbHM9Yygnbm90LndvcmsnLCAncGFydHRpbWUnLCAnZnVsbHRpbWUnKSl9KQ0KDQojJyAjIyBGaXQgbW9kZWwgd2l0aCBtYWluIGVmZmVjdHMNCndsZi5tdWx0aW5vbSA8LSBtdWx0aW5vbShwYXJ0aWMgfiBoaW5jb21lICsgY2hpbGRyZW4sIA0KICAgICAgICAgICAgICAgICAgICAgICAgIGRhdGE9V29tZW5sZiwgSGVzcz1UUlVFKQ0KIycgIyMgU3VtbWFyaWVzDQpzdW1tYXJ5KHdsZi5tdWx0aW5vbSwgV2FsZD1UUlVFKQ0KDQojJyBgYnJvb206OnRpZHkoKWAgcHV0cyB0aGlzIGluIGEgbW9yZSBjb252ZW5pZW50IGZvcm1hdA0KYnJvb206OnRpZHkod2xmLm11bHRpbm9tKQ0KDQojJyAjIyBBbm92YSB0ZXN0cw0KQW5vdmEod2xmLm11bHRpbm9tKQ0KDQojJyAjIyBPdmVyYWxsIHRlc3QgQiA9IDANCiMnIFRlc3QgdGhlIGh5cG90aGVzaXMgdGhhdCBhbGwgY29lZmZpY2llbnRzID0gMC4NCiMnIEZvciBgbGluZWFySHlwb3RoZXNpcygpYCwgbmVlZCB0byBjb25zdHJ1Y3QgdGhlIGNvbWJpbmF0aW9ucyBvZg0KIycgdGhlIGRpY2hvdG9taWVzIGFuZCB0aGUgcHJlZGljdG9ycy4NCihIIDwtIG91dGVyKGMoInBhcnR0aW1lOiIsICJmdWxsdGltZToiKSwNCiAgICAgICAgICAgYygiY2hpbGRyZW5wcmVzZW50IiwgImhpbmNvbWUiKSwNCiAgICAgICAgICAgcGFzdGUwKSkNCg0KY2FyOjpsaW5lYXJIeXBvdGhlc2lzKHdsZi5tdWx0aW5vbSwgYXMudmVjdG9yKEgpKSANCg0KDQojJyAjIyBFeGFtaW5lIGNvZWZmaWNpZW50cw0KY29lZih3bGYubXVsdGlub20pDQpleHAoY29lZih3bGYubXVsdGlub20pKQ0KDQoNCiMnICMjIFdhbGQgdGVzdHMgJiBwLXZhbHVlcw0Kc3RhdHMgPC0gc3VtbWFyeSh3bGYubXVsdGlub20sIFdhbGQ9VFJVRSkNCnogPC0gc3RhdHMkV2FsZC5yYXRpb3MNCnAgPC0gMiAqICgxIC0gcG5vcm0oYWJzKHopKSkNCnphcHNtYWxsKHApDQoNCg0KIycgIyMgUGxvdCBmaXR0ZWQgdmFsdWVzIGluIHR3byBwYW5lbHMNCm9wIDwtIHBhcihtZnJvdz1jKDEsMikpDQoNCnByZWRpY3RvcnMgPC0gZXhwYW5kLmdyaWQoaGluY29tZT0xOjQ1LCBjaGlsZHJlbj1jKCdhYnNlbnQnLCAncHJlc2VudCcpKQ0KcC5maXQgPC0gcHJlZGljdCh3bGYubXVsdGlub20sIHByZWRpY3RvcnMsIHR5cGU9J3Byb2JzJykNCkhpbmMgPC0gMTptYXgocHJlZGljdG9ycyRoaW5jb21lKQ0KZm9yICgga2lkcyBpbiBjKCJhYnNlbnQiLCAicHJlc2VudCIpICkgew0KCWRhdGEgPC0gc3Vic2V0KGRhdGEuZnJhbWUocHJlZGljdG9ycywgcC5maXQpLCBjaGlsZHJlbj09a2lkcykNCglwbG90KCByYW5nZShIaW5jKSwgYygwLDEpLCB0eXBlPSJuIiwNCgkJeGxhYj0iSHVzYmFuZCdzIEluY29tZSIsIHlsYWI9J0ZpdHRlZCBQcm9iYWJpbGl0eScsDQoJCW1haW4gPSBwYXN0ZSgiQ2hpbGRyZW4iLCBraWRzKSkNCglsaW5lcyhIaW5jLCBkYXRhWywgJ25vdC53b3JrJ10sIGx3ZD0zLCBjb2w9ImJsYWNrIiwgbHR5PTEpDQoJbGluZXMoSGluYywgZGF0YVssICdwYXJ0dGltZSddLCBsd2Q9MywgY29sPSJibHVlIiwgIGx0eT0yKQ0KCWxpbmVzKEhpbmMsIGRhdGFbLCAnZnVsbHRpbWUnXSwgbHdkPTMsIGNvbD0icmVkIiwgICBsdHk9MykNCiAgaWYgKGtpZHM9PSJhYnNlbnQiKSB7DQogIGxlZ2VuZCg1LCAwLjk3LCBsdHk9MTozLCBsd2Q9MywgY29sPWMoImJsYWNrIiwgImJsdWUiLCAicmVkIiksDQogICAgbGVnZW5kPWMoJ25vdCB3b3JraW5nJywgJ3BhcnQtdGltZScsICdmdWxsLXRpbWUnKSkNCiAgICB9DQp9DQpwYXIob3ApDQoNCg0KDQojJyAjIyBFZmZlY3QgcGxvdHMNCiMnIA0Kd2xmLmVmZmVjdHMgPC0gYWxsRWZmZWN0cyh3bGYubXVsdGlub20pDQpwbG90KHdsZi5lZmZlY3RzLCBzdHlsZT0nc3RhY2tlZCcpDQoNCnBsb3QoRWZmZWN0KGMoImhpbmNvbWUiLCAiY2hpbGRyZW4iKSwgd2xmLm11bHRpbm9tKSwgDQogICAgIHN0eWxlID0gInN0YWNrZWQiLCBrZXkuYXJncz1saXN0KHg9Ljc1LCB5PS4yNSksDQogICAgIGNvbG9ycyA9IGMoZ3JleSguODUpLCAicGluayIsICJsaWdodGJsdWUiKQ0KICAgICApDQoNCg0KIycgIyMgcGxvdHRpbmcgcHJvYmFiaWxpdGllcw0KDQojJyBnZXQgZml0dGVkIHByb2JhYmlsaXRpZXMNCm9wdGlvbnMoZGlnaXRzPTMpDQpwcmVkaWN0b3JzIDwtIGV4cGFuZC5ncmlkKGhpbmNvbWU9MTo1MCwgY2hpbGRyZW49YygnYWJzZW50JywgJ3ByZXNlbnQnKSkNCmZpdCA8LSBkYXRhLmZyYW1lKHByZWRpY3RvcnMsIA0KICAgICAgICAgICAgICAgICAgcHJlZGljdCh3bGYubXVsdGlub20sIHByZWRpY3RvcnMsIHR5cGU9J3Byb2JzJykpDQoNCg0KIycgIyMgUmUtc2hhcGUgZm9yIHBsb3R0aW5nDQpsaWJyYXJ5KHRpZHlyKQ0KDQojJyB0aWR5cjo6Z2F0aGVyKCkNCnBsb3RkYXQgPC0gZml0IHw+DQogIGdhdGhlcihrZXk9IkxldmVsIiwgdmFsdWU9IlByb2JhYmlsaXR5Iiwgbm90Lndvcms6ZnVsbHRpbWUpIA0KDQoNCnBsb3RkYXQgPC0gZml0IHw+DQogIHBpdm90X2xvbmdlcihub3Qud29yazpmdWxsdGltZSwgDQogICAgICAgICAgICAgICBuYW1lc190byA9ICJMZXZlbCIsDQogICAgICAgICAgICAgICB2YWx1ZXNfdG8gPSAiUHJvYmFiaWxpdHkiKSANCmhlYWQocGxvdGRhdCkNCg0KbGlicmFyeShkaXJlY3RsYWJlbHMpDQpnZyA8LQ0KZ2dwbG90KHBsb3RkYXQsIGFlcyh4ID0gaGluY29tZSwgeSA9IFByb2JhYmlsaXR5LCBjb2xvdXIgPSBMZXZlbCkpICsgDQogIGdlb21fbGluZShzaXplPTEuNSkgKyB0aGVtZV9idygpICsgDQogIGZhY2V0X2dyaWQofiBjaGlsZHJlbiwgbGFiZWxsZXIgPSBsYWJlbF9ib3RoKSArDQogIHRoZW1lX2J3KGJhc2Vfc2l6ZSA9IDE0KQ0KZGlyZWN0LmxhYmVsKGdnLCBsaXN0KCJ0b3AuYnVtcHR3aWNlIiwgZGwudHJhbnMoeSA9IHkgKyAwLjIpKSkNCg0KDQo=