Arthritis data: Proportional odds model

Load packages and data

library(vcd)
library(car)         # for Anova()
library(MASS)        # for polr()
library(effects)     # effect plots

data(Arthritis, package = "vcd")

Fit the proportional odds model

arth.polr <- polr(Improved ~ Sex + Treatment + Age, data=Arthritis)
summary(arth.polr)

## Call:
## polr(formula = Improved ~ Sex + Treatment + Age, data = Arthritis)
## 
## Coefficients:
##                     Value Std. Error t value
## SexMale          -1.25168    0.54636  -2.291
## TreatmentTreated  1.74529    0.47589   3.667
## Age               0.03816    0.01842   2.072
## 
## Intercepts:
##             Value   Std. Error t value
## None|Some    2.5319  1.0571     2.3952
## Some|Marked  3.4309  1.0912     3.1443
## 
## Residual Deviance: 145.4579 
## AIC: 155.4579

Anova(arth.polr)      # Type II tests

## Analysis of Deviance Table (Type II tests)
## 
## Response: Improved
##           LR Chisq Df Pr(>Chisq)    
## Sex         5.6880  1  0.0170812 *  
## Treatment  14.7095  1  0.0001254 ***
## Age         4.5715  1  0.0325081 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arth.effects <- allEffects(arth.polr, xlevels=list(age=seq(15,45,5)) )
plot(arth.effects)  # visualize all main effects

plot terms not in model

plot(effect("Sex:Age", arth.polr))

plot(effect("Treatment:Age", arth.polr))