r anova mixed-models resampling random-effects

Nonparametric way to perform ANOVA of linear mixed model and power calcualtion

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. 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 methods. Could anyone help me with the following questions? It is really appreciated.

Is there a method and available R package to perform bootstrap anova (and post-hoc comparisons) of linear mixed model?
Is it still necessary to calculate the power of the bootstrap or nonparametric anova? If so, how to calculate the power?
I was able to use simr with methods anova to calculate power of testing group, exam, and the interaction for a lme model object. Also, can simr also be used to find power of post hoc pairwise comparisons or emmeans should be used? Thanks

Solution

This is slightly deep statistically. There are two basic approaches to resampling-based inference:

bootstrapping (parametric or nonparametric) simulates the sampling distribution of desired quantities, from which you can compute confidence intervals. (You can also compute two-tailed p-values by inverting the confidence interval, which comes out to prob_lt_0 = mean(bootsamp<0); p_val = 2*min(prob_lt_0, 1-prob_lt_0).) Post-hoc comparisons are a little more complicated because they involve not just tests of comparisons, but also some kind of multiple comparisons correction (Tukey, etc.)
permutation tests simulate the sampling distributions of test statistics under the null hypothesis, from which we can get p-values.

Parametric bootstrapping is a convenient way of resampling, especially if we have complicated grouping structures (spatial/temporal or crossed random effects) which make resampling independent groups difficult, or if we are assessing a GLMM where residual bootstrapping doesn't work. However, it assumes that the model is 'correct', which fails here because you don't want to rely on the Normality of the conditional distribution.

library(lmerTest)
library(lmeresampler)

ss <- transform(sleepstudy, oddsub = factor((as.numeric(Subject) %% 2 == 1)))
fm1 <- lmer(Reaction ~ Days + oddsub + (Days | Subject), ss)

set.seed(101)
boot1 <- reb_bootstrap(fm1,
                       .f = fixef,
                       B = 1000,
                       reb_type = 2)

## not sure why reb_bootstrap seems to ignore the .f specification?
## but we get something useful anyway

confint(boot1, type = "perc")
## 1 beta1   253.   232.    273.  perc   0.95
## 2 beta2    10.5    7.24   13.5 perc   0.95
## 3 beta3    -3.56 -32.3    26.0 perc   0.95

 
pval <- function(x) { plt0 <- mean(x<0); 2*min(plt0, 1-plt0) }
apply(boot1$replicates[,1:3], 2, pval)
## beta1 beta2 beta3 
## 0.000 0.000 0.802

You can do a permutation test, but you have to do the test one term at a time, by comparing nested models (e.g. the full model vs. a model without the effect of Days):

library(predictmeans)
fm0 <- update(fm1, . ~ . - Days)
p1 <- permlmer(fm0, fm1)

Data: ss
Models:
lmer0: Reaction ~ oddsub + (Days | Subject)
lmer1: Reaction ~ Days + oddsub + (Days | Subject)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq) Perm-p
lmer0    6 1787.4 1806.6 -887.70   1775.4                            
lmer1    7 1765.9 1788.2 -875.94   1751.9 23.537  1 1.2256e-06  0.003

For power, I would probably just do the analysis assuming normality and then treat it as a slightly optimistic estimate; power analyses are always guesses/approximations anyway.

PS part of the reason that people are suggesting other options (robust LMM, ordinal models, etc.) is that these are kind of a nuisance.