# [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.