[R] Validating Minitab's "Expanded Gage R&R Study" using R and lme4

[Matt Jacob]

I'm trying to validate the results of an "Expanded Gage R&R Study" in
Minitab using R and lme4, but I can't get the numbers to match up in
certain situations. I can't tell whether my model is wrong, my data is
bad, or something else is going on.

For instance, here's some data for which the results don't match:

https://i.stack.imgur.com/5PCgm.png

After running the gage study, these are the results according to
Minitab:

                                Study Var  %Study Var  %Tolerance
Source             StdDev (SD)   (6 * SD)       (%SV)  (SV/Toler)
Total Gage R&R         1.76277    10.5766      100.00       14.36
  Repeatability        0.00000     0.0000        0.00        0.00
  Reproducibility      1.76277    10.5766      100.00       14.36
    B                  0.00000     0.0000        0.00        0.00
    A*B                1.76277    10.5766      100.00       14.36
Part-To-Part           0.00000     0.0000        0.00        0.00
  A                    0.00000     0.0000        0.00        0.00
Total Variation        1.76277    10.5766      100.00       14.36

But when I mimic Minitab's results by parsing the output from lmer() and
doing the arithmetic in Excel, this is what I see:

https://i.stack.imgur.com/EGg9F.png

The raw output from lmer() was:

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ 1 + (1 | A) + (1 | B) + (1 | A:B)
   Data: d

REML criterion at convergence: -100.1

Scaled residuals: 
       Min         1Q     Median         3Q        Max 
-1.308e-07 -1.308e-07 -1.308e-07 -6.541e-08  1.308e-07 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 A:B      (Intercept) 1.333e+00 1.154e+00
 B        (Intercept) 7.066e-04 2.658e-02
 A        (Intercept) 2.260e-03 4.754e-02
 Residual             2.655e-14 1.629e-07
Number of obs: 8, groups:  A:B, 4; B, 2; A, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)    52.17       0.57   91.53
convergence code: 0
Model failed to converge with max|grad| = 0.422755 (tol = 0.002,
component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

And the R code that produced that output is:

library(lme4)
A <- factor(c(1, 1, 2, 2, 2, 1, 2, 1))
B <- factor(c(1, 2, 1, 2, 1, 2, 2, 1))
y <- c(51.356124843620798, 51.356124843620798, 54.8816618912481,
51.356124843620798, 54.8816618912481, 51.356124843620798,
51.356124843620798, 51.356124843620798)
d <- data.frame(y, A, B)
fm <- lmer(y ~ 1 + (1|A) + (1|B) + (1|A:B), d)
summary(fm)

For a different measurement with a different response, it's a completely
different situation! Given the following data:

https://i.stack.imgur.com/cH0bO.png

The resulting table from Minitab is:

                                Study Var  %Study Var  %Tolerance
Source             StdDev (SD)   (6 * SD)       (%SV)  (SV/Toler)
Total Gage R&R        0.193649    1.16190       55.90        1.00
  Repeatability       0.093541    0.56125       27.00        0.48
  Reproducibility     0.169558    1.01735       48.95        0.88
    B                 0.132288    0.79373       38.19        0.68
    A*B               0.106066    0.63640       30.62        0.55
Part-To-Part          0.287228    1.72337       82.92        1.49
  A                   0.287228    1.72337       82.92        1.49
Total Variation       0.346410    2.07846      100.00        1.79

And after plugging my R results into Excel, I get exactly the same
thing:

https://i.stack.imgur.com/jUEAP.png

Which was produced by this R code:

library(lme4)
A <- factor(c(1, 1, 2, 2, 2, 1, 2, 1))
B <- factor(c(1, 2, 1, 2, 1, 2, 2, 1))
y <- c(-49.4, -49.8, -50.1, -50.1, -50.0, -49.9, -50.2, -49.6)
d <- data.frame(y, A, B)
fm <- lmer(y ~ 1 + (1|A) + (1|B) + (1|A:B), d)
summary(fm)

That generated the following lmer() summary:

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ 1 + (1 | A) + (1 | B) + (1 | A:B)
   Data: d

REML criterion at convergence: -3.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.7705 -0.6853 -0.1039  0.4379  1.4151 

Random effects:
 Groups   Name        Variance Std.Dev.
 A:B      (Intercept) 0.01125  0.10607 
 B        (Intercept) 0.01750  0.13229 
 A        (Intercept) 0.08250  0.28723 
 Residual             0.00875  0.09354 
Number of obs: 8, groups:  A:B, 4; B, 2; A, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept) -49.8875     0.2322  -214.9

Is the difference attributable to the warnings produced by lmer() about
the model failing to converge and being nearly unidentifiable? What
could Minitab be doing differently when the measurement data contains
only two distinct values?

Matt

This question is cross-posted to
http://stats.stackexchange.com/questions/262170/how-can-i-validate-minitabs-expanded-gage-rr-study-using-open-source-tools