Some symmetric and heavy-tailed distributions, such as Student's-t and Slash, have been suggested for robust inference in linear mixed models. These robust models are characterized by the degrees of freedom of these distributions, and include the normal distribution when the degrees of freedom approach infinity. This simulation study investigated joint estimation of degrees of freedom for the residual and all other genetic and non-genetic parameters in the Slash distributed residual datasets. Bivariate data with heavy-tailed distributed residuals were generated using Slash distributions with 4 or 12 degrees of freedom. Models with bivariate Student's-t, Slash and normal residuals were fitted to each dataset using a hierarchical Bayesian approach. Predictive log-likelihood values strongly favoured the bivariate Student's-t and Slash models over the normal models for simulated heavy-tailed datasets. Posterior mean estimates of degrees of freedom parameters seemed to be accurate and unbiased. Estimates of sire and herd variances were similar, if not identical, across fitted models. Posterior mean and 95% posterior probability interval estimates of error variances in simulated datasets were found to be similar. Reliable estimates of degrees of freedom were obtained in all simulated datasets. The predictive log-likelihood was able to distinguish the correct model among the models fitted to heavy-tailed datasets.