Can you use TukeyHSD as a post hoc test for a mixed model?

I completed a mixed model with repeated measures. I did an ANOVA comparing my model to a null model and got a significant value. Now I need to identify where those significant differences are in my data. What post hoc test can I use? Will Tukey HSD work? Many of the other test suggested such as wald.test and 1-pchisq don't work for me. If you need more info on the tests I ran let me know. Thanks!

• Can you elaborate on your situation & your analyses (etc)? What do you mean they didn't "work for you"? – gung - Reinstate Monica Feb 23 '16 at 19:24
• The answer is yes, provided you have a nice balanced design and you use the appropriate standard error of differences derived from the covariance structure of the mixed model. But if by this you mean to compute the HSD via the formula in Ch 2 or so of your design book, the answer is a definite no. If you're using R, the one-step (default) test in the multcomp package computes the equivalent of Tukey-adjusted p values, even in some messier cases; as does the lsmeans package with adjust="mvt" – rvl Feb 23 '16 at 20:00
• @rvl, can you elaborate on that & turn it into an official answer? I think that would be a nice contribution. – gung - Reinstate Monica Feb 23 '16 at 20:29
• I'll try. Won't add the R-specific stuff unless it's needed. – rvl Feb 23 '16 at 22:21

Consider a balanced design under the simplest assumptions, and let $\bar y_{i\cdot\cdot}$ denote the mean at the $i$th level of one of the fixed factors. With these utopian assumptions, the $\bar y_{i\cdot\cdot}$ are either independent, multivariate normal with a compound-symmetric covariance structure whereby any pair of them has the same correlation $\rho$. It can then be shown that the set of all pairwise pairwise comparisons among the $\bar y_{i\cdot\cdot}$ has the same correlation matrix regardless of the value of $\rho$; and hence that the distribution of the Studentized range holds. Accordingly, the Tukey HSD method is valid as long as the correct standard error estimate is used. The HSD at level $\alpha$ is given by $$HSD_{ij} = \frac{q_\alpha(k,\nu)}{\sqrt 2}\cdot SE(\bar y_{i\cdot\cdot} - \bar y_{j\cdot\cdot})$$ where $q_\alpha(k,\nu)$ is the critical value for the Studentized range, $k$ is the number of different levels $i$, and $\nu$ is the appropriate degrees of freedom for estimating the $SE$. (The tabled value is divided by $\sqrt 2$ because the Studentized range is defined in terms of the SE of one mean rather than the SE of a difference of two of them.) The HSD method protects the overall significance level of the family of comparisons by assuring that in the worst case (comparing the maximum with the minimum), the probability of a type I error is $\alpha$.
In a more complicated situation such as unbalanced data and fancy covariance structures, the Studentized-range tables no longer apply because the joint distribution of the comparisons is more complicated. However, with appropriate software, an approximate critical value for for $\max(\bar y_{i\cdot\cdot} - \bar y_{j\cdot\cdot})$ can be obtained using the multivariate $t$ distribution with the covariance structure found in the data; then that may be used in place of $\frac{q_\alpha(k,\nu)}{\sqrt 2}$ in the formula.