[R-sig-ME] Random effects with custom design matrices

[Arnaud Le Rouzic]

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.