What to do with a lmer model that violates the assumptions of normal residuals?

I have the following lmer model:

    ball6lmer=lmer(Width_mm~Buried+(Buried|Chamber), data=ballData,control=lmerControl(optimizer="optimx",optCtrl=list(method="nlminb")), REML=FALSE)


Width_mm is the diameter of dung balls made by dung beetles Buried is whether or not a dung ball was buried by a beetle (0=no, 1= yes). 26 chambers were included in the experiment and multiple dung balls came from the same chamber, so I have included the random term ~1|chamber to account for the non-independence of data. Furthermore, log likelihood tests showed that the model was improved by allowing the slopes of the Buried term to vary over chambers (random=Buried|Chamber). The output looks like this:

    Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Width_mm ~ Buried + (Buried | Chamber)
Data: ballData
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

AIC      BIC   logLik deviance df.resid
668.4    689.9   -328.2    656.4      260

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.7717 -0.3993 -0.1455  1.0212  2.4765

Random effects:
Groups   Name        Variance Std.Dev. Corr
Chamber  (Intercept) 0.009828 0.09914
Buried1     0.182687 0.42742  1.00
Residual             0.625589 0.79094
Number of obs: 266, groups:  Chamber, 26

Fixed effects:
Estimate Std. Error t value
(Intercept) 11.14463    0.06689  166.60
Buried1     -0.92226    0.13534   -6.81

Correlation of Fixed Effects:
(Intr)
Buried1 -0.270
>


My residual plots do not look so good though:

The residuals are clearly not normal. I am now stuck as to how to proceed with this analysis. I could do a simple Mann–Whitney U test to test for a difference between Buried, however this would ignore the non-independence of dung balls from the same chamber.

The next option that I am considering is to bootstrap the regression coefficients, and obtain 95% confidence intervals for the regression coefficients. However I have never done this before and would not know where to start. So I am looking for a bit of guidance on how to proceed. Specifically:

1) Would a simple Mann–Whitney U test be an appropriate alternative or is it incorrect to ignore the non-independence of the data?

2) Is trying to bootstrap this model something that a novice like me could try or is it something best left to the professionals? If it is a relatively straight forward thing to do, could anyone recommend any basic resources to get me started?

3) Are there any other options for dealing with lmer models that violate the assumptions of normal residuals? My response variable only has 5 possible variables (width of 9,10,11,12 or 13mm) so I don’t think transforming the response variable will help.

Any guidance would be much appreciated. Thanks.

AN UPDATE:

John as correctly pointed out that I have been using the wrong predictor variable. I have now respecified the model to predict ball width based on whether or not an egg was laid in a dung ball (EggLaid, a 2 level predictor) The model is now:

    lmer(Width_mm~EggLaid+(1|Chamber), data=ballData, REML=FALSE)


With the following output:

    Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Width_mm ~ EggLaid + (1 | Chamber)
Data: ballData

AIC      BIC   logLik deviance df.resid
712.8    727.1   -352.4    704.8      262

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.35659 -0.57566 -0.06791  0.79089  2.62157

Random effects:
Groups   Name        Variance Std.Dev.
Chamber  (Intercept) 0.1059   0.3254
Residual             0.7636   0.8739
Number of obs: 266, groups:  Chamber, 26

Fixed effects:
Estimate Std. Error t value
(Intercept)  11.0697     0.1078  102.65
EggLaid1     -0.5620     0.1119   -5.02

Correlation of Fixed Effects:
(Intr)
EggLaid1 -0.593


The QQ plot of residuals look much better, however I still have an odd looking residual vs fitted plot, which I assume reflects the 5 values of the response variable:

Is this residual plot "passable" given the 5 values of the reponse variable, or is this model not OK and I need to bootstrap the coefficients? Also to clarify, the values for width (the response) are rounded to the nearest mm.

A FURTHER UPDATE: I have had a go at bootstrapping the regression coefficients (using parametric bootsrap) using some helpful info here. My code and output is as follows:

    ball2=lmer(Width_mm~EggLaid+(1|Chamber), data=ballData, REML=FALSE)
#This is the same model as the one above (output is above)

#With 5000 bootstrapped replicates
b_par<-bootMer(x=ball2,FUN=fixef,nsim=5000, use.u = FALSE, type="parametric")

#Bootstrapped CI for the intercept
boot.ci(b_par,type="perc",index=1)

#Bootstrapped CI for the EggLaid1
boot.ci(b_par,type="perc",index=2)

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates

CALL :
boot.ci(boot.out = b_par, type = "perc", index = 1)

Intervals :
Level     Percentile
95%   (10.86, 11.28 )
Calculations and Intervals on Original Scale

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates

CALL :
boot.ci(boot.out = b_par, type = "perc", index = 2)

Intervals :
Level     Percentile
95%   (-0.7781, -0.3469 )
Calculations and Intervals on Original Scale


I am a complete newbie to bootstrapping, so if I have made an error in the code, any suggestions would be graciously welcomed

• Why are you predicting width with burial? Simply based on order wouldn't one expect the burial to be dependent in width and not the other way around? But if you really want to predict width with burial are you sure there aren't some measurement issues? You have a large gap in your fitted values and it seems like there are 5 groupings of measures. I'd check my data coding at least. – John Oct 24 '16 at 2:36
• @John you raise a very good point. The female buries a ball, and then hollows out a chamber to lay her egg into the ball. The hollowing out by the female (i.e. egg laying) changes the size of the ball. I have updated the analysis accordingly. The QQplot looks alot better, but the residuals v predicted are still a bit squiffy (see above). To clarify, there are only 5 values of the reponse variable because the measurements were rounded to the nearest mm – JeanDrayton Oct 24 '16 at 4:46

The funny residuals almost certainly relate to the quantization of the response to the set $\left\{9, 10, 11, 12, 13 \right\}$. A normal distribution, conditional on covariates, can't do that--the event has probability zero. As you've said in your update, this quantization is due to rounding (ie, the actual values are not exactly $9.0$ or $13.0$ but were just rounded to those values).

I think this is a bit of a case of "good news/bad news." On the one hand, if the outcome has been rounded, that's really a complicated form of interval censoring, so non-parametric inference won't help. I hypothesize that the point estimates will be at least somewhat biased due to the censoring, and you would need to use a method that models this in order to reduce this bias. FWIW, I don't see any readily available software that does this for mixed models. The R package grouped fits these models for vanilla GLMs.

On the other hand, the rounding may not be a huge issue. In the first residuals vs fitted plot in which you regress Buried there seems there was a tendency to avoid prediction of non-integer values, which makes me worry that the model could be optimistic with its standard error estimates due to the quantization. This doesn't seem to be so true in the second model where EggLaid was the predictor. It might be worth running a simulation study to bound the bias, because I don't have a great intuition about how small or big this bias might be.

In conclusion, I would say yes, you should bootstrap, and that ideally you should use a semi-parametric bootstrap that jitters the response variables with a bit of noise. It looks like this happen already if you use the function bootMer in lme4. From the docs:

If ‘use.u’ is ‘FALSE’ and ‘type’ is ‘"parametric"’, each simulation generates new values of both the “spherical” random effects u and the i.i.d. errors epsilon, using ‘rnorm()’ with parameters corresponding to the fitted model ‘x’.

These iid errors epsilon will mean that the bootstrapped sample is no longer quantized. So as long the the initial residual variance is approximately correct, the bootstrap will do something reasonable.

• The width measurements were measured to the nearest mm, so the quantization is due to rounding (for example a ball could have been 10.3mm, but was recorded as 10mm). I have updated the residual plots after realising that I used the wrong predictor. The updated plot is above. Based on this plot, would you recommend bootstrapping? Or can I use the current model? – JeanDrayton Oct 24 '16 at 5:01
• Hi Andrew, I have had a crack at the semi parametric bootstrapping as you suggested and my code is above. Could you please clarify what you mean by "jitters the response variables with a bit of noise. The 1 mm rounding suggest that you should add at least .25 mm since that is within the measurement limits"? Thank you for your suggestions so far, they have been very useful. – JeanDrayton Oct 24 '16 at 6:47
• I took a closer look at the lmer docs and the "parametric" bootstrap actually nearly does what I suggested regarding the jittering. I also revised the answer to reflect a bit more of my uncertainty about exactly how big of an effect the rounding could have. If I were you, I would run a simulation to better understand the impact. – Andrew M Oct 24 '16 at 17:43
• Thanks Andrew. I have altered the code above to perform a parametric bootstrap. Thank you for your time to look over my code and offer your thoughts, it was very helpful and much appreciated. – JeanDrayton Oct 28 '16 at 22:47
• The Bayesian R package brms has a brm function that can include censoring. This might address some of your comments. Of course, it's slower than lmer -- MCMC sampling and all -- but more flexible as well. – Wayne Oct 28 '16 at 23:25