I'm entering this discussion late so I may be discussing issues that have already been addressed.
As I understand it, Federico, you began by describing a model for data in which two factors have a fixed set of levels and one factor has an extensible, or "random", set of levels and you wanted to fit a model that you described as y ~ effect1 * effect2 * effect3 The problem is that this specification is not complete. An interaction of factors with fixed levels and a factor with random levels can mean, in the lmer specification, lmer(y ~ effect1 * effect2 + (1| effect3) + (1|effect1:effect2:effect3), ...) or lmer(y ~ effect1 * effect2 + (effect1*effect2 | effect3), ...) or other variations. When you specify a random effect or an random interaction term you must, either explicitly or implicitly, specify the form of the variance-covariance matrix associated with those random effects. The "advantage" that other software may provide for you is that it chooses the model for you but that, of course, means that you only have the one choice. If you can describe how many variance components you think should be estimated in your model and what they would represent then I think it will be easier to describe how to fit the model. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.