1
$\begingroup$

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!

Notations:

  • $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$

Likelihood:

$\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$

$\endgroup$

closed as off-topic by Tim, gung, Xi'an, John, Nick Cox Jul 20 '15 at 22:13

  • This question does not appear to be about statistics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\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
1
$\begingroup$

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".

$\endgroup$

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