Manuscript: Robustness of linear mixed-effects models to violations of distributional assumptions
Review by: Henrik Singmann

Summary: The authors perform a large set of simulation studies for testing the robustness of linear and some selected generalized linear mixed models when a number of model assumptions are violated. Results show that estimates from mixed models are relatively robust with respect to the violations introduced in the simulations.

Evaluation: I am not aware of a similar simulation study and thus find the results extremely interesting and important. I have read the manuscript with great interest and believe many researchers that want to utilize such models will do so similarly. My only real issue with the manuscript is that I wished it was longer and contained more simulations and/or more results. This does not mean I do not think this manuscript is ready for publication as is, I believe it is. But I hope the authors will extend their work in the future to encompass more situations.
In particular, she simulated situations only concern random intercept models and not the (at least in my field) more common random slopes models. This is acknowledged in the discussion, but this reservation should also be mentioned in the last paragraph (e.g., "We conclude that mixed effects models are largely robust even to quite severe violations of model assumptions at least in cases in which the model only contains random intercepts.").
Furthermore, the only thing I feel is missing from the results is what the effect of the violation is on the estimated degrees of freedoms from either the Satterthwaite method (as implemented in lmerTest) or the Kenward-Roger method (as implemented in pbkrtest).
As said above, I think neither of this points should prevent the publication in the present form as the manuscript already contains a very large and interesting set of results and is written in a clear and concise manner.

In addition to these general comments, I have a few specific comments I hope the authors can address in a revision.

• I have problems understanding the logic behind and expression for restoring the variance in the case of the bimodal and heteroskedastic distribution. In particular, I do not fully understand how exactly the variability in the square root terms (i.e., the terms set to 0.129 and 0.064, respectively) is added. Is this done after simulating fitting the simulated data sets on top of the model estimates (i.e., the variability is not actually part of the simulation and model fitting)? I am sure I am misunderstanding what was actually done here, but if the noise were not part of the simulation but added post-hoc, I would not understand how these results are actually informative to the question at hand. A better way to me seems to scale either the baseline or group variances based on the theoretical variances for the data generating distributions, at least in the bimodal case. For example, there is an analytical expression for the variance of a bimodal distribution that is a mixture of two normals, see: https://stats.stackexchange.com/a/354566/442
  According to this formula, in your case the theoretical true variance should be 2.75, which also seems to hold via simulation:
  nsim <- 1e6
  var(c(rnorm(nsim/2, 1.5, sqrt(0.5)), rnorm(nsim/2, -1.5, sqrt(0.5))))
• p. 10: I find it difficult to understand the data generating model for the missing random effects section. Either already refer to the graphical representation or make the nesting more explicit via another formal definition. For example, maybe this can be written in mathematical terms.
• p. 11, l. 226: "Balanced sampling with data missing at random:" The explanation of the sampling plan does not describe any missing data and thus seems incomplete.
• I think it would be preferable if the figure caption where in proximity to the figures in the supplemental material and not only at the beginning (this point is of course moot, if the supplemental material will be typeset by the journal as well).
• The text uses the term "simulation set" whereas the supplemental material and figures use the term "scenario". Consistency would be preferable.
• l. 151 & 152: Hat symbol above sigmas is distorted in the PDF.
• Table 1: I find the definition of the term "prediction error" unclear. It is: "Square root of the mean squared difference between the estimated and the true (simulated) value." The parenthetical '(simulated)' seems to imply that the true value changes in each simulation (as it is simulated). However, as far as I understand it, the true value is fixed for all simulations (i.e., the value represented by the red triangle in for example Figure 1).
• Table 1: BLUPs are also known as conditional modes in lme4 terminology, see: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#confidence-intervals-on-conditional-meansblupsrandom-effects
• l. 259 - 261: "The only notable exception was bias in the estimate of the group variance when the underlying distribution was bimodal, which resulted in slight upward biased (Figure S2)." Two issues: (1) The figure with bimodal distribution is S1. (2) I found it very difficult to see this "slight upward bias". As far as I can tell it refers to Scenario B3 (row 3) and column 5. Maybe this panel can also be printed in dark to make this clear.
• l. 322: "though coverage of the confidence intervals with somewhat too low", the word 'with' seems wrong here.