# Fit the proportional odds model, p.179, in McCullagh and Nelder (1989)
data(pneumo)
pneumo = transform(pneumo, let=log(exposure.time))
(fit = vglm(cbind(normal, mild, severe) ~ let,
cumulative(parallel=TRUE, reverse=TRUE), pneumo))
fit@y # Sample proportions
weights(fit, type="prior") # Number of observations
coef(fit, matrix=TRUE)
constraints(fit) # Constraint matrices
# Check that the model is linear in let
fit2 = vgam(cbind(normal, mild, severe) ~ s(let, df=2),
cumulative(reverse=TRUE), pneumo)
plot(fit2, se=TRUE, overlay=TRUE, lcol=1:2, scol=1:2)
# Check the proportional odds assumption with a likelihood ratio test
(fit3 = vglm(cbind(normal, mild, severe) ~ let,
cumulative(parallel=FALSE, reverse=TRUE), pneumo))
1 - pchisq(2*(logLik(fit3)-logLik(fit)),
df=length(coef(fit3))-length(coef(fit)))
Run the code above in your browser using DataLab