8.2

library(knitr)
library(MASS)
library(car)
## Loading required package: carData
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:car':
## 
##     logit
library(rms)
## Loading required package: Hmisc
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
## 
## Attaching package: 'rms'
## The following objects are masked from 'package:VGAM':
## 
##     calibrate, lrtest
## The following objects are masked from 'package:car':
## 
##     Predict, vif
library(reshape2)
library(ggplot2)
library(directlabels)
library(effects)
## Use the command
##     lattice::trellis.par.set(effectsTheme())
##   to customize lattice options for effects plots.
## See ?effectsTheme for details.
data("housing", package="MASS")
head(housing$Sat)
## [1] Low    Medium High   Low    Medium High  
## Levels: Low < Medium < High

a

PO.mod <- polr(Sat ~ Type + Infl + Cont, data = housing, weights = Freq)
summary(PO.mod)
## 
## Re-fitting to get Hessian
## Call:
## polr(formula = Sat ~ Type + Infl + Cont, data = housing, weights = Freq)
## 
## Coefficients:
##                 Value Std. Error t value
## TypeApartment -0.5724    0.11924  -4.800
## TypeAtrium    -0.3662    0.15517  -2.360
## TypeTerrace   -1.0910    0.15149  -7.202
## InflMedium     0.5664    0.10465   5.412
## InflHigh       1.2888    0.12716  10.136
## ContHigh       0.3603    0.09554   3.771
## 
## Intercepts:
##             Value   Std. Error t value
## Low|Medium  -0.4961  0.1248    -3.9739
## Medium|High  0.6907  0.1255     5.5049
## 
## Residual Deviance: 3479.149 
## AIC: 3495.149

Log odds coefficient can be interpreted as: for 1 unit inc in the predictor, the response variable is expected to change by coefficient while the other variables in the model are held constant

b

Anova(PO.mod)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Sat
##      LR Chisq Df Pr(>Chisq)    
## Type   55.910  3  4.391e-12 ***
## Infl  108.239  2  < 2.2e-16 ***
## Cont   14.306  1  0.0001554 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

All predictor variables are significant in the model since p-values are all < 0.001 for response variable Satisfaction

c

Assuming odds are proportional

houspo <- vglm(Sat ~ Type + Infl + Cont, data = housing, weights = Freq, family=cumulative(parallel = TRUE))
houspo
## 
## Call:
## vglm(formula = Sat ~ Type + Infl + Cont, family = cumulative(parallel = TRUE), 
##     data = housing, weights = Freq)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2 TypeApartment    TypeAtrium   TypeTerrace 
##    -0.4961346     0.6907087     0.5723499     0.3661866     1.0910147 
##    InflMedium      InflHigh      ContHigh 
##    -0.5663937    -1.2888184    -0.3602845 
## 
## Degrees of Freedom: 144 Total; 136 Residual
## Residual deviance: 3479.149 
## Log-likelihood: -1739.575

Assuming odds are not parallel

housnpo <- vglm(Sat ~ Type + Infl + Cont, data = housing, weights = Freq, family=cumulative(parallel = FALSE))
housnpo
## 
## Call:
## vglm(formula = Sat ~ Type + Infl + Cont, family = cumulative(parallel = FALSE), 
##     data = housing, weights = Freq)
## 
## 
## Coefficients:
##   (Intercept):1   (Intercept):2 TypeApartment:1 TypeApartment:2 
##      -0.4461683       0.6466011       0.6011532       0.5379260 
##    TypeAtrium:1    TypeAtrium:2   TypeTerrace:1   TypeTerrace:2 
##       0.1911686       0.4830450       1.0790012       1.1189824 
##    InflMedium:1    InflMedium:2      InflHigh:1      InflHigh:2 
##      -0.5924647      -0.5497942      -1.2191597      -1.3077420 
##      ContHigh:1      ContHigh:2 
##      -0.4304931      -0.2956644 
## 
## Degrees of Freedom: 144 Total; 130 Residual
## Residual deviance: 3470.579 
## Log-likelihood: -1735.289

Comparison

coef(houspo, matrix=TRUE)
##               logit(P[Y<=1]) logit(P[Y<=2])
## (Intercept)       -0.4961346      0.6907087
## TypeApartment      0.5723499      0.5723499
## TypeAtrium         0.3661866      0.3661866
## TypeTerrace        1.0910147      1.0910147
## InflMedium        -0.5663937     -0.5663937
## InflHigh          -1.2888184     -1.2888184
## ContHigh          -0.3602845     -0.3602845
coef(housnpo, matrix=TRUE)
##               logit(P[Y<=1]) logit(P[Y<=2])
## (Intercept)       -0.4461683      0.6466011
## TypeApartment      0.6011532      0.5379260
## TypeAtrium         0.1911686      0.4830450
## TypeTerrace        1.0790012      1.1189824
## InflMedium        -0.5924647     -0.5497942
## InflHigh          -1.2191597     -1.3077420
## ContHigh          -0.4304931     -0.2956644
op <- par(mfrow=c(1,3))
plot.xmean.ordinaly(Sat ~ Type + Infl + Cont, data = housing, lwd= 2, pch= 16, subn= FALSE)

The proportional odds assumption assumes the slope of the line is fixed and does not change for each new value. If the slope changes for each new value, then it is NPO. The difference in slopes between the 2 coefficient columns is negligible in NPO. Hence the condition for proportional odd holds

d

plot(Effect(c("Infl","Type"), PO.mod))
## 
## Re-fitting to get Hessian

Observing the plots, it can said that when Infl is low, chance of Sat being low is high and when Infl is high, chance of Sat being high is also high.