Robustness of linear mixed-effects model to violations of multuvariate normality

Linear mixed-effects models (LMMs) are widely used to analyze nested longitudinal data by simultaneously modeling within-person dynamics and between-person differences. LMMs assume that random effects follow a multivariate normal (MVN) distribution, an assumption that could be suspicious in practice. I conducted a Monte Carlo simulation study to evaluate the robustness of LMMs when true random effects do not follow an MVN distribution. Data were generated under a two-level longitudinal design in which a linear fixed effect of time was of primary inferential interest, crossing two sample-size configurations (3 occasions × 50 subjects; 10 × 300) with two residual structures (white noise; AR(1)). Following Sklar's theorem, random effects were generated from joint distributions defined by copulas (Gaussian, t) and variance-matched marginals (normal, t, triangular), isolating distributional shape from location and scale. Each dataset was analyzed with four models using REML estimation: ordinary regression, random-intercept-only, random-intercept-and-slope with white noise residual, and random-intercept-and-slope with AR(1) residuals. Fixed-effect estimates remained essentially unbiased under all violations, with near-nominal coverage whenever a random slope was specified, where instead of non-normality, omitting the random slope was the primary threat to inference (coverage ≈ 50%). In contrast, heavy-tailed marginals produced substantial underestimation of random-effect standard deviations (relative bias ≈ −20%) regardless of copula or sample size. The AR(1) parameter was unrecoverable with only three occasions, going undetected by AIC/BIC. We discuss implications for model specification and for variance-component-focused applications such as hierarchical meta-analysis.