2 factor factorial model with random factors

Question

Using lme4, how does one define a full 2-factor factorial model with both factors (and their interaction) being random?

Specifically I am trying to recreate the results from Montgomery's Design of Experiments book (7th edition) example 13.2. In this example there are 2 random factors and I want to include the interaction in the model as Montgomery tests for significance in the full model first. I've tried several things but cannot recreate the results in R. I would think something like what's given below would work, but it does not.

lmer(y ~ (1|Parts) + (1|Operators) + (1|Parts:Operators))

The results I'm seeing are provided below given the model specification above. The results show that the variance component estimate for Parts:Operators is 0 suggesting that the interaction model isn’t specified correctly, or that R is not handling the model correctly.

Random effects:
 Groups          Name        Variance Std.Dev.
 Parts:Operators (Intercept)  0.00000 0.0000  
 Parts           (Intercept) 10.25127 3.2018  
 Operators       (Intercept)  0.01063 0.1031  
 Residual                     0.88316 0.9398

In contrast, the results I expect to see using REML are shown below (these are the variance component estimates shown by Montgomery that I can reproduce in JMP but not R).

Random effects:
 Groups          Name        Variance Std.Dev.
 Parts:Operators (Intercept) -0.13991 0.1219  
 Parts           (Intercept) 10.27983 3.3738  
 Operators       (Intercept)  0.01491 0.0330  
 Residual                     0.99167 0.1811

y ~ A + B + A*B where A and B are both random factors... a 2-factor full factorial model where both factors are random — Jake, Dec 18 '14 at 13:25
And what is the problem with recreating results? remember that: (a) Montgomery used different model (ANOVA) and (d) different software, so there would be nothing strange if the results were not exactly the same. — Tim, Dec 18 '14 at 14:12
So different software, yes... I can use JMP or Minitab and get the same results, I'm justing trying to figure out how I can also do this in R, just a learning exercise. As for the different method, I should get the same results using ANOVA or REML because it is a completely balanced design. As a matter of fact, he shows the same results (same variance components) using both methods. — Jake, Dec 18 '14 at 14:16
This seems like a reasonable specification. Can you provide more details, i.e. results vs. expected results? — Ben Bolker, Dec 19 '14 at 03:36
@BenBolker, I added my results vs. expected results to the original post. — Jake, Jan 06 '15 at 16:50

score 3 · Accepted Answer · edited Apr 13 '17 at 12:44

The reason for this difference is that lme4 does not allow negative variances. The Parts:Operators variance component given by Montgomery is negative. It is not hard to come up with situations where the method-of-moments estimator of a particular variance component is negative (e.g. when the observed variance among group means is less that $\sigma^2_r/n$, where $\sigma^2_r$ is the observed within-group variance and $n$ is the number of observations per group). There have been extensive discussions on the r-sig-mixed-models@r-project.org mailing list (e.g. here) about why one would like to exclude negative variances (e.g., negative variances are mathematically impossible) or include them (e.g., they are sometimes a convenient shortcut). The 'modern' way to handle these situations is to fit a compound-symmetry model where the within-group correlations are allowed to be negative; this is implemented in lme, but so far only in the experimental branch of lme4. A similar but not identical situation is addressed (showing an lme example) in this CV question (I would be curious what SAS PROC MIXED does in this situation ...)

A minor note: it looks like the "Standard deviation" column that you quote from Montgomery is giving the standard errors of the variance components. In lme4 the "Standard deviation" column gives the standard deviation of the random effect, i.e. just the square root of the variances. (There is another can of worms here, about whether the standard error of variance components is, or is not, a good indicator of the uncertainty ... a reasonable starting point is here.)

2 factor factorial model with random factors

1 Answers1