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!
 A: 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 the above formula, the SE is based on a combination of one or more variance components, depending on which ones are involved in the determining the variance of the difference. A between-subjects comparison involves more variances than a within-subjects one. You can't just use the formula shown in a textbook in the context of a one-way ANOVA. Some software can output these standard errors; in other software, it isn't easy to get.
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.
