[R] lmer model specification

[Ben Bolker]

Roy <jroyrobertson <at> comcast.net> writes:

> 
> I am running version 2.15.2 64 bit version on 64 bit Windows 7. I have a 
> data set with the following structure:
> Fixed Effect: locationFact
> Random Effects: datefact, timefact nested in datefact, interactions of 
> datefact and timefact with locationFact

  I think that corresponds to

  locationFact + (locationFact|datefact/timefact)

(if you want to allow for correlations between the main 
effect and interaction term of locationFact among datefact levels and
among timefact levels) or more likely (assuming the
interactions are independent)

  locationFact + (1|(1+locationFact):(datefact/timefact))

which should be equivalent to 

  locationFact + (1|datefact) + (1|datefact:timefact) +
        (1|locationFact:datefact)+  (1|locationFact:datefact:timefact)

> I fit the model with the latest version of lme4.
> 
> The formula is:   Thick2 ~ locationFact + (1 | datefact) + (1 | 
> datefact/timefact) +      (1 | locationFact:datefact) + (1 | 
> datefact/locationFact:timefact)

  this expands to

 locationFact + 
 (1|datefact) + (1|datefact + datefact:timefact) +
 (1|locationFact:datefact) + (1 + datefact + datefact:locationFact:timefact)

  You can see that you do indeed have datefact included three
times here ...

> 
> Other elements of  output object are:
> 
> Linear mixed model fit by REML
> 
> Random effects:
>   Groups                                            Name Variance Std.Dev.
>   locationFact:timefact:datefact  (Intercept)  34.614   5.8833
>   timefact:datefact                        (Intercept)  96.795 9.8385
>   locationFact:datefact                 (Intercept)  56.375   7.5083
>   datefact                                       (Intercept) 20.341   
> 4.5101
>   datefact                                       (Intercept) 20.338   
> 4.5098
>   datefact                                       (Intercept) 20.340   
> 4.5100
>   Residual 447.252  21.1483
> I tested this model using rmel with another software package and found 
> that the datefact variance is the sum of the 3 datefact variance 
> estimates above. Is there a way to specify the model so I do not get the 
> 3 datefact estimates? The same applies to the datefact BLUPs. I have to 
> add them to get the actual BLUP.

One way to explore how this works is to use model.matrix, e.g.

> df <- factor(1:2)
> lf <- factor(1:2)
> tf <- factor(1:2)
> colnames(model.matrix(~df/lf:tf))
 [1] "(Intercept)" "df2"         "df1:lf1:tf1" "df2:lf1:tf1" "df1:lf2:tf1"
 [6] "df2:lf2:tf1" "df1:lf1:tf2" "df2:lf1:tf2" "df1:lf2:tf2" "df2:lf2:tf2"