library(car) # for data and Anova()
data(Womenlf)
some(Womenlf)
## partic hincome children region
## 29 not.work 17 present Prairie
## 33 not.work 15 present Ontario
## 36 not.work 19 absent Ontario
## 51 parttime 10 present Prairie
## 58 not.work 23 present Atlantic
## 104 not.work 15 absent BC
## 128 parttime 13 present Prairie
## 143 not.work 11 present Prairie
## 208 fulltime 11 absent Quebec
## 235 not.work 7 present Quebec
Womenlf <- within(Womenlf,{
working <- recode(partic, " 'not.work' = 'no'; else = 'yes' ")
fulltime <- recode (partic,
" 'fulltime' = 'yes'; 'parttime' = 'no'; 'not.work' = NA")})
some(Womenlf)
## partic hincome children region fulltime working
## 31 not.work 13 present Ontario <NA> no
## 69 not.work 19 present Ontario <NA> no
## 82 parttime 15 present Ontario no yes
## 122 not.work 23 present Atlantic <NA> no
## 128 parttime 13 present Prairie no yes
## 137 parttime 14 present Ontario no yes
## 143 not.work 11 present Prairie <NA> no
## 204 not.work 28 present Quebec <NA> no
## 207 not.work 11 present Quebec <NA> no
## 216 parttime 17 present Quebec no yes
Womenlf <- within(Womenlf, contrasts(children)<- 'contr.treatment')
mod.working <- glm(working ~ 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.6767 -0.8652 -0.7768 0.9292 1.9970
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.33583 0.38376 3.481 0.0005 ***
## hincome -0.04231 0.01978 -2.139 0.0324 *
## childrenpresent -1.57565 0.29226 -5.391 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.73
##
## Number of Fisher Scoring iterations: 4
mod.fulltime <- glm(fulltime ~ hincome + children, family=binomial, data=Womenlf)
summary(mod.fulltime)
##
## Call:
## glm(formula = fulltime ~ hincome + children, family = binomial,
## data = Womenlf)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4047 -0.8678 0.3949 0.6213 1.7641
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.47777 0.76711 4.534 5.80e-06 ***
## hincome -0.10727 0.03915 -2.740 0.00615 **
## childrenpresent -2.65146 0.54108 -4.900 9.57e-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.49
##
## Number of Fisher Scoring iterations: 5
Anova(mod.working)
## Analysis of Deviance Table (Type II tests)
##
## Response: working
## LR Chisq Df Pr(>Chisq)
## hincome 4.8264 1 0.02803 *
## children 31.3229 1 2.185e-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 8.981 1 0.002728 **
## children 32.136 1 1.437e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Prepare for plotting
attach(Womenlf)
# get fitted values for both sub-models
predictors <- expand.grid(hincome=1:45, children=c('absent', 'present'))
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)
op <- par(mfrow=c(1,2))
plotdata <- data.frame(predictors, p.full, p.part, p.not)
Hinc <- 1:max(hincome)
for ( kids in c("absent", "present") ) {
data <- subset(plotdata, children==kids)
plot( range(hincome), c(0,1), type="n",
xlab="Husband's Income", ylab='Fitted Probability',
main = paste("Children", kids))
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(5, 0.97, lty=1:3, lwd=3, col=c("black", "blue", "red"),
legend=c('not working', 'part-time', 'full-time'))
}
}
par(op)
detach(Womenlf)