3.4 PER MILLION
MSA in Two Dimensions
Capturing measurement variation
when two variables are
by Joseph D. Conklin
Introducing measurement system analysis (MSA) with single variable examples is good for communicating basic concepts of the subject.1 After practitioners are comfortable with the single variable case, graduating to the multiple variable one is the next logical step. Understanding measurement variation in two or more dimensions simultaneously is a regularly occurring need.
For a flavor of the multiple dimensional case, consider the example of a metal fabricator producing rectangular components with tolerances for length and width to yield a desired surface area. The variation in surface area is a function of the joint variation of length and width, and the measurement variation of surface area is a function of the measurement variations of length and width, and their interplay.
Background and assumptions
Assume the components have equal tolerance for length and width of 30 mm +/- 2 mm. As is common in the single variable case, we will assess measurement variation by a combination of control charts and by a comparison of the observed variation with the component tolerances.
The control charts indicate the stability of the measuring process. The magnitude of the measuring variation determines the size of change the measuring system can realistically detect.
Control charts accommodating multivariate data are required. Instead of looking at the variance of a single variable, we work with the variance-covariance matrix of length and width variables. One traditional technique that lends itself to multivariate measurement capability uses a quantity called Hotelling’s T2 statistic.2
Harold Hotelling was a statistician who developed the technique. T2 refers to the fact that in some cases, this statistic can be expressed as the square of a t-distribution.3 When assessing multivariate measurement capability, using Hotelling’s T2 statistic entails several assumptions regarding the data:
- The data follow a multivariate normal distribution. While strong, this assumption over time has proven a reasonable one for some applications and is a good starting point when introducing measurement variation of multiple variables. We will adopt this assumption for the example covered here. Tests exist for seeing whether data follows a multivariate normal.4 These are beyond the scope of this example and will not be covered.
- The true variance-covariance matrix is unknown and must be estimated from the available data.
- Independent, replicated measurements are available. In other words, we have to be able to measure a single object multiple times, one measurement immediately following the other, such that the results of the current measurement do not influence the results of the next and the object itself is unchanged after measurement.
For control chart purposes, a set of replicated measurements acts as a rational subgroup.5 The number of replicates or subgroups in an actual study would be determined by sample size formulas.6 For purposes of illustration, we will work with 10 subgroups of 16 replicates each.
Summary statistics and Hotelling’s T2 control chart
The necessary summary statistics for length are listed in Online Tables 1-2. They are based on the data in Online Table 3. The summary statistics based on the width are found in Online Table 4. Four tools will assess the measuring system. First, a multivariate control chart sees whether the composite variation of length and width changes in a small, random manner.
We also will compare the composite variation to the tolerances so we can compute a measurement capability ratio. For supplemental information, we also will construct separate control charts for length and width. This last step is usually needed to determine what is causing a multivariate chart to go out of control.
Three things are required for computing Hotelling’s T2 statistic: the matrix of variances and covariances, the inverse of this matrix and a quadratic form involving this inverse, as well as a vector of deviations from the sample means of the variables. The variance-covariance matrix is calculated for each subgroup.7
To form the variance-covariance matrix, construct a 2 x 2 matrix with the variances along the principal diagonal and the covariances in the off-diagonal positions. In the more general case of p variables, the variance of the ith variable is placed in the ith element of the principal diagonal. The ijth cell is the covariance of variable i with variable j. The order of the matrix is p x p, and it is symmetric.
Online Table 5 is the variance-covariance matrixes by subgroup for our example data. Cell 1, 1 is the variance of the length in the subgroup. Cells 1, 2 and 2, 1 are the covariance of length and width. Cell 2, 2 is the variance of width in the subgroup. In our example, the days are subgroups.
After we have the variance-covariance matrixes for each subgroup, the next step is to form the average variance-covariance matrix. The ijth cell in the average matrix, as the name implies, is the average of the ijth cells across the subgroup matrixes. For the data in our exercise, the average variance-covariance matrix works out to be:
The inverse of this average variance-covariance matrix, denoted A, is that unique matrix, denoted A-1, that when multiplied by A, becomes the identity matrix with ones along the principal diagonal and zeroes elsewhere. The inverse A-1 for the example data works out to be:
In practice, the matrixes are of higher order, and the calculations are carried out via computer programs supporting this function.8 The last information for computing the T2 statistic requires the vector of deviations from the mean. The vector of deviations is formed by subtracting the mean of each variable from the raw data. So here, we subtract the grand mean of the length measurements from the mean of the length data for each day. We perform a similar calculation for the width data.
The grand mean of all the length data is 30.82 mm, and the grand mean of all the width data is 30.81 mm (see Online Table 6).
The values for Hotelling’s T2 statistic are found by computing the quadratic form (n)(d)(A-1)(d’). The vector of deviations is d, A-1 is the inverse matrix, d’ is the transpose of d, and n is the number of observations in a single subgroup.
The vector of deviations d is a 2 x 1 vector in this example. The inverse matrix is the 2 x 2 inverse calculated earlier. d’ is a 1 x 2 vector, the transpose of d. n, the number of observations in a single subgroup, is 16, the number of replicate measurements per day.9
Hotelling’s T2 value results are computed as follows:
- Premultiply the inverse by d to produce a 1 x 2 vector.
- Multiply the resulting 1 x 2 vector by the 2 x 1 vector d’.
- Multiply the resulting scalar quantity by n for Hotelling’s T2 statistic. The calculations of Hotelling’s T2 statistics for the days forming various subgroups are shown in Online Table 7.
A multivariate control chart based on Hotelling’s T2 statistic plots the T2 values for each subgroup in time order between the lower and upper control limits. Because the values are the results of a quadratic form, the smallest possible value is zero. The formula for the lower control limit is a function of the F distribution:
[p(n-1)]/(n-p) x F[(1-(a/2)), p, n-p]
in which n = the number of values in a single subgroup, p = the number of variables, and F[(1-( a/2)), p, n-p] is the value of an F distribution with p and n-p degrees of freedom with upper tail area equal to (1-( a//2)).10
The formula for the upper control limit is also a function of the F distribution:
[p(n-1)]/(n-p) x F(a/2, p, n-p)
in which F(a/2, p, n-p) is the value of an F distribution with p and n-p degrees of freedom with upper tail area equal to a/2.
For this example, we have n = 16 and p = 2. For a, let’s use the level implied by the traditional three sigma control limits on a Shewhart X̅ chart: 0.0027. The F distribution value needed for the lower control limit is F(1-0.00135, 2, 16-2) or F(0.99865, 2, 14). The value for the upper control limit is F(0.00135, 2, 14). F(0.099865, 2, 14) = 0.001351. F(0.00135, 2, 14) = 10.99.
The lower and upper control limits for the multivariate control chart, therefore, are [2(16-1)]/(16-2) x 0.001351 = 0.0029 and (30/14) x 10.99 = [2(16-1)]/(16-2) x 10.99 = (30/14) x 0.001351 = 23.55. Figure 1 is the multivariate control chart for our data.
The random movement between the control limits of the T2 chart suggests the measurement process is stable and predictable.
Separate control charts for length and width
As an additional check, Figures 2-5 are X––and standard deviation (s) control charts for the length and width data. The data and calculations for the charts are shown in Online Tables 8-9.11
The length and width data move randomly inside the limits of their respective control charts. This supports the conclusion of Hotelling’s T2 chart.
The individual control charts will agree with the T2 chart when the variables are independent of one another. In other words, a change in one variable does not influence how another variable might change.
In real life, variables can be correlated. In general, the larger the number of variables, the greater the chance of some correlation occurring among at least some of them.
When variables are correlated, the signals from the individual charts might disagree with the T2 chart. In that case, the signal from the T2 should govern because it takes into account the correlations between the variables.
The T2 chart is not the only multivariate control chart available, but it also works well in a variety of situations. Information on other multivariate control charts can be found by searching the quality control literature.12
With the control chart analysis complete, we can move to the last step: calculating the total variability and a multivariate measurement capability ratio.
Total measurement variation and capability ratio
The total measurement variation and capability ratio for the multiple variable case is more complicated and less standardized than for the single variable case. Decisions are required about what probability distribution is appropriate and how much of the probability distribution should be reflected in the total figure for measurement variation.
If we say the multivariate result for measurement variation should capture the traditional probability in the common single variable analysis—99.73%—the multivariate normal for our example data leads to a calculation involving the area of an ellipse. The area of the ellipse works out to be 10.7193 mm2.13 See the sidebar, "Capturing 99.73% of the Probability for a Bivariate Normal" for more details.
Capturing 99.73% of the
Probability for a Bivariate Normal
The distribution assumed for the measurement variation is a bivariate normal. Interpret the total measurement variation as the cumulative probability encompassed by an ellipse of constant Mahalanobis distance r2. The cumulative probability of interest is 0.9973.
Using the formula at the bottom of p. 4 of "N-Dimensional Cumulative Function, and Other Useful Facts About Gaussians and Normal Densities,"1 the desired r2 equals (-2 x ln(1 - 0.9973) = (-2 x ln(0.0027)) = (-2 x -5.9145) = 11.8290. By the formula for Mahalanobis distance on p. 2 of "N-Dimensional Cumulative Function, and Other Useful Facts about Gaussians and Normal Densities,"2 compute r2 for our example data is presented as:
==> 3.20435L2 - 0.46962LW- 0.46962LW + 3.75054W2 = 11.8290, where L = length and W = width.
The area of the ellipse spanning the required Malahanobis distance is π x a x b, in which a is the length of the semimajor axis and b is the length of the semiminor axis. These can be found by setting l or w in the earlier expression to 0 and solving.
Setting L to 0 reduces the expression to:
3.75054W2 = 11.8290 ==>
W2 = 3.1539 ==>
W = 1.7759.
Setting W to 0 reduces the expression to:
3.20435 L2 = 11.8290 ==>
L2 = 3.6915 ==>
L = 1.9213.
The length of the semimajor axis is 1.9213 mm, and the length of the semiminor axis is 1.7759 mm.
The area of the ellipse works out to be 3.1416 x 1.9213 mm x 1.7759 mm = 10.7193 mm2, the result for the measurement variation corresponding to the desired probability of 0.9973. —J.D.C.
- Michael Bensimhoun, "N-Dimensional Cumulative Function and Other Useful
Facts About Gaussians and Normal Densities," June 2009, https://upload.wikimedia.org/wikipedia/commons/a/a2/Cumulative_
function_n_dimensional_Gaussians_12.2013.pdf, p. 4.
- Ibid, p. 2.
The area of the ellipse, 10.7193 mm2, is compared to the surface area formed by the tolerances for length and width around their respective targets. The surface area is 16 mm2, as shown in Figure 6.
The measurement capability ratio follows as:
10.7193 mm2/16mm2 = 0.6700 = 67%.
This measurement system is unacceptable for assuring conformance to the customer requirements.14 The best way to set the performance standard for a measuring system is to decide how large of a change matters to the data users. Then at least the results of the measurement variation study will show the ability of the system to detect such a change.
Addressing the measuring system
With the stable, predictable behavior shown on the control charts and with a measurement capability ratio at a level generally regarded as unacceptable, the data support redesign of the measuring system. More in-depth study will help decide whether this should be done by using new equipment, developing new operating procedures, enacting new maintenance and training efforts—or some combination of these.
In cases in which the control charts show a lack of control, there is also evidence of potential problems with the measuring system. At this point, the performance standards of the measurement system should be changed. If this is not an option, resources must be invested to improve it.
Measurement results for data that were obtained during a period when the measurement system was not performing as required should be checked against some other system known to be acceptable. If a check shows the original data did not accurately represent the process at the time of measurement, corrective action should be considered and implemented.
References and Notes
- For illustrations of measurement system studies for single variable examples, see Automotive Industry Action’s (AIAG) Measurement Systems Analysis, fourth edition, AIAG, 2010.
- For more information about variance-covariance matrixes and Hotelling’s T2 distribution, see Kanti V. Mardia, J.T. Kent and J. M. Bibby’s Multivariate Analysis, Academic Press, 1979.
- The t-distribution is a fundamental one in the statistics of single variables. For more information, see George W. Snedecor and William G. Cochran’s Statistical Methods, eighth edition, Iowa State University Press, 1989.
- For a technical discussion of tests for multivariate normality, see D.R. Cox and N.J.H. Small’s "Testing Multivariate Normality," Biometrika, Vol. 65, No. 2, 1978, pp. 263-272.
- For details of rational subgroups and other aspects of control charts, see Eugene L. Grant and Richard S. Leavenworth’s Statistical Quality Control, seventh edition, McGraw-Hill, 1996.
- For some of the theory behind the sample size formulas for measurement system analysis (MSA), see the chapter on random effects models in Snedecor and Cochran’s Statistical Methods, reference 3. For examples of sample size formulas in more general contexts, see William G. Cochran’s Sampling Techniques, third edition, John Wiley & Sons Inc., 1977.
- Variance-covariances matrixes, their inverses and quadratic forms are part of the context of linear models, another body of theory contributing to MSA. For more detailed coverage of these and other important and related concepts, see Shayle R. Searle’s Linear Models, John Wiley & Sons Inc., 1971.
- Two statistical programs useful for this purpose are SAS, a product of SAS Institute and R, a product of the R Project for Statistical Computing, available for free download at www.r-project.org.
- For the context behind these calculations, see Searle, Linear Models, reference 7.
- For more information about the F distribution and its uses, see Snedecor and Cochran’s Statistical Methods, reference 3. Tail areas of the F distribution are available in tables published in this and similar books. They also may be obtained from the statistics programs mentioned in reference 8.
- For more information on the calculations shown in Online Table 4, including tables of factors for limits on x̅ and s control charts, see Grant and Leavenworth’s Statistical Quality Control, reference 5.
- For an example, see Cynthia A. Lowry and Douglas C. Montgomery’s "A Review of Multivariate Control Charts," IIE Transactions, Vol. 27, Issue 6, 1995, pp. 800-810.
- The calculations in the online sidebar "Capturing 99.73% of the Probability for a Bivariate Normal" are based on the discussion in "N-Dimensional Cumulative Function, and Other Useful Facts About Gaussians and Normal Densities" at http://tinyurl.com/gaussians.
- The percentage calculated here is a smaller than the better criterion. For more in depth discussion of acceptable performance levels for measurement systems, see AIAG’s Measurement Systems Analysis, reference 1. What is acceptable depends on the importance of the quality characteristic and the maximum tolerable risk of wrong conclusions about the measurement system.
Joseph D. Conklin is a mathematical statistician in Washington, D.C. He earned a master’s degree in statistics from Virginia Tech in Blacksburg and is a senior member of ASQ. Conklin is also an ASQ-certified quality manager, engineer, auditor, reliability engineer, and Six Sigma Black Belt.