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.2517     0.5464   -2.29
## TreatmentTreated  1.7453     0.4759    3.67
## Age               0.0382     0.0184    2.07
## 
## Intercepts:
##             Value  Std. Error t value
## None|Some    2.532  1.057      2.395 
## Some|Marked  3.431  1.091      3.144 
## 
## Residual Deviance: 145.46 
## AIC: 155.46

Anova(arth.polr)      # Type II tests

## Analysis of Deviance Table (Type II tests)
## 
## Response: Improved
##           LR Chisq Df Pr(>Chisq)    
## Sex           5.69  1    0.01708 *  
## Treatment    14.71  1    0.00013 ***
## Age           4.57  1    0.03251 *  
## ---
## 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))