[R-sig-ME] Model representation for mixed models

[Paolo Innocenti]

Hi,

I am fitting a mixed model with lme4, specifically:

model <- lmer(Y ~ female.type*male.type+(1|female.pop)+(1|male.pop))

with female.pop nested in female.type and
male.pop nested in male.type:

 > model.frame(model)
               Y female.type male.type female.pop male.pop
M1M1-1 3.851804        Mono      Mono      Mono1    Mono1
M1M1-2 3.628088        Mono      Mono      Mono1    Mono1
M1P1-1 3.619024        Mono      Prom      Mono1    Prom1
M1P1-2 3.676578        Mono      Prom      Mono1    Prom1
M2M2-1 3.561432        Mono      Mono      Mono2    Mono2
M2M2-2 3.982403        Mono      Mono      Mono2    Mono2
M2P2-1 3.728523        Mono      Prom      Mono2    Prom2
M2P2-2 3.754911        Mono      Prom      Mono2    Prom2
M3M3-1 3.868304        Mono      Mono      Mono3    Mono3
M3M3-2 3.875949        Mono      Mono      Mono3    Mono3
M3P3-1 3.873787        Mono      Prom      Mono3    Prom3
M3P3-2 3.573503        Mono      Prom      Mono3    Prom3
M4M4-1 3.691399        Mono      Mono      Mono4    Mono4
M4M4-2 3.774606        Mono      Mono      Mono4    Mono4
M4P4-1 3.901790        Mono      Prom      Mono4    Prom4
M4P4-2 3.882356        Mono      Prom      Mono4    Prom4
P1M1-1 3.710269        Prom      Mono      Prom1    Mono1
P1M1-2 3.690842        Prom      Mono      Prom1    Mono1
P1P1-1 3.735968        Prom      Prom      Prom1    Prom1
P1P1-2 3.575529        Prom      Prom      Prom1    Prom1
P2M2-1 3.521016        Prom      Mono      Prom2    Mono2
P2M2-2 3.888614        Prom      Mono      Prom2    Mono2
P2P2-1 3.566240        Prom      Prom      Prom2    Prom2
P2P2-2 3.799283        Prom      Prom      Prom2    Prom2
P3M3-1 3.947283        Prom      Mono      Prom3    Mono3
P3M3-2 3.508592        Prom      Mono      Prom3    Mono3
P3P3-1 3.670293        Prom      Prom      Prom3    Prom3
P3P3-2 3.817314        Prom      Prom      Prom3    Prom3
P4M4-1 3.624735        Prom      Mono      Prom4    Mono4
P4M4-2 3.892305        Prom      Mono      Prom4    Mono4
P4P4-1 3.864690        Prom      Prom      Prom4    Prom4
P4P4-2 3.604446        Prom      Prom      Prom4    Prom4

and "contr.sum" as contrasts.

I need to write this model in a rigorous statistical way, but my 
knowledge of stats/math/algebra is really poor, so I find it quite 
difficult. I can guess it is something like:

y_(ijkhw) = u + F_i + M_j + I_(ij) + a_(ik) + b_(jh) + e_(ijkhw)

i = 1,2
j = 1,2
k = 1,...,4
h = 1,...,4
w = 1,...,32

a (female.pop) ~ N(0, (s^2)_1)
b (male.pop) ~ N(0, (s^2)_2)
e (error) ~ N(0, s^2)

I read in ch. 2.1 of the Pinheiro and Bates about the "Laird and Ware 
(1982)" formulation based on two matrices, one for the fixed and the 
random effects, something like:

y_i = X_i Beta + Z_i b_i + e_i

but I cannot really figure out how to apply it to my dataset (which 
should have partially crossed random factors).

My questions:

1 - which formulation is most appropriate/clear?

2 - where can I learn how to write the models properly (and which 
information I have to give to explain my model fully)?

3 - the "model" object contains all the information necessary to produce 
the formal structure of the model in a "pretty" (LaTeX maybe?) format. 
Is there a function that takes the model object and spit out the formal 
representation of the model? I think it would be extremely helpful (and 
pedagogic) for non-statistician users of R.

Best,
paolo

 > model.matrix(model)
       [,1] [,2] [,3] [,4]
  [1,]    1   -1   -1    1
  [2,]    1   -1   -1    1
  [3,]    1   -1    1   -1
  [4,]    1   -1    1   -1
  [5,]    1   -1   -1    1
  [6,]    1   -1   -1    1
  [7,]    1   -1    1   -1
  [8,]    1   -1    1   -1
  [9,]    1   -1   -1    1
[10,]    1   -1   -1    1
[11,]    1   -1    1   -1
[12,]    1   -1    1   -1
[13,]    1   -1   -1    1
[14,]    1   -1   -1    1
[15,]    1   -1    1   -1
[16,]    1   -1    1   -1
[17,]    1    1   -1   -1
[18,]    1    1   -1   -1
[19,]    1    1    1    1
[20,]    1    1    1    1
[21,]    1    1   -1   -1
[22,]    1    1   -1   -1
[23,]    1    1    1    1
[24,]    1    1    1    1
[25,]    1    1   -1   -1
[26,]    1    1   -1   -1
[27,]    1    1    1    1
[28,]    1    1    1    1
[29,]    1    1   -1   -1
[30,]    1    1   -1   -1
[31,]    1    1    1    1
[32,]    1    1    1    1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$female.type
[1] "contr.sum"

attr(,"contrasts")$male.type
[1] "contr.sum"




-- 
Paolo Innocenti
Department of Animal Ecology, EBC
Uppsala University
Norbyvägen 18D
75236 Uppsala, Sweden