6.4 - Fully Nested Random Effects Model: Quality Control Example

6.4 - Fully Nested Random Effects Model: Quality Control Example

The temperature of a process in a manufacturing industry is critical to quality control. The researchers want to characterize the sources of this variability. They choose 4 plants, 4 operators within each plant, look at 4 shifts for each operator, and then measure temperature for each of the three batches used in production.

Collected data was read into SAS and proc mixed procedure was used to obtain the ANOVA model.

data fullnest;
input Temp Plant Operator Shift Batch;
datalines;
477   1   1   1   1
472   1   1   1   2
481   1   1   1   3
478   1   1   2   1
475   1   1   2   2
474   1   1   2   3
472   1   1   3   1
475   1   1   3   2
468   1   1   3   3
482   1   1   4   1
477   1   1   4   2
474   1   1   4   3
471   1   2   1   1
474   1   2   1   2
470   1   2   1   3
479   1   2   2   1
482   1   2   2   2
477   1   2   2   3
470   1   2   3   1
477   1   2   3   2
483   1   2   3   3
480   1   2   4   1
473   1   2   4   2
478   1   2   4   3
475   1   3   1   1
472   1   3   1   2
470   1   3   1   3
460   1   3   2   1
469   1   3   2   2
472   1   3   2   3
477   1   3   3   1
483   1   3   3   2
475   1   3   3   3
476   1   3   4   1
480   1   3   4   2
471   1   3   4   3
465   1   4   1   1
464   1   4   1   2
471   1   4   1   3
477   1   4   2   1
475   1   4   2   2
471   1   4   2   3
481   1   4   3   1
477   1   4   3   2
475   1   4   3   3
470   1   4   4   1
475   1   4   4   2
474   1   4   4   3
484   2   1   1   1
477   2   1   1   2
481   2   1   1   3
477   2   1   2   1
482   2   1   2   2
481   2   1   2   3
479   2   1   3   1
477   2   1   3   2
482   2   1   3   3
477   2   1   4   1
470   2   1   4   2
479   2   1   4   3
472   2   2   1   1
475   2   2   1   2
475   2   2   1   3
472   2   2   2   1
475   2   2   2   2
470   2   2   2   3
472   2   2   3   1
477   2   2   3   2
475   2   2   3   3
482   2   2   4   1
477   2   2   4   2
483   2   2   4   3
485   2   3   1   1
481   2   3   1   2
477   2   3   1   3
482   2   3   2   1
483   2   3   2   2
485   2   3   2   3
477   2   3   3   1
476   2   3   3   2
481   2   3   3   3
479   2   3   4   1
476   2   3   4   2
485   2   3   4   3
477   2   4   1   1
475   2   4   1   2
476   2   4   1   3
476   2   4   2   1
471   2   4   2   2
472   2   4   2   3
475   2   4   3   1
475   2   4   3   2
472   2   4   3   3
481   2   4   4   1
470   2   4   4   2
472   2   4   4   3
475   3   1   1   1
470   3   1   1   2
469   3   1   1   3
477   3   1   2   1
471   3   1   2   2
474   3   1   2   3
469   3   1   3   1
473   3   1   3   2
468   3   1   3   3
477   3   1   4   1
475   3   1   4   2
473   3   1   4   3
470   3   2   1   1
466   3   2   1   2
468   3   2   1   3
471   3   2   2   1
473   3   2   2   2
476   3   2   2   3
478   3   2   3   1
480   3   2   3   2
474   3   2   3   3
477   3   2   4   1
471   3   2   4   2
469   3   2   4   3
466   3   3   1   1
465   3   3   1   2
471   3   3   1   3
473   3   3   2   1
475   3   3   2   2
478   3   3   2   3
471   3   3   3   1
469   3   3   3   2
471   3   3   3   3
475   3   3   4   1
477   3   3   4   2
472   3   3   4   3
469   3   4   1   1
471   3   4   1   2
468   3   4   1   3
473   3   4   2   1
475   3   4   2   2
473   3   4   2   3
477   3   4   3   1
470   3   4   3   2
469   3   4   3   3
463   3   4   4   1
471   3   4   4   2
469   3   4   4   3
484   4   1   1   1
477   4   1   1   2
480   4   1   1   3
476   4   1   2   1
475   4   1   2   2
474   4   1   2   3
475   4   1   3   1
470   4   1   3   2
469   4   1   3   3
481   4   1   4   1
476   4   1   4   2
472   4   1   4   3
469   4   2   1   1
475   4   2   1   2
479   4   2   1   3
482   4   2   2   1
483   4   2   2   2
479   4   2   2   3
477   4   2   3   1
479   4   2   3   2
475   4   2   3   3
472   4   2   4   1
476   4   2   4   2
479   4   2   4   3
470   4   3   1   1
481   4   3   1   2
481   4   3   1   3
475   4   3   2   1
470   4   3   2   2
475   4   3   2   3
469   4   3   3   1
477   4   3   3   2
482   4   3   3   3
485   4   3   4   1
479   4   3   4   2
474   4   3   4   3
469   4   4   1   1
473   4   4   1   2
475   4   4   1   3
477   4   4   2   1
473   4   4   2   2
471   4   4   2   3
470   4   4   3   1
468   4   4   3   2
474   4   4   3   3
483   4   4   4   1
477   4   4   4   2
476   4   4   4   3
;
proc mixed data=fullnest covtest method=type3;
class Plant Operator Shift Batch;
model temp=;
random plant operator(plant) shift(plant operator) ;
run;

In the SAS code, notice that there are no terms on the right-hand side of the model statement. This is because SAS uses the model statement to specify fixed effects only. The random statement is used to specify the random effects. The proc mixed procedure will perform the fully nested random effects model as specified above, and produces the following output:

 
Type 3 Analysis of Variance
Source DF Sum of Squares Mean Square Expected Mean Square Error Term Error DF F Value Pr > F
Plant 3 731.515625 243.838542 Var(Residual) + 3 Var(Shift(Plant*Operato)) + 12 Var(Operator(Plant)) + 48 Var(Plant) MS(Operator(Plant)) 12 5.85 0.0106
Operator(Plant) 12 499.812500 41.651042 Var(Residual) + 3 Var(Shift(Plant*Operato)) + 12 Var(Operator(Plant)) MS(Shift(Plant*Operato)) 48 1.30 0.2483
Shift(Plant*Operato) 48 1534.916667 31.977431 Var(Residual) + 3 Var(Shift(Plant*Operato)) MS(Residual) 128 2.58 <.0001
Residual 128 1588.000000 12.406250 Var(Residual) . . . .
 
Covariance Parameter Estimates
Cov Parm Estimate Standard
Error
Z Value Pr Z
Plant 4.2122 4.1629 1.01 0.3116
Operator(Plant) 0.8061 1.5178 0.53 0.5953
Shift(Plant*Operato) 6.5237 2.2364 2.92 0.0035
Residual 12.4063 1.5508 8.00 <.0001

 

The largest (and significant) variance components are: (1) the shift within a plant \(\times\) operator combination and (2) the batch-to-batch variation within the shift (the residual).

Note that the Covariance Parameter Estimates here are in fact the variance components. SAS does not express the variance components as percentages in this procedure, but by summing the variance components for all sources to serve as the denominator, each source can be expressed as a percentage.

Because this type of model is so commonly employed, SAS also offers two other procedures to obtain the variance components results: proc varcomp (which stands for variance components) and proc nested.

The equivalent code for these procedures is as follows:

The proc varcomp:

proc varcomp data=fullnest;
class  Plant Operator Shift Batch;
model temp= plant operator(plant) shift(plant operator);
run;

Note that the model statement for proc varcomp differs from the mixed procedure in that proc varcomp assumes that the factors listed in the model statement are random effects.

Partial Output:

 
MIVQUE(0) Estimates
Variance Component Temp
Var(Plant) 4.21224
Var(Operator(Plant)) 0.80613
Var(Shift(Plant*Operato)) 6.52373
Var(Error) 12.40625

Note that, even in this procedure we will have to use the sum for a total and calculate the percentages ourselves. On the other hand, the proc nested procedure will provide the full output including the percentages:

proc nested data=fullnest;
class plant operator shift;
var temp;
run;

Partial Output:

 
Nested Random Effects Analysis of Variance for Variable Temp
Variance Source DF Sum of Squares F Value Pr > F Error Term Mean Square Variance Component Percent of Total
Total 191 4354.244792       22.797093 23.948351 100.0000
Plant 3 731.515625 5.85 0.0106 Operator 243.838542 4.212240 17.5889
Operator 12 499.812500 1.30 0.2483 Shift 41.651042 0.806134 3.3661
Shift 48 1534.916667 2.58 <.0001 Error 31.977431 6.523727 27.2408
Error 128 1588.000000       12.406250 12.406250 51.8042

Calculation of the Variance Components

From the SAS output, we get the EMS coefficients. We can use those to compute the estimated variance components.

Source MS EMS Variance Components % Variation
Plant 243.84 \(\sigma_{\epsilon}^{2}+3\sigma_{\gamma}^{2}+12\sigma_{\beta}^{2}+48\sigma_{\alpha}^{2}\) 4.21 17.58
Operator(Plant) 41.65 \(\sigma_{\epsilon}^{2}+3\sigma_{\gamma}^{2}+12\sigma_{\beta}^{2}\) 0.806 3.37
Shift(Plant × Operator) 31.98 \(\sigma_{\epsilon}^{2}+3\sigma_{\gamma}^{2}\) 6.52 27.24
Residual 12.41 \(\sigma_{\epsilon}^{2}\) 12.41 51.80
    Total 23.95  

One can show that MS is an unbiased estimator for EMS (using the properties of Method of Moments estimates). With that, we can algebraically solve for each variance component. Start at the bottom of the table and work up the hierarchy.

First of all, the estimated variance component for the Residuals is given:

\(\mathbf{12.41} = \hat{\sigma}_{\text{error}}^{2}\) = \(\hat{\sigma}_{\epsilon}^{2}\)

Then we can use this information and subtract it from the Shift(Plant \(\times\) Operator) MS to get:

\(\begin{align}31.98 &= \hat{\sigma}_{\epsilon}^{2}+3\hat{\sigma}_{\gamma \text{ or Shift(Plant*Operator)}}^{2} \\[10pt] \hat{\sigma}_{\gamma}^{2}&=\frac{31.98-12.41}{3}=\mathbf{6.52}\end{align}\)

Similarly, we use what we know for Error and Shift(Plant \(\times\) Operator) and subtract it from the Operator(Plant) MS to get:

\(\begin{align}41.65&=\hat{\sigma}_{\epsilon}^{2}+3\hat{\sigma}_{\gamma}^{2} +12\hat{\sigma}_{\beta \text{ or Operator(Plant)}}^{2}\\[10pt] &=31.98+12\hat{\sigma}_{\beta}^{2}\\[10pt] \hat{\sigma}_{\beta}^{2}&=\frac{41.65-31.98}{12}\\[10pt] &=\mathbf{0.806}\end{align}\)

And finally, we use what we know for Error and Shift(Plant \(\times\) Operator) and Operator(Plant), subtract it from the Plant MS to get:

\(\begin{align}243.84&=\hat{\sigma}_{\epsilon}^{2}+3\hat{\sigma}_{\gamma}^{2} +12\hat{\sigma}_{\beta}^{2}+48\hat{\sigma}_{\alpha \text{ or Plant}}^{2}\\[10pt] &=41.65+48\hat{\sigma}_{\alpha}^{2}\\[10pt] \hat{\sigma}_{\alpha}^{2}&=\frac{243.84-41.65}{48}\\[10pt] &=\mathbf{4.21}\end{align}\)

Our total = 12.41 + 6.52 + 0.806 + 4.21 = 23.95

Then, dividing each variance component by the total gives the % values shown in the output from SAS proc nested.

Minitab has a separate program just for this type of analysis for our example (Quality Data ), under:

Stat > ANOVA > Fully Nested ANOVA

and you specify the model in the boxes provided:

Minitab Fully Nested ANOVA window with C4 shift chosen

The output you get is very comprehensive and includes the variance components expressed as percentages.

Nested ANOVA: Temp versus Plant, Operator, Shift

Analysis of Variance for Temp

Source DF SS MS F P
Plant 3 731.5156 243.8385 5.854 0.011
Operator 12 499.8125 41.6510 1.303 0.248
Shift 48 1534.9167 31.9774 2.578 0.000
Error 128 1588.0000 12.4062    
Total 191 4354.2448      

Variance Components

Source Var Comp. # of Total StDev
Plant 4.212 17.59 2.052
Operator 0.806 3.37 0.898
Shift 6.524 27.24 2.554
Error 12.406 51.80 3.522
Total 23.948   4.894

Expected Mean Squares

1 Plant 1.00(4) + 3.00(3) + 12.00(2) + 48.00(1)
2 Operator 1.00(4) + 3.00(3) + 12.00(2)
3 Shift 1.00(4) + 3.00(3)
4 Error 1.00(4)

First, load and attach the quality data in R.

setwd("~/path-to-folder/")
quality_data<-read.table("quality_data.txt",header=T,sep ="\t")
attach(quality_data)

We then fit the fully nested random effects model. There are no fixed effects, so we display only the LRT table for the random effects.

library(lme4,lmerTest)
lm_rand = lmer(Temp ~ (1 | Plant) + (1 | Plant:Operator) + (1 | Plant:(Operator:Shift)))
ranova(lm_rand)
ANOVA-like table for random-effects: Single term deletions

Model:
Temp ~ (1 | Plant) + (1 | Plant:Operator) + (1 | Plant:(Operator:Shift))
                             npar  logLik    AIC     LRT Df Pr(>Chisq)    
                                5 -548.59 1107.2                          
(1 | Plant)                     4 -551.03 1110.1  4.8755  1    0.02724 *  
(1 | Plant:Operator)            4 -548.77 1105.5  0.3530  1    0.55240    
(1 | Plant:(Operator:Shift))    4 -557.36 1122.7 17.5317  1  2.826e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

NTo examine the variance components further we print the estimates.


vc = VarCorr(lm_rand, parm="theta_", oldNames = FALSE)
print(vc,comp=c("Variance"))
 Groups                 Name        Variance
 Plant:(Operator:Shift) (Intercept)  6.52371
 Plant:Operator         (Intercept)  0.80614
 Plant                  (Intercept)  4.21227
 Residual                           12.40626

We can compute the variance percentages manually.

terms = data.frame(vc)$grp
vars = data.frame(vc)$vcov
tot_var = sum(vars)
percent_var = (vars/tot_var)*100
data.frame(cbind(terms,vars,percent_var))
                   terms              vars      percent_var
1 Plant:(Operator:Shift)  6.52370980810034 27.2407289699471
2         Plant:Operator 0.806136528703026 3.36614400964098
3                  Plant  4.21226520994903 17.5889452947895
4               Residual  12.4062556768176 51.8041817256224

Note R does not produce the EMS coefficients, so if those are of interest another software may be useful.


Legend
[1]Link
Has Tooltip/Popover
 Toggleable Visibility