I would like to fit a basic "animal model" (a kind of linear mixed model, see below) using the latest version of lme4 on CRAN (v1.7). To check the results, I am fitting it on simulated data, and comparing the results with the rrBLUP package. However, the results of lme4 are really strange and I don't succeed in understanding why (even after reading the vignette and this answer). My R code is available as a GitHub Gist here. Any help would be much appreciated!


  • $N$: number of animals, known

  • $\mu$: global mean, unknown

  • $P$: number of fixed effects, known

  • $X$: $N \times P$ design matrix of fixed effects, known

  • $\boldsymbol{b}$: $P$-dimensional vector of fixed effects, unknown

  • $W = [ \boldsymbol{1} \; X]$: $N \times (P+1)$ matrix, known

  • $a = [\mu \; \boldsymbol{b}']'$: $(P+1)$ vector, unknown

  • $Q$: number of "genetics" random effects (in this document, $Q=N$), known

  • $Z$: $N \times Q$ design matrix of "genetics" random effects, known

  • $\sigma_u^2$: variance component of the "genetics" random effects, unknown

  • $A$: $Q \times Q$ matrix of additive relationships (obtained from pedigree or from molecular markers), known

  • $G = \sigma_u^2 A$: $Q \times Q$ covariance matrix of the "genetic" random effects, unknown

  • $\boldsymbol{u}$: $Q$-dimensional vector (often called "breeding values"), unknown

  • $\sigma^2$: variance component of the errors, unknown

  • $R = \sigma^2 I_N$: $N \times N$ covariance matrix of the errors, unknown

  • $\boldsymbol{e}$: $N$-dimensional vector of errors, unknown

  • $\mathcal{N}_N$: multivariate Normal distribution of dimension $N$


$\boldsymbol{y} = W \boldsymbol{a} + Z \boldsymbol{u} + \boldsymbol{e}$ where $\boldsymbol{u} \sim \mathcal{N}_N(\boldsymbol{0}, G)$ and $\boldsymbol{e} \sim \mathcal{N}_N(\boldsymbol{0}, R)$

We also assume zero covariance between $\boldsymbol{u}$ and $\boldsymbol{e}$:

$\begin{pmatrix} \boldsymbol{y} \\ \boldsymbol{u} \\ \boldsymbol{e} \end{pmatrix} = \mathcal{N} \begin{pmatrix} W \boldsymbol{a}, & V & ZG & R \\ \boldsymbol{0}, & GZ' & G & \boldsymbol{0} \\ \boldsymbol{0}, & R & \boldsymbol{0} & R \end{pmatrix} $

where $V = Var(\boldsymbol{y}) = ZGZ' + R$

Primary goal: estimate $\sigma_u^2$ and $\sigma^2$

  • $\begingroup$ hi, sorry I didn't answer your question via e-mail. I will take a crack at this if I get a chance. $\endgroup$ – Ben Bolker Jan 22 '15 at 1:46
  • $\begingroup$ the only thing I can think of after a brief look at the code is that the translation from the internal theta parameters to the estimated variance-covariance matrix via VarCorr is not working properly/as expected (in general, one can't rely on the standard accessor methods working normally if the internal structure of the model is very different from what's expected). Is there a scaling factor missing? I would take a look at lme4:::VarCorr.merMod and see if it's doing what you think it should be doing ... $\endgroup$ – Ben Bolker Jan 22 '15 at 4:19
  • $\begingroup$ @BenBolker no need to apologize! I perfectly understand you're busy. I'll look at lme4:::VarCorr.merMod in the mean time. $\endgroup$ – tflutre Jan 22 '15 at 7:41
  • 4
    $\begingroup$ I am voting to close this thread since OP declared that the problem was due the bug in his/her code, it was not a statistical problem. $\endgroup$ – Tim Jul 20 '15 at 12:05
  • $\begingroup$ @Tim I completely understand and am sorry for it, but does "closing the thread" means "removing it"? Because other people may be interested in keeping this example somewhere. For instance, someone up-voted the question. $\endgroup$ – tflutre Jul 20 '15 at 12:19

My bad, there was a bug in my script (available as a gist). On line 97, instead of relmat <- list(animal=G), I should have written relmat <- list(animal=A). This is now corrected in the gist (v2). So the good news is that lme4 (v1.7) indeed works quite well to fit a basic "animal model".

| cite | improve this answer | |

Not the answer you're looking for? Browse other questions tagged or ask your own question.