R Markdown

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.

Model below shows that all terms are significant, including the interaction.

library(vcdExtra)
## Loading required package: vcd
## Loading required package: grid
## Loading required package: gnm
data("Vietnam", package="vcdExtra")
str(Vietnam)
## 'data.frame':    40 obs. of  4 variables:
##  $ sex     : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
##  $ year    : int  1 1 1 1 2 2 2 2 3 3 ...
##  $ response: Factor w/ 4 levels "A","B","C","D": 1 2 3 4 1 2 3 4 1 2 ...
##  $ Freq    : int  13 19 40 5 5 9 33 3 22 29 ...
library(MASS)
Vietnam$response = ordered(Vietnam$response)
Vietnam_polr = polr(response ~ year * sex, data = Vietnam, weights = Freq, Hess = TRUE)
summary(Vietnam_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
library(car)
## Loading required package: carData
## 
## Attaching package: 'carData'
## The following object is masked from 'package:vcdExtra':
## 
##     Burt
Anova(Vietnam_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
  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. *Since the model are not the same they aren’t proportional

Comparison shows that Vietnam_NPO fits better.

library(VGAM)
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:car':
## 
##     logit
Vietnam_PO = vglm(response ~ year * sex, data = Vietnam, weights = Freq, family = cumulative(parallel = TRUE))
Vietnam_NPO = vglm(response ~ year * sex, data = Vietnam, weights = Freq, family = cumulative(parallel = FALSE))
lrtest(Vietnam_NPO, Vietnam_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

(c)Fit the multinomial logistic model, also allowing an interaction. Use car::Anova () to assess the model terms.

library(nnet)
(Vietnam_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
## 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
## 
## Residual Deviance: 7677.71 
## AIC: 7701.71
Anova(Vietnam_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

(d)Produce an effect plot for this model and describe the nature of the interaction.

library(effects)
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
plot(allEffects(Vietnam_polr), style = "stacked")

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.

Anova result shows that there’s no significant difference between the two models because p value is larger than 0.05.

Vietnam = within(Vietnam, yearMale <- year * (sex == "Male"))
Vietnam_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(Vietnam_mtn, Vietnam_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