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.

library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 3.4.4
## Loading required package: vcd
## Warning: package 'vcd' was built under R version 3.4.4
## Loading required package: grid
## Loading required package: gnm
## Warning: package 'gnm' was built under R version 3.4.4
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
  1. Fit the proportional odds model to these data, allowing an interaction of year and sex.
library(MASS)
## Warning: package 'MASS' was built under R version 3.4.4
library(car)
## Warning: package 'car' was built under R version 3.4.4
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.4.4
## 
## Attaching package: 'carData'
## The following object is masked from 'package:vcdExtra':
## 
##     Burt
Prop.mod=polr(response~year+sex,data = Vietnam)
Anova(Prop.mod)
## Analysis of Deviance Table (Type II tests)
## 
## Response: response
##      LR Chisq Df Pr(>Chisq)
## year        0  1          1
## sex         0  1          1

The proportional odds model shows that there is not any relevant predictors based on response type. ??? (b) 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.

library(VGAM)
## Warning: package 'VGAM' was built under R version 3.4.4
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:car':
## 
##     logit
PO.mod = vglm(response ~ year * sex, data = Vietnam, weights = Freq, family = cumulative(parallel = TRUE))
## Warning in eval(slot(family, "initialize")): response should be ordinal---
## see ordered()
NPO.mod = vglm(response ~ year * sex, data = Vietnam, weights = Freq, family = cumulative(parallel = FALSE))
## Warning in eval(slot(family, "initialize")): response should be ordinal---
## see ordered()
lrtest(PO.mod, NPO.mod)
## Likelihood ratio test
## 
## Model 1: response ~ year * sex
## Model 2: response ~ year * sex
##   #Df  LogLik Df Chisq Pr(>Chisq)    
## 1 114 -3878.5                        
## 2 108 -3839.7 -6 77.57  1.134e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Analysis indicates that coefficients differ for PO and NPO. It can be concluded that the intial proportional odds model does not hold for this data set. ??? (c) Fit the multinomial logistic model, also allowing an interaction. Use car::Anova () to assess the model terms.

library(nnet)
## Warning: package 'nnet' was built under R version 3.4.4
Multi.model = multinom(response ~ year * sex, data = Vietnam)
## # weights:  20 (12 variable)
## initial  value 55.451774 
## final  value 55.451774 
## converged
Anova(Multi.model)
## Analysis of Deviance Table (Type II tests)
## 
## Response: response
##          LR Chisq Df Pr(>Chisq)
## year            0  3          1
## sex             0  3          1
## year:sex        0  3          1

The multinomial logistic modelshows that there are relevant predictors based on response type.

  1. Produce an effect plot for this model and describe the nature of the interaction.
library(effects)
## Warning: package 'effects' was built under R version 3.4.4
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
plot(allEffects(Multi.model), style = "stacked", main = "Responses in relation to year and sex")

  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"))
Simp.Multi.model = multinom(response ~ sex + yearMale, data = Vietnam)
## # weights:  16 (9 variable)
## initial  value 55.451774 
## final  value 55.451774 
## converged
anova(Multi.model,Simp.Multi.model)
## Likelihood ratio tests of Multinomial Models
## 
## Response: response
##            Model Resid. df Resid. Dev   Test    Df LR stat. Pr(Chi)
## 1 sex + yearMale       111   110.9035                              
## 2     year * sex       108   110.9035 1 vs 2     3        0       1

Anova test indicate that there is not a relevant difference between the two models; since p value is larger than 0.05.