[R] Validating Minitab's "Expanded Gage R&R Study" using R and lme4
Matt Jacob
matt at jacobmail.org
Thu Feb 16 02:14:39 CET 2017
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
More information about the R-help
mailing list