# NOT RUN {
set.seed(1)
pneumo <- transform(pneumo, let = log(exposure.time),
x3 = runif(nrow(pneumo)))
fit1 <- vglm(cbind(normal, mild, severe) ~ let , propodds, pneumo)
fit2 <- vglm(cbind(normal, mild, severe) ~ let + x3, propodds, pneumo)
fit3 <- vglm(cbind(normal, mild, severe) ~ let , cumulative, pneumo)
# Various equivalent specifications of the LR test for testing x3
(ans1 <- lrtest(fit2, fit1))
ans2 <- lrtest(fit2, 2)
ans3 <- lrtest(fit2, "x3")
ans4 <- lrtest(fit2, . ~ . - x3)
c(all.equal(ans1, ans2), all.equal(ans1, ans3), all.equal(ans1, ans4))
# Doing it manually
(testStatistic <- 2 * (logLik(fit2) - logLik(fit1)))
(pval <- pchisq(testStatistic, df = df.residual(fit1) - df.residual(fit2),
lower.tail = FALSE))
(ans4 <- lrtest(fit3, fit1)) # Test PO (parallelism) assumption
# }
Run the code above in your browser using DataLab