Exercise 8.4 Refer to Homework #5, second question for a description of the Vietnam data set in vcdExtra. The goal here is to fit models for the polytomous response variable in relation to year and sex.

  1. Fit the proportional odds model to these data, allowing an interaction of year and sex.
  2. Is there evidence that the proportional odds assumption does not hold for this data set? Use the methods described in “Ordinal response: Proportional odds model” to assess this.
  3. Fit the multinomial logistic model, also allowing an interaction. Use car::Anova () to assess the model terms.
  4. Produce an effect plot for this model and describe the nature of the interaction.
  5. Fit the simpler multinomial model in which there is no effect of year for females and the effect of year is linear for males (on the logit scale). Test whether this model is significantly worse than the general multinomial model with interaction.
data("Vietnam", package = "vcdExtra")
head(Vietnam)
##      sex year response Freq
## 1 Female    1        A   13
## 2 Female    1        B   19
## 3 Female    1        C   40
## 4 Female    1        D    5
## 5 Female    2        A    5
## 6 Female    2        B    9
Vietnam$response = ordered(Vietnam$response)
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
vietpo = vglm(response ~ year +sex , family=cumulative(parallel = TRUE),data = Vietnam)

vietpo
## 
## Call:
## vglm(formula = response ~ year + sex, family = cumulative(parallel = TRUE), 
##     data = Vietnam)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3          year       sexMale 
## -1.098612e+00 -6.879800e-16  1.098612e+00 -8.881784e-17  1.256074e-16 
## 
## Degrees of Freedom: 120 Total; 115 Residual
## Residual deviance: 110.9035 
## Log-likelihood: -55.45177
vietnpo = vglm(response ~ year +sex , family=cumulative(parallel = FALSE),data = Vietnam)

vietnpo
## 
## Call:
## vglm(formula = response ~ year + sex, family = cumulative(parallel = FALSE), 
##     data = Vietnam)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3        year:1        year:2 
## -1.098612e+00 -6.879800e-16  1.098612e+00  4.634791e-16  6.326856e-17 
##        year:3     sexMale:1     sexMale:2     sexMale:3 
##  8.599751e-17 -3.918388e-16 -1.133916e-15 -8.918715e-16 
## 
## Degrees of Freedom: 120 Total; 111 Residual
## Residual deviance: 110.9035 
## Log-likelihood: -55.45177
coef(vietpo ,matrix=TRUE)
##             logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3])
## (Intercept)  -1.098612e+00  -6.879800e-16   1.098612e+00
## year         -8.881784e-17  -8.881784e-17  -8.881784e-17
## sexMale       1.256074e-16   1.256074e-16   1.256074e-16
coef(vietnpo ,matrix=TRUE)
##             logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3])
## (Intercept)  -1.098612e+00  -6.879800e-16   1.098612e+00
## year          4.634791e-16   6.326856e-17   8.599751e-17
## sexMale      -3.918388e-16  -1.133916e-15  -8.918715e-16
library(nnet)
viet.multinom = multinom(response ~ year +sex ,data = Vietnam, Hess= TRUE)
## # weights:  16 (9 variable)
## initial  value 55.451774 
## final  value 55.451774 
## converged
summary(viet.multinom )
## Call:
## multinom(formula = response ~ year + sex, data = Vietnam, Hess = TRUE)
## 
## Coefficients:
##   (Intercept) year sexMale
## B           0    0       0
## C           0    0       0
## D           0    0       0
## 
## Std. Errors:
##   (Intercept)      year   sexMale
## B    1.140175 0.3162278 0.8944272
## C    1.140175 0.3162278 0.8944272
## D    1.140175 0.3162278 0.8944272
## 
## Residual Deviance: 110.9035 
## AIC: 128.9035
library(car)
## 
## Attaching package: 'car'
## The following object is masked from 'package:VGAM':
## 
##     logit
Anova(viet.multinom)
## Analysis of Deviance Table (Type II tests)
## 
## Response: response
##      LR Chisq Df Pr(>Chisq)
## year        0  3          1
## sex         0  3          1