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.
data(Vietnam)
Vtn_tab = xtabs(Freq ~ ., data = Vietnam)
cotabplot( ~ year + response | sex, data = Vtn_tab, gp = shading_Marimekko(Vtn_tab[1,,]))

Above graphs illustrate a linear relationship of year and response. The situation for male students(choose C and D) is more clear than female students from the graphs.

library(MASS)
Vietnam$response = ordered(Vietnam$response)
Vtn_polr = polr(response ~ year * sex, data = Vietnam, weights = Freq, Hess = TRUE)
summary(Vtn_polr)
## Call:
## polr(formula = response ~ year * sex, data = Vietnam, weights = Freq, 
##     Hess = TRUE)
## 
## Coefficients:
##                Value Std. Error t value
## year          0.1008    0.05481   1.840
## sexMale      -1.4367    0.22301  -6.442
## year:sexMale  0.2333    0.06025   3.873
## 
## Intercepts:
##     Value   Std. Error t value
## A|B -1.3533  0.2059    -6.5737
## B|C -0.2538  0.2041    -1.2437
## C|D  2.1879  0.2104    10.3968
## 
## Residual Deviance: 7757.059 
## AIC: 7769.059
Anova(Vtn_polr)
## Analysis of Deviance Table (Type II tests)
## 
## Response: response
##          LR Chisq Df Pr(>Chisq)    
## year      166.510  1  < 2.2e-16 ***
## sex        58.473  1  2.061e-14 ***
## year:sex   15.041  1  0.0001052 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Thus, all terms are significant, including the interaction term with sex

  1. 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.
Vtn_PO = vglm(response ~ year * sex, data = Vietnam, weights = Freq, family = cumulative(parallel = TRUE))
Vtn_NPO = vglm(response ~ year * sex, data = Vietnam, weights = Freq, family = cumulative(parallel = FALSE))
lrtest(Vtn_NPO, Vtn_PO)
## Likelihood ratio test
## 
## Model 1: response ~ year * sex
## Model 2: response ~ year * sex
##   #Df  LogLik Df Chisq Pr(>Chisq)    
## 1 108 -3839.7                        
## 2 114 -3878.5  6 77.57  1.134e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The above results show that the unconstrained model fits better than the proportional odds model for this data set.

  1. Fit the multinomial logistic model, also allowing an interaction. Use car::Anova () to assess the model terms.
library(nnet)
Vtn_mtn = multinom(response ~ year * sex, data = Vietnam, weights = Freq)
## # weights:  20 (12 variable)
## initial  value 4362.668354 
## iter  10 value 3927.047674
## final  value 3838.854903 
## converged
summary(Vtn_mtn)
## Call:
## multinom(formula = response ~ year * sex, data = Vietnam, weights = Freq)
## 
## Coefficients:
##   (Intercept)          year    sexMale year:sexMale
## B   0.3928759 -0.0005171856 -0.9719739    0.1469216
## C   1.1452464  0.1475132943 -1.8898266    0.1708179
## D  -1.3100155  0.1908925232 -2.0083782    0.4550651
## 
## Std. Errors:
##   (Intercept)       year   sexMale year:sexMale
## B   0.3914932 0.11020875 0.4105455   0.11627237
## C   0.3368612 0.09373079 0.3572585   0.09985147
## D   0.5677546 0.15097871 0.6163987   0.16217684
## 
## Residual Deviance: 7677.71 
## AIC: 7701.71
Anova(Vtn_mtn)
## Analysis of Deviance Table (Type II tests)
## 
## Response: response
##          LR Chisq Df Pr(>Chisq)    
## year      175.884  3    < 2e-16 ***
## sex       137.627  3    < 2e-16 ***
## year:sex    7.751  3    0.05145 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. Produce an effect plot for this model and describe the nature of the interaction.
library(effects)
## Loading required package: carData
## 
## Attaching package: 'carData'
## The following objects are masked from 'package:car':
## 
##     Guyer, UN, Vocab
## The following object is masked from 'package:vcdExtra':
## 
##     Burt
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
plot(allEffects(Vtn_polr), style = "stacked")

The majority of female students prefer more peaceful options, whereas a majority of male students prefer more violent options.

  1. 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.
Vietnam = within(Vietnam, yearMale <- year * (sex == "Male"))
Vtn_mtn2 = multinom(response ~ sex + yearMale, data = Vietnam, weights = Freq)
## # weights:  16 (9 variable)
## initial  value 4362.668354 
## iter  10 value 3894.170452
## final  value 3841.568268 
## converged
Anova(Vtn_mtn2)
## Analysis of Deviance Table (Type II tests)
## 
## Response: response
##          LR Chisq Df Pr(>Chisq)    
## sex        242.80  3  < 2.2e-16 ***
## yearMale   178.21  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Vtn_mtn, Vtn_mtn2)
## Likelihood ratio tests of Multinomial Models
## 
## Response: response
##            Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 sex + yearMale       111   7683.137                                
## 2     year * sex       108   7677.710 1 vs 2     3  5.42673 0.1430872

The simpler model is not significant worse than the general model with interaction.