Assumptions of linear mixed model not met I have a repeated measures dataset with which I'm testing if individuals are consistent in their boldness scores (continuous variable) over time (trials). Towards this, I generated linear mixed effects model with boldness scores as the dependent variable, trials as fixed effects and individual IDs as random effects. However, when I checked the residuals of the fitted model, the assumption of normality is not met. I have attached the qqplot of the residuals.
Before transformation
 
I then transformed the boldness scores by boxcox method and again ran the model. The residuals seem to appear more normal now when compared to the residuals from the model with the non- transformed variable. However, even after transformation, shapiro tests show that the residuals are not normally distributed.
After transformation
 
The skeweness of the residuals has also reduced from +1.111 to -0.24 after transformation.
Can I proceed with the model ignoring the assumption of normality of residuals? (see blogpost here which suggests linear mixed models are robust even if this assumption is not met).  
 A: The untransformed residuals appear to be bimodal. If so, this may indicate that there are two clusters within the dataset, and a transformation is very unlikely to result in a well-fitting model, and if it did there is a high chance that it would be over-fitted.
It is evident that the Box-Cox-transformed model residuals are not plausibly normal, as we would expect if the data are bimodal/clustered.
I would advise caution before proceeding.  Typically this situation arises because the model is mis-specified. Some possible reasons include:  


*

*there could be one or more missing predictor variables that
"explain" the bimodal structure (that is, they explain some
separation in the data that the bimodality).  

*there could be one or more missing higher order terms of existing explanatory variables that "explain" the bimodal structure.

*the underlying data generation model may not be approximated by a linear model in the range of the data observed


Possible solutions are:  


*

*Add other explanatory variables, if available, ideally informed by
theory of the data generation process

*Again informed by the underlying theory, try adding higher order
terms.  

*Consider performing a cluster analysis or latent class analysis to
identify clusters and introduce a new dummy variable which identifies
which cluster an observation belongs to.  

*Consider a non-linear model, again informed by the underlying theory

*Consider splitting the data so that one cluster is in one dataset and the other is in the other.
