[R-sig-ME] Random effects with custom design matrices
Arnaud Le Rouzic
Arnaud.Le-Rouzic at legs.cnrs-gif.fr
Fri Feb 18 17:15:04 CET 2011
Hi the list,
Sorry if this is a FAQ, I could not find anything corresponding in the
archive of the list.
Along with some colleagues, we were wondering if there was a simple way
to customize the design matrix of random effects in lmer. The particular
context was that of a "diallel" experiment in quantitative genetics,
which consists in crossing N populations each other to measure the
hybrid phenotype. The model is
y_ijk = mu + g_i + g_j + d_ij + e_ijk
and the "funny" thing is that g_i and g_j are 1/2 the effect of the two
parental populations, i.e. two occurrences of the same random effect.
d_ij is the interaction effect, and it is non-0 only if i and j are
different.
It is quite straightforward to write the design matrices. Let's consider
a simple data set, like:
father mother phen
1 1 1 1.56
2 1 1 -0.09
3 1 2 2.18
4 1 2 0.73
5 1 3 0.87
6 1 3 -0.06
If the model is written Y = X1 G + X2 D, matrix X1 is
1 1 0 0
2 1 0 0
3 0.5 0.5 0
4 0.5 0.5 0
5 0.5 0 0.5
6 0.5 0 0.5
etc.,
and matrix X2 is
1 0 0 0
2 0 0 0
3 1 0 0
4 1 0 0
5 0 1 0
6 0 1 0
etc
However, lmer crashes when providing the matrices as such:
lmer(phen~1+(1|X1) + (1|X2), data=dd)
There seems to be a workaround in both cases, for X1 each line can be
duplicated in the data and the parent coded as a regular factor, and for
X2 the 0 lines can be removed both from the data vector and from the
design matrix, so that, again, it can be coded as a vector of factors.
However, I cannot find a way to fit both effects at the same time, which
is quite embarrassing. An alternative would be to fit an animal model
with MCMCglmm, but it looks like killing a fly with a bazooka (and on
the top of it, I can't find any automatic way to include dominance in
the animal model).
I am quite sure I am missing something, but I don't know if this is a
pure R issue or if there is something in the underlying statistics that
makes it impossible to handle "non standard" design matrices.
Thanks for your help!
Arnaud.
More information about the R-sig-mixed-models
mailing list