##       partic hincome children   region
## 3   not.work      45  present  Ontario
## 56  not.work      17   absent Atlantic
## 94  parttime       9  present  Ontario
## 98  fulltime      15   absent  Ontario
## 108 not.work      19  present  Ontario
## 171 fulltime       7  present  Prairie
## 177 not.work      15  present  Ontario
## 252 not.work      23   absent   Quebec

Recode to create dichotomies

Create binary variables working (working or not working), and fulltime (only for those working).

Womenlf <- Womenlf |>
  mutate(working = ifelse(partic=="not.work", 0, 1)) |>
  mutate(fulltime = case_when(
    working & partic == "fulltime" ~ 1,
    working & partic == "parttime" ~ 0)
  )
some(Womenlf, 8)
##       partic hincome children   region working fulltime
## 4   not.work      23  present  Ontario       0       NA
## 9   not.work      15  present  Ontario       0       NA
## 30  fulltime      14  present  Prairie       1        1
## 114 not.work       1  present Atlantic       0       NA
## 136 fulltime      11  present  Ontario       1        1
## 154 parttime      28  present       BC       1        0
## 177 not.work      15  present  Ontario       0       NA
## 180 fulltime      13   absent       BC       1        1

Fit models for each dichotomy

Womenlf <- within(Womenlf, contrasts(children)<- 'contr.treatment')
mod.working <- glm(working ~ hincome + children, family=binomial, data=Womenlf)
mod.fulltime <- glm(fulltime ~ hincome + children, family=binomial, data=Womenlf)

summary(mod.working)
## 
## Call:
## glm(formula = working ~ hincome + children, family = binomial, 
##     data = Womenlf)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.677  -0.865  -0.777   0.929   1.997  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       1.3358     0.3838    3.48   0.0005 ***
## hincome          -0.0423     0.0198   -2.14   0.0324 *  
## childrenpresent  -1.5756     0.2923   -5.39    7e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 356.15  on 262  degrees of freedom
## Residual deviance: 319.73  on 260  degrees of freedom
## AIC: 325.7
## 
## Number of Fisher Scoring iterations: 4
summary(mod.fulltime)
## 
## Call:
## glm(formula = fulltime ~ hincome + children, family = binomial, 
##     data = Womenlf)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.405  -0.868   0.395   0.621   1.764  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.4778     0.7671    4.53  5.8e-06 ***
## hincome          -0.1073     0.0392   -2.74   0.0061 ** 
## childrenpresent  -2.6515     0.5411   -4.90  9.6e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 144.34  on 107  degrees of freedom
## Residual deviance: 104.49  on 105  degrees of freedom
##   (155 observations deleted due to missingness)
## AIC: 110.5
## 
## Number of Fisher Scoring iterations: 5
lmtest::coeftest(mod.working)
## 
## z test of coefficients:
## 
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       1.3358     0.3838    3.48   0.0005 ***
## hincome          -0.0423     0.0198   -2.14   0.0324 *  
## childrenpresent  -1.5756     0.2923   -5.39    7e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lmtest::coeftest(mod.fulltime)
## 
## z test of coefficients:
## 
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.4778     0.7671    4.53  5.8e-06 ***
## hincome          -0.1073     0.0392   -2.74   0.0061 ** 
## childrenpresent  -2.6515     0.5411   -4.90  9.6e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod.working)  
## Analysis of Deviance Table (Type II tests)
## 
## Response: working
##          LR Chisq Df Pr(>Chisq)    
## hincome      4.83  1      0.028 *  
## children    31.32  1    2.2e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod.fulltime)
## Analysis of Deviance Table (Type II tests)
## 
## Response: fulltime
##          LR Chisq Df Pr(>Chisq)    
## hincome       9.0  1     0.0027 ** 
## children     32.1  1    1.4e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prepare for plotting

Create a grid of values of hincome and children. Generate predicted values from the models

attach(Womenlf)
predictors <- expand.grid(hincome=1:45, children=c('absent', 'present'))
# get fitted values for both sub-models
p.work     <- predict(mod.working, predictors, type='response')
p.fulltime <- predict(mod.fulltime, predictors, type='response')

calculate unconditional probs for the three response categories

p.full <- p.work * p.fulltime
p.part <- p.work * (1 - p.fulltime)
p.not <- 1 - p.work

NB: the plot looks best if the plot window is made wider than tall

op <- par(mfrow=c(1,2))
# children absent panel
plot(c(1,45), c(0,1), 
    type='n', xlab="Husband's Income", ylab='Fitted Probability',
    main="Children absent")
lines(1:45, p.not[1:45],  lty=1, lwd=3, col="black")
lines(1:45, p.part[1:45], lty=2, lwd=3, col="blue")
lines(1:45, p.full[1:45], lty=3, lwd=3, col="red")

legend(5, 0.97, lty=1:3, lwd=3, col=c("black", "blue", "red"),
    legend=c('not working', 'part-time', 'full-time'))  

# children present panel
plot(c(1,45), c(0,1), 
    type='n', xlab="Husband's Income", ylab='Fitted Probability',
    main="Children present")
lines(1:45, p.not[46:90],  lty=1, lwd=3, col="black")
lines(1:45, p.part[46:90], lty=2, lwd=3, col="blue")
lines(1:45, p.full[46:90], lty=3, lwd=3, col="red")

par(op)

A more general way to make the plot

op <- par(mfrow=c(1,2))
plotdata <- data.frame(predictors, p.full, p.part, p.not)
Hinc <- 1:max(plotdata$hincome)
for ( kids in c("absent", "present") ) {
    data <- subset(plotdata, children==kids)
    plot( range(data$hincome), c(0,1), type="n",
        xlab="Husband's Income", ylab='Fitted Probability',
        main = paste("Children", kids),
        cex.lab = 1.3)
    lines(Hinc, data$p.not,  lwd=3, col="black", lty=1)
    lines(Hinc, data$p.part, lwd=3, col="blue",  lty=2)
    lines(Hinc, data$p.full, lwd=3, col="red",   lty=3)
  if (kids=="absent") {
  legend(15, 0.97, lty=1:3, lwd=3, col=c("black", "blue", "red"),
    legend=c('not working', 'part-time', 'full-time'))
    }
}

par(op)


detach(Womenlf)
IycgLS0tDQojJyB0aXRsZTogIlBvbHl0b21vdXMgRGF0YTogbmVzdGVkIGRpY2hvdG9taWVzIg0KIycgYXV0aG9yOiAiTWljaGFlbCBGcmllbmRseSINCiMnIGRhdGU6ICJgciBmb3JtYXQoU3lzLkRhdGUoKSlgIg0KIycgb3V0cHV0Og0KIycgICBodG1sX2RvY3VtZW50Og0KIycgICAgIHRoZW1lOiByZWFkYWJsZQ0KIycgICAgIGNvZGVfZG93bmxvYWQ6IHRydWUNCiMnIC0tLQ0KDQoNCiMrIGVjaG89RkFMU0UNCmtuaXRyOjpvcHRzX2NodW5rJHNldCgNCiAgd2FybmluZyA9IEZBTFNFLCAgICMgYXZvaWQgd2FybmluZ3MgYW5kIG1lc3NhZ2VzIGluIHRoZSBvdXRwdXQNCiAgbWVzc2FnZSA9IEZBTFNFDQopDQoNCmxpYnJhcnkoY2FyKSAgICMgZm9yIGRhdGEgYW5kIEFub3ZhKCkNCmxpYnJhcnkoZHBseXIpDQpkYXRhKFdvbWVubGYsIHBhY2thZ2UgPSAiY2FyRGF0YSIpDQpzb21lKFdvbWVubGYsIDgpDQoNCg0KIycgIyMgUmVjb2RlIHRvIGNyZWF0ZSBkaWNob3RvbWllcw0KIycgQ3JlYXRlIGJpbmFyeSB2YXJpYWJsZXMgYHdvcmtpbmdgICh3b3JraW5nIG9yIG5vdCB3b3JraW5nKSwgYW5kIGBmdWxsdGltZWAgKG9ubHkgZm9yIHRob3NlIHdvcmtpbmcpLg0KV29tZW5sZiA8LSBXb21lbmxmIHw+DQogIG11dGF0ZSh3b3JraW5nID0gaWZlbHNlKHBhcnRpYz09Im5vdC53b3JrIiwgMCwgMSkpIHw+DQogIG11dGF0ZShmdWxsdGltZSA9IGNhc2Vfd2hlbigNCiAgICB3b3JraW5nICYgcGFydGljID09ICJmdWxsdGltZSIgfiAxLA0KICAgIHdvcmtpbmcgJiBwYXJ0aWMgPT0gInBhcnR0aW1lIiB+IDApDQogICkNCnNvbWUoV29tZW5sZiwgOCkNCg0KDQojJyAjIyBGaXQgbW9kZWxzIGZvciBlYWNoIGRpY2hvdG9teQ0KV29tZW5sZiA8LSB3aXRoaW4oV29tZW5sZiwgY29udHJhc3RzKGNoaWxkcmVuKTwtICdjb250ci50cmVhdG1lbnQnKQ0KbW9kLndvcmtpbmcgPC0gZ2xtKHdvcmtpbmcgfiBoaW5jb21lICsgY2hpbGRyZW4sIGZhbWlseT1iaW5vbWlhbCwgZGF0YT1Xb21lbmxmKQ0KbW9kLmZ1bGx0aW1lIDwtIGdsbShmdWxsdGltZSB+IGhpbmNvbWUgKyBjaGlsZHJlbiwgZmFtaWx5PWJpbm9taWFsLCBkYXRhPVdvbWVubGYpDQoNCnN1bW1hcnkobW9kLndvcmtpbmcpDQpzdW1tYXJ5KG1vZC5mdWxsdGltZSkNCg0KbG10ZXN0Ojpjb2VmdGVzdChtb2Qud29ya2luZykNCmxtdGVzdDo6Y29lZnRlc3QobW9kLmZ1bGx0aW1lKQ0KDQpBbm92YShtb2Qud29ya2luZykgIA0KQW5vdmEobW9kLmZ1bGx0aW1lKQ0KDQojJyAjIyBQcmVwYXJlIGZvciBwbG90dGluZw0KIycgQ3JlYXRlIGEgZ3JpZCBvZiB2YWx1ZXMgb2YgYGhpbmNvbWVgIGFuZCBgY2hpbGRyZW5gLiBHZW5lcmF0ZSBwcmVkaWN0ZWQgdmFsdWVzIGZyb20gdGhlIG1vZGVscw0KYXR0YWNoKFdvbWVubGYpDQpwcmVkaWN0b3JzIDwtIGV4cGFuZC5ncmlkKGhpbmNvbWU9MTo0NSwgY2hpbGRyZW49YygnYWJzZW50JywgJ3ByZXNlbnQnKSkNCiMgZ2V0IGZpdHRlZCB2YWx1ZXMgZm9yIGJvdGggc3ViLW1vZGVscw0KcC53b3JrICAgICA8LSBwcmVkaWN0KG1vZC53b3JraW5nLCBwcmVkaWN0b3JzLCB0eXBlPSdyZXNwb25zZScpDQpwLmZ1bGx0aW1lIDwtIHByZWRpY3QobW9kLmZ1bGx0aW1lLCBwcmVkaWN0b3JzLCB0eXBlPSdyZXNwb25zZScpDQoNCiMnIGNhbGN1bGF0ZSB1bmNvbmRpdGlvbmFsIHByb2JzIGZvciB0aGUgdGhyZWUgcmVzcG9uc2UgY2F0ZWdvcmllcw0KcC5mdWxsIDwtIHAud29yayAqIHAuZnVsbHRpbWUNCnAucGFydCA8LSBwLndvcmsgKiAoMSAtIHAuZnVsbHRpbWUpDQpwLm5vdCA8LSAxIC0gcC53b3JrDQoNCiMnIE5COiB0aGUgcGxvdCBsb29rcyBiZXN0IGlmIHRoZSBwbG90IHdpbmRvdyBpcyBtYWRlIHdpZGVyIHRoYW4gdGFsbA0Kb3AgPC0gcGFyKG1mcm93PWMoMSwyKSkNCiMgY2hpbGRyZW4gYWJzZW50IHBhbmVsDQpwbG90KGMoMSw0NSksIGMoMCwxKSwgDQogICAgdHlwZT0nbicsIHhsYWI9Ikh1c2JhbmQncyBJbmNvbWUiLCB5bGFiPSdGaXR0ZWQgUHJvYmFiaWxpdHknLA0KICAgIG1haW49IkNoaWxkcmVuIGFic2VudCIpDQpsaW5lcygxOjQ1LCBwLm5vdFsxOjQ1XSwgIGx0eT0xLCBsd2Q9MywgY29sPSJibGFjayIpDQpsaW5lcygxOjQ1LCBwLnBhcnRbMTo0NV0sIGx0eT0yLCBsd2Q9MywgY29sPSJibHVlIikNCmxpbmVzKDE6NDUsIHAuZnVsbFsxOjQ1XSwgbHR5PTMsIGx3ZD0zLCBjb2w9InJlZCIpDQoNCmxlZ2VuZCg1LCAwLjk3LCBsdHk9MTozLCBsd2Q9MywgY29sPWMoImJsYWNrIiwgImJsdWUiLCAicmVkIiksDQogICAgbGVnZW5kPWMoJ25vdCB3b3JraW5nJywgJ3BhcnQtdGltZScsICdmdWxsLXRpbWUnKSkgIA0KDQojIGNoaWxkcmVuIHByZXNlbnQgcGFuZWwNCnBsb3QoYygxLDQ1KSwgYygwLDEpLCANCiAgICB0eXBlPSduJywgeGxhYj0iSHVzYmFuZCdzIEluY29tZSIsIHlsYWI9J0ZpdHRlZCBQcm9iYWJpbGl0eScsDQogICAgbWFpbj0iQ2hpbGRyZW4gcHJlc2VudCIpDQpsaW5lcygxOjQ1LCBwLm5vdFs0Njo5MF0sICBsdHk9MSwgbHdkPTMsIGNvbD0iYmxhY2siKQ0KbGluZXMoMTo0NSwgcC5wYXJ0WzQ2OjkwXSwgbHR5PTIsIGx3ZD0zLCBjb2w9ImJsdWUiKQ0KbGluZXMoMTo0NSwgcC5mdWxsWzQ2OjkwXSwgbHR5PTMsIGx3ZD0zLCBjb2w9InJlZCIpDQpwYXIob3ApDQoNCiMnICMjIEEgbW9yZSBnZW5lcmFsIHdheSB0byBtYWtlIHRoZSBwbG90DQpvcCA8LSBwYXIobWZyb3c9YygxLDIpKQ0KcGxvdGRhdGEgPC0gZGF0YS5mcmFtZShwcmVkaWN0b3JzLCBwLmZ1bGwsIHAucGFydCwgcC5ub3QpDQpIaW5jIDwtIDE6bWF4KHBsb3RkYXRhJGhpbmNvbWUpDQpmb3IgKCBraWRzIGluIGMoImFic2VudCIsICJwcmVzZW50IikgKSB7DQoJZGF0YSA8LSBzdWJzZXQocGxvdGRhdGEsIGNoaWxkcmVuPT1raWRzKQ0KCXBsb3QoIHJhbmdlKGRhdGEkaGluY29tZSksIGMoMCwxKSwgdHlwZT0ibiIsDQoJCXhsYWI9Ikh1c2JhbmQncyBJbmNvbWUiLCB5bGFiPSdGaXR0ZWQgUHJvYmFiaWxpdHknLA0KCQltYWluID0gcGFzdGUoIkNoaWxkcmVuIiwga2lkcyksDQoJCWNleC5sYWIgPSAxLjMpDQoJbGluZXMoSGluYywgZGF0YSRwLm5vdCwgIGx3ZD0zLCBjb2w9ImJsYWNrIiwgbHR5PTEpDQoJbGluZXMoSGluYywgZGF0YSRwLnBhcnQsIGx3ZD0zLCBjb2w9ImJsdWUiLCAgbHR5PTIpDQoJbGluZXMoSGluYywgZGF0YSRwLmZ1bGwsIGx3ZD0zLCBjb2w9InJlZCIsICAgbHR5PTMpDQogIGlmIChraWRzPT0iYWJzZW50Iikgew0KICBsZWdlbmQoMTUsIDAuOTcsIGx0eT0xOjMsIGx3ZD0zLCBjb2w9YygiYmxhY2siLCAiYmx1ZSIsICJyZWQiKSwNCiAgICBsZWdlbmQ9Yygnbm90IHdvcmtpbmcnLCAncGFydC10aW1lJywgJ2Z1bGwtdGltZScpKQ0KICAgIH0NCn0NCnBhcihvcCkNCg0KDQpkZXRhY2goV29tZW5sZikNCiAgDQo=