I have a logistic regression model which has one (significant) explanatory term and seven levels. I am trying to work out how to test for significant differences between the levels. It seems I have two options:

1) Looking at the p-values associated with each of the level coefficents, and re-running the model seven times changing the reference level each time. I can then construct a table of all pairwise comparisons. As these wouldn't really be separate tests (as I'm only changing the reference level, not the model), I don't think I'd need to correct for multiple comparisons?

2) A Tukey HSD i.e. summary(glht(sp2model, mcp(factor="Tukey")))

I'd be grateful if people could advise on whether option 1 is correct and/or what seems reasonable. I am working in R.


1 Answer 1


Let's say your model includes only the binary (?) response and a single predictor variable with 7 levels, so that you are modelling the log odds of your binary response taking the value 1 (rather than 0) as a function of the predictor variable.

If you have no a priori contrasts (or comparisons) that you are interested in testing based on your model, then it makes sense to consider all pairwise comparisons of levels of the predictor variable so that you can compare the log odds of the binary response variable taking the value 1 (rather than 0) for pairs of levels of the predictor variable.

If you choose your option 1), the p-values associated with the contrasts of log odds among pairs of levels will not be adjusted for multiplicity, even though they should. The multiplicity arises from performing multiple comparisons using the same data set. As per http://home.cc.umanitoba.ca/~kesel/bookchap.pdf, "when multiple tests of significance are computed, the probability that at least one will be significant by chance alone increases with the number of tests examined."

If you choose your option 2), the p-values can be adjusted for multiplicity. For example:

model.glht <- glht(model, linfct = mcp(factor = "Tukey"))     
summary(model.glht, test = adjusted("Westfall"))

Try summary(model.glht) as well and check whether there is any mention of single-step adjustment for the p-values reported in its output.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.