I have a small data where there are 3 groups (A,B,C) and 5 participants from each group. All of those participants are measured 6 times on each of 7 different exams, so each participant get 6*7=42 scores in total. A simple linear mixed model was built mylmm<-lmer(score ~1+group+exam+group*exam+(1|participant), data = mydata). I could obtain the anova results and post-hoc pairwise comparison for group, exam, and interaction of them using anova(mylmm) and multiple comparison function.

However, the data is very small (only 5 participants) and residual of mylmm is not normal, so the power is insufficient. The score is continuous and I really need to keep it in the original scale if possible (but tranformations including box-cox don't work). I am aware of robust mixed model using robustlmm and residual bootstrap mixed model using lmeresampler. However, I am unable to perform anova and multiple comparisons using these packages. Could anyone help me with the following questions? It is really appreciated.

  1. Is there a method and available R package to perform bootstrap anova (and post-hoc comparisons) of linear mixed model?
  2. Is it still necessary to calculate the power of the bootstrap or nonparametric anova? If so, how to calculate the power?
  3. I am also aware of aligned ranks transformation ANOVA. Does this method work here and especially for small data? Also, how can I calculate for the power?
  4. There is also another related but not so important question. For the model built directly from lme (without bootstrapping), I was able to use simr with methods anova to calculate power of testing group, exam, and the interaction. Can simr also be used to find power of post hoc pairwise comparisons? Thanks.
  • $\begingroup$ When you say the residuals are not normal, are you referring to the within-participant residuals? Can you show us the graph from lattice::qqmath(mylmm)? Also see this handy thread about residuals. stats.stackexchange.com/questions/111010/… $\endgroup$
    – Erik Ruzek
    Jan 3 at 20:51

1 Answer 1


Aligned-rank nonparametric procedures are essentially ad hoc. I suggest using a random effects semiparametric ordinal regression model such as the proportional odds model. The R ordinal package will do this. To accurately handle the random effects take care to not use the default numerical integration procedure, as briefly discussed here. See here for resources for ordinal regression. The proportional odds model generalizes the Wilcoxon and Kruskal-Wallis tests.

Note that a Bayesian proportional odds model will provide even more accurate inference due to not using numerical integration approximations.

  • $\begingroup$ Thank you @FrankHarrell. Sorry it wasn't so clear. The score is not ordinal but continuous and I need to interpret the results in the original scale. $\endgroup$
    – ksing
    Jan 3 at 14:30
  • $\begingroup$ With ordinal regression you can restate the results on any scale, e.g., get estimated medians and means. Examples are here. Ordinal regression works great for continuous Y. $\endgroup$ Jan 4 at 7:43

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.