Heavy tailed residuals in linear mixed model I'm a new user of linear mixed models and I'm experiencing some troubles with that. I have a dataset with 680000 measures of milk production from 2017 to 2020, from a population of almost 37000 cows in 175 farms. I'm just trying to delucidate if there is an interaction between the thermal stress of the animals with their milk production. For that, I have a thermal stress index (THI) transformed in a factor variable with four levels (low stress, medium stress, high stress, extreme stress) and the daily milk production for each cow in Kg/day. In the following image you can see that Kg/day has a right skewed distribution.

I have also added some other important variables to my model like:

*

*number of calvings  (NP): factor variable with 5 levels,

*number of days since calving (DEL): factor, 10 levels

*season of the year (spring, summer,...) when the calving has ocurred (TP): factor, 4 levels

*number of milkings (NM) (2 daily milkings, 3 milkings, milking robot): factor, 3 levels.

I use the lmer fuction from the lme4 package in R to fit the model, with the THI, NP, DEL, TP and NM as fixed effects and cow and farm identification as random effects.
f3.3 <- lmer(KGDAY~1+THI+DEL*NP+TP+DEL*NM + (1|COW)+(1|FARM), data=(data), REML=T,control = lmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))

I have chosen this fixed effect structure by comparing the AIC of the models fitted by ML. Regarding the random part, I have attempt to use a random slope for each cow but I obtain a "boundary singular fit" (wich indicated my ramdom part is too complex, so I avoided it).
The problem begans when I analyzed the residuals structure. As you can see in the below images, I think they are heterocedastic
plot(fitted(f3.3), residuals(f3.3))
abline(h=0, lty=2)


and with heavy tailed distribution

The milk production data has a bottom limit of 4kg/day and a upper limit of 99kg/day due to the collection data routine. So I think this is the reason to the diagonal line in the
fittedVsresiduals plot.
To solve this, I have tried to transform my production data with a logaritmic, a sqrt and a Box-cox transform, but the residuals didn't vary significantly. Then I tried to change the variance and correlation structure through the nlme package but the fitting process took too much time (days) and therefore I override this option.After that, I have also attempted to use the robuslmm package but new problems arised, since my computer doen't have enough memory capacity and an error message emerges during computing (I'm using an Intel(R) Core(TM) i5-6500T, 2.50GHz, 8GB RAM and the x64bit version).
My questions are:

*

*Is it possible that the high amount of data is given my an excessive amount of significant fixed effects?

*Should I be concerned about the non-normality of the residuals, since I just want to check the effect of the THI in milk production? I will not use the model to forecast, for example.

*If yes, how can I tackle this heavy tailed distribution? Should I use a generalized mixed model or should I keep trying with the robust linear mixed model?

 A: 

*

*Is it possible that the high amount of data is given my an excessive amount of significant fixed effects?


Yes, your 680000 measurements are not independent but are made within the same cows and the same farms, and possibly also within similar small time-frames (I count on average about 18 data points per cow, but possibly some of the data points can be made close to one another and be effectively a single sample counted twice).
You do include random intercepts, but you may have all sorts of additional correlations that you did not account for and these may increase the apparent significance.
If you want an assessment of independence that is independent from the statistical model assumptions, then you might want to use part of the data as a test set to find out how well your model performs in making predictions for different cows, farms, and times.



*Should I be concerned about the non-normality of the residuals, since I just want to check the effect of the THI in milk production? I will not use the model to forecast, for example.


The non-normality is not extreme. The qq-plot indicates heavy tails but they are not of a type that makes the estimate of the standard error wrong.
There are different type of heavy tails. You have for instance the tails that relate to infinite variance, and that is not your case (since the milk production is bounded).
What you have is just a discrepancy from the standard normal approximation. This occurs possibly because the variance of the residuals is not homogeneous. But that does not invalidate the estimates of the regression and neither the estimates of the error in the regression coefficients (those coefficients are a linear sum of your observations and will follow approximately a normal distribution)



*If yes, how can I tackle this heavy tailed distribution? Should I use a generalized mixed model or should I keep trying with the robust linear mixed model?


I wouldn't be so much worried about point 2, but more about point 1.
Instead of using 10 levels of DEL I would prefer to treat DEL as a continuous variable and model it with some function (that possibly involves much less than 10 coefficients).
If you use a small data set, like one year from a few farms, to explore the data and behavior of the relationships, then you could come up with some sound model. That would reduce potential problems like clipping (a production one day later does not suddenly make a jump according to the model and create weird plots of residuals), and also you could make the model more flexible (maybe you need to include more interaction terms and more random effects). You could also use existing literature to come up with a model and reduce the amount of data that you want to use (waste) on coming up with a potential model.
Then fit coefficients based on another part of the data set, and finally measure performance with a test set.
