Heavy-tailed errors in mixed-effects model I'm relatively new to statistical modelling and `R', so please let me know If I should provide any further information/plots. I did originally post this question here, but unfortunately have not received any responses yet.
I am using the lme() function from nlme in R to test the significance of fixed
effects of a repeated measures design. My experiment involves subjects listening to a pair of sounds and adjusting the sound level (in decibels) until both are equally loud. This is done for a 40 different pairs of stimuli, with both orders tested (A/B and B/A). There are a
total of 160 observations per subject (this includes 1 replication for every condition), and 14 participants in all.
The model:
LevelDifference ~ pairOfSounds*order, random = (1|Subject/pairOfSounds/order)
I have built the model up by AIC/BIC and likelihood ratio tests (method =
"ML"). Residual plots for the within-group errors are shown below:

The top left plot shows standardised residuals vs fitted values. I don't see any
systematic pattern in the residuals, so I assume that the constant variation assumption
is valid, although further inspection of the subject-by-subject residuals do
show some unevenness. In conjunction with the top right plot, I have no reason to suspect
non-linearities. 
My main concern is the lower left qqplot which reveals that the residuals are
heavy-tailed. I'm not sure where to go from here. From reading Pinheiro and Bates
(2000, p. 180), the
fixed-effects tests tend to be more conservative when the tails are
symmetrically distributed. So perhaps I'm OK if the p-values are very low?
The level two and three random effects show a similar departure from normality.
Basically:


*

*How do such heavy-tailed residuals effect the inference for fixed-effects?

*Would robust regression be an appropriate alternative to check?

 A: Looking at models based on the t-distribution is potentially helpful as others wrote. However one reason to use the Gaussian assumption is that the Gaussian distribution minimises the Fisher-Information for given variance. This means that Gaussian parameters cannot be as precisely estimated as parameters of other distributions given the knowledge of the distribution. This is a good thing in the sense that using the Gaussian distribution is normally a conservative choice, as is also shown in the answer by Eoin. The larger uncertainty indicated by the Gaussian model may actually be realistic given that we don't know what the true distribution is.
Note that this argument technically doesn't apply to the $t_1$- and $t_2$-distribution, because these don't have existing variances. However your residuals do not look quite as evil as residuals of these tend to look.
If you run regressions, fitting Gaussians to heavy tailed distributions (even $t_1$ or $t_2$) will normally make confidence intervals and p-values larger, so if you have small p-values, these should still be reliable. It is far more dangerous to have outliers/heavy tails in the explanatory variables ("leverage points"), but as far as I understand your experiment, you don't have those.
