2 factor factorial model with random factors 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  

 A: 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.)
