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