load data

Use data from the carData package. Canadian Women’s Labour-Force Participation. The Womenlf data frame has 263 rows and 4 columns. The data are from a 1977 survey of the Canadian population.

library(carData)
data("Womenlf")

# Labour-Force Participation
summary(Womenlf$partic)
## fulltime not.work parttime 
##       66      155       42

model

Model partic as a function of hincome (husband’s income) and children (presence of children in household).

library(nnet)
m <- multinom(partic ~ hincome + children, data = Womenlf, 
              Hess = TRUE)
## # weights:  12 (6 variable)
## initial  value 288.935032 
## iter  10 value 211.441198
## final  value 211.440963 
## converged

Notice the fitted values are probabilities that sum to 1. The first subject would be classified as “not.work” since that category has the highest probability.

head(fitted(m))
##      fulltime  not.work  parttime
## 1 0.093327569 0.7136254 0.1930471
## 2 0.111423199 0.7014260 0.1871508
## 3 0.005281102 0.7464419 0.2482770
## 4 0.044637912 0.7429858 0.2123763
## 5 0.064856573 0.7316822 0.2034612
## 6 0.184773776 0.6490593 0.1661669

using a matrix on lhs

Assume that instead of selecting a category, the subjects had allocated percentages to the categories, similar to what Jim describes. Create a left hand side matrix.

Y <- round(fitted(m) * 100)
head(Y)
##   fulltime not.work parttime
## 1        9       71       19
## 2       11       70       19
## 3        1       75       25
## 4        4       74       21
## 5        6       73       20
## 6       18       65       17

Now model the matrix as a function of hincome and children.

m2 <- multinom(Y ~ hincome + children, 
               data = Womenlf,
               Hess = TRUE)
## # weights:  12 (6 variable)
## initial  value 28837.473965 
## iter  10 value 21061.859596
## final  value 21061.611503 
## converged

The fitted values are still probabilities that sum to 1.

head(fitted(m2))
##      fulltime  not.work  parttime
## 1 0.091931590 0.7149080 0.1931604
## 2 0.109910197 0.7027802 0.1873096
## 3 0.005100011 0.7471615 0.2477385
## 4 0.043727966 0.7439781 0.2122940
## 5 0.063708944 0.7328136 0.2034774
## 6 0.183023680 0.6505448 0.1664315

visualize model

Probability of “preferring” fulltime decreases as husband’s income increases. Likewise, probability of “preferring” not.work increases as husband’s income increases

library(ggeffects)
plot(ggpredict(m2, terms = "hincome"))