A quick tour of qcc

Luca Scrucca

10 Jul 2017

Introduction

qcc is a contributed R package for statistical quality control charts which provides:

This document gives a quick tour of qcc (version 2.7) functionalities. It was written in R Markdown, using the knitr package for production.

Further details are provided in the following paper:

Scrucca, L. (2004) qcc: an R package for quality control charting and statistical process control. R News 4/1, 11-17.

For a nice blog post discussing the qcc package, in particular how to implement the Western Eletric Rules (WER), see http://blog.yhathq.com/posts/quality-control-in-r.html.

library(qcc)

Shewhart charts

x-bar chart

data(pistonrings)
diameter = with(pistonrings, qcc.groups(diameter, sample))
head(diameter)
##     [,1]   [,2]   [,3]   [,4]   [,5]
## 1 74.030 74.002 74.019 73.992 74.008
## 2 73.995 73.992 74.001 74.011 74.004
## 3 73.988 74.024 74.021 74.005 74.002
## 4 74.002 73.996 73.993 74.015 74.009
## 5 73.992 74.007 74.015 73.989 74.014
## 6 74.009 73.994 73.997 73.985 73.993
q1 = qcc(diameter[1:25,], type="xbar", newdata=diameter[26:40,])

plot(q1, chart.all=FALSE)

plot(q1, add.stats=FALSE)

q1 = qcc(diameter[1:25,], type="xbar", newdata=diameter[26:40,],
         confidence.level=0.99)

Add warning limits at 2 std. deviations:

q1 = qcc(diameter[1:25,], type="xbar", newdata=diameter[26:40,], plot=FALSE)
(warn.limits = limits.xbar(q1$center, q1$std.dev, q1$sizes, 2))
##       LCL      UCL
##  73.99242 74.00993
plot(q1, restore.par = FALSE)
abline(h = warn.limits, lty = 3, col = "chocolate")

R chart

q2 = qcc(diameter[1:25,], type="R")

summary(q2)
## 
## Call:
## qcc(data = diameter[1:25, ], type = "R")
## 
## R chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00800 0.01800 0.02100 0.02276 0.02600 0.03900 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  0.02276
## Standard deviation:  0.009785039 
## 
## Control limits:
##  LCL        UCL
##    0 0.04812533
q3 = qcc(diameter[1:25,], type="R", newdata=diameter[26:40,])

summary(q3)
## 
## Call:
## qcc(data = diameter[1:25, ], type = "R", newdata = diameter[26:40,     ])
## 
## R chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00800 0.01800 0.02100 0.02276 0.02600 0.03900 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  0.02276
## Standard deviation:  0.009785039 
## 
## Summary of group statistics in diameter[26:40, ]:
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## 0.01400000 0.01900000 0.02500000 0.02453333 0.02750000 0.04400000 
## 
## Group sample size:  5
## Number of groups:  15 
## 
## Control limits:
##  LCL        UCL
##    0 0.04812533

S chart

q4 = qcc(diameter[1:25,], type="S")

summary(q4)
## 
## Call:
## qcc(data = diameter[1:25, ], type = "S")
## 
## S chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##        Min.     1st Qu.      Median        Mean     3rd Qu.        Max. 
## 0.002863564 0.007314369 0.008467585 0.009240037 0.011928956 0.016177144 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  0.009240037
## Standard deviation:  0.009829977 
## 
## Control limits:
##  LCL        UCL
##    0 0.01930242
q5 = qcc(diameter[1:25,], type="S", newdata=diameter[26:40,])

summary(q5)
## 
## Call:
## qcc(data = diameter[1:25, ], type = "S", newdata = diameter[26:40,     ])
## 
## S chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##        Min.     1st Qu.      Median        Mean     3rd Qu.        Max. 
## 0.002863564 0.007314369 0.008467585 0.009240037 0.011928956 0.016177144 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  0.009240037
## Standard deviation:  0.009829977 
## 
## Summary of group statistics in diameter[26:40, ]:
##        Min.     1st Qu.      Median        Mean     3rd Qu.        Max. 
## 0.005310367 0.007367603 0.010329569 0.009761757 0.011232319 0.016546903 
## 
## Group sample size:  5
## Number of groups:  15 
## 
## Control limits:
##  LCL        UCL
##    0 0.01930242

Variable control limits

out = c(9, 10, 30, 35, 45, 64, 65, 74, 75, 85, 99, 100)
diameter2 = with(pistonrings, qcc.groups(diameter[-out], sample[-out]))
summary(qcc(diameter2[1:25,], type="xbar"))

## 
## Call:
## qcc(data = diameter2[1:25, ], type = "xbar")
## 
## xbar chart for diameter2[1:25, ] 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 73.99020 73.99680 73.99980 74.00076 74.00425 74.01020 
## 
## Summary of group sample sizes:                
##   sizes  3 4  5
##   counts 4 4 17
## 
## Number of groups:  25
## Center of group statistics:  74.00075
## Standard deviation:  0.009856701 
## 
## Control limits:
##          LCL      UCL
##     73.98753 74.01398
##     73.98368 74.01782
## ...                  
##     73.98753 74.01398
summary(qcc(diameter2[1:25,], type="R"))

## 
## Call:
## qcc(data = diameter2[1:25, ], type = "R")
## 
## R chart for diameter2[1:25, ] 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00800 0.01600 0.02100 0.02168 0.02600 0.03900 
## 
## Summary of group sample sizes:                
##   sizes  3 4  5
##   counts 4 4 17
## 
## Number of groups:  25
## Center of group statistics:  0.02230088
## Standard deviation:  0.009856701 
## 
## Control limits:
##     LCL        UCL
##       0 0.04785198
##       0 0.04857007
## ...               
##       0 0.04785198

p and np charts

data(orangejuice)
q1 = with(orangejuice, 
          qcc(D[trial], sizes=size[trial], type="p"))

summary(q1)
## 
## Call:
## qcc(data = D[trial], type = "p", sizes = size[trial])
## 
## p chart for D[trial] 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0800000 0.1600000 0.2100000 0.2313333 0.2950000 0.4800000 
## 
## Group sample size:  50
## Number of groups:  30
## Center of group statistics:  0.2313333
## Standard deviation:  0.421685 
## 
## Control limits:
##            LCL       UCL
##     0.05242755 0.4102391
##     0.05242755 0.4102391
## ...                     
##     0.05242755 0.4102391
q2 = with(orangejuice, 
          qcc(D[trial], sizes=size[trial], type="np"))

summary(q2)
## 
## Call:
## qcc(data = D[trial], type = "np", sizes = size[trial])
## 
## np chart for D[trial] 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  4.00000  8.00000 10.50000 11.56667 14.75000 24.00000 
## 
## Group sample size:  50
## Number of groups:  30
## Center of group statistics:  11.56667
## Standard deviation:  2.981763 
## 
## Control limits:
##       LCL      UCL
##  2.621377 20.51196

Remove out-of-control points (see help(orangejuice) for the reasons):

inc = setdiff(which(orangejuice$trial), c(15,23))
q2 = with(orangejuice, 
          qcc(D[inc], sizes=size[inc], type="p",
              newdata=D[!trial], newsizes=size[!trial]))

data(orangejuice2)
q1 = with(orangejuice2,
          qcc(D[trial], sizes=size[trial], type="p", 
              newdata=D[!trial], newsizes=size[!trial]))

summary(q1)
## 
## Call:
## qcc(data = D[trial], type = "p", sizes = size[trial], newdata = D[!trial],     newsizes = size[!trial])
## 
## p chart for D[trial] 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0400000 0.0800000 0.1100000 0.1108333 0.1200000 0.2400000 
## 
## Group sample size:  50
## Number of groups:  24
## Center of group statistics:  0.1108333
## Standard deviation:  0.3139256 
## 
## Summary of group statistics in D[!trial]:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.020   0.080   0.100   0.109   0.140   0.220 
## 
## Group sample size:  50
## Number of groups:  40 
## 
## Control limits:
##     LCL       UCL
##       0 0.2440207
##       0 0.2440207
## ...              
##       0 0.2440207

c and u charts

data(circuit)
q1 = with(circuit, 
          qcc(x[trial], sizes=size[trial], type="c"))

summary(q1)
## 
## Call:
## qcc(data = x[trial], type = "c", sizes = size[trial])
## 
## c chart for x[trial] 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  5.00000 16.00000 19.00000 19.84615 24.00000 39.00000 
## 
## Group sample size:  100
## Number of groups:  26
## Center of group statistics:  19.84615
## Standard deviation:  4.454902 
## 
## Control limits:
##       LCL      UCL
##  6.481447 33.21086

Remove out-of-control points (see help(circuit) for the reasons)

inc = setdiff(which(circuit$trial), c(6,20))
q2 = with(circuit, 
          qcc(x[inc], sizes=size[inc], type="c", labels=inc, 
              newdata=x[!trial], newsizes=size[!trial], newlabels=which(!trial)))

summary(q2)
## 
## Call:
## qcc(data = x[inc], type = "c", sizes = size[inc], labels = inc,     newdata = x[!trial], newsizes = size[!trial], newlabels = which(!trial))
## 
## c chart for x[inc] 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 10.00000 16.00000 19.00000 19.66667 24.00000 31.00000 
## 
## Group sample size:  100
## Number of groups:  24
## Center of group statistics:  19.66667
## Standard deviation:  4.434712 
## 
## Summary of group statistics in x[!trial]:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.00   15.75   18.50   18.30   21.00   28.00 
## 
## Group sample size:  100
## Number of groups:  20 
## 
## Control limits:
##       LCL     UCL
##  6.362532 32.9708
q3 = with(circuit, 
          qcc(x[inc], sizes=size[inc], type="u", labels=inc, 
              newdata=x[!trial], newsizes=size[!trial], newlabels=which(!trial)))

summary(q3)
## 
## Call:
## qcc(data = x[inc], type = "u", sizes = size[inc], labels = inc,     newdata = x[!trial], newsizes = size[!trial], newlabels = which(!trial))
## 
## u chart for x[inc] 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.1000000 0.1600000 0.1900000 0.1966667 0.2400000 0.3100000 
## 
## Group sample size:  100
## Number of groups:  24
## Center of group statistics:  0.1966667
## Standard deviation:  0.4434712 
## 
## Summary of group statistics in x[!trial]:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0900  0.1575  0.1850  0.1830  0.2100  0.2800 
## 
## Group sample size:  100
## Number of groups:  20 
## 
## Control limits:
##         LCL      UCL
##  0.06362532 0.329708
data(pcmanufact)
q1 = with(pcmanufact, 
          qcc(x, sizes=size, type="u"))

summary(q1)
## 
## Call:
## qcc(data = x, type = "u", sizes = size)
## 
## u chart for x 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    1.40    2.00    1.93    2.25    3.20 
## 
## Group sample size:  5
## Number of groups:  20
## Center of group statistics:  1.93
## Standard deviation:  1.389244 
## 
## Control limits:
##         LCL      UCL
##  0.06613305 3.793867

Continuous one-at-time data

# viscosity data (Montgomery, pag. 242)
x = c(33.75, 33.05, 34, 33.81, 33.46, 34.02, 33.68, 33.27, 33.49, 33.20,
       33.62, 33.00, 33.54, 33.12, 33.84)
q1 = qcc(x, type="xbar.one")

summary(q1)
## 
## Call:
## qcc(data = x, type = "xbar.one")
## 
## xbar.one chart for x 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 33.00000 33.23500 33.54000 33.52333 33.78000 34.02000 
## 
## Group sample size:  1
## Number of groups:  15
## Center of group statistics:  33.52333
## Standard deviation:  0.4261651 
## 
## Control limits:
##       LCL      UCL
##  32.24484 34.80183
q2 = qcc(x, type="xbar.one", std.dev = "SD")

summary(q2)
## 
## Call:
## qcc(data = x, type = "xbar.one", std.dev = "SD")
## 
## xbar.one chart for x 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 33.00000 33.23500 33.54000 33.52333 33.78000 34.02000 
## 
## Group sample size:  1
## Number of groups:  15
## Center of group statistics:  33.52333
## Standard deviation:  0.3415928 
## 
## Control limits:
##       LCL      UCL
##  32.49855 34.54811

Standardized p chart

In this example we show how to extend the package by defining a new control chart, i.e. a standardized p chart (type = "p.std").

Function to compute group statistics and center:

stats.p.std = function(data, sizes)
{
  data = as.vector(data)
  sizes = as.vector(sizes)
  pbar = sum(data)/sum(sizes)
  z = (data/sizes - pbar)/sqrt(pbar*(1-pbar)/sizes)
  list(statistics = z, center = 0)
}

Function to compute within-group standard deviation:

sd.p.std = function(data, sizes, ...) { return(1) }

Function to compute control limits based on normal approximation:

limits.p.std = function(center, std.dev, sizes, conf)
{
  if(conf >= 1) 
    { lcl = -conf
      ucl = +conf 
  }
  else
    { if(conf > 0 & conf < 1)
        { nsigmas = qnorm(1 - (1 - conf)/2)
          lcl = -nsigmas
          ucl = +nsigmas }
      else stop("invalid 'conf' argument.") 
  }
  limits = matrix(c(lcl, ucl), ncol = 2)
  rownames(limits) = rep("", length = nrow(limits))
  colnames(limits) = c("LCL", "UCL")
  return(limits)
}

Example with simulated data:

# set unequal sample sizes
n = c(rep(50,5), rep(100,5), rep(25, 5))
# generate randomly the number of successes
x = rbinom(length(n), n, 0.2)
# plot the control chart with variable limits
summary(qcc(x, type="p", size=n))

## 
## Call:
## qcc(data = x, type = "p", sizes = n)
## 
## p chart for x 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.100   0.160   0.200   0.204   0.250   0.320 
## 
## Summary of group sample sizes:                   
##   sizes  25 50 100
##   counts  5  5   5
## 
## Number of groups:  15
## Center of group statistics:  0.208
## Standard deviation:  0.4058768 
## 
## Control limits:
##            LCL       UCL
##     0.03580105 0.3801990
##     0.03580105 0.3801990
## ...                     
##     0.00000000 0.4515261
# plot the standardized control chart
summary(qcc(x, type="p.std", size=n))

## 
## Call:
## qcc(data = x, type = "p.std", sizes = n)
## 
## p.std chart for x 
## 
## Summary of group statistics:
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.8815445 -0.7859403 -0.0985521 -0.0527780  0.7145025  1.3797289 
## 
## Summary of group sample sizes:                   
##   sizes  25 50 100
##   counts  5  5   5
## 
## Number of groups:  15
## Center of group statistics:  0
## Standard deviation:  1 
## 
## Control limits:
##  LCL UCL
##   -3   3

Operating Characteristic Curves

An operating characteristic curve graphically provides information about the probability of not detecting a shift in the process. The function oc.curves() is a generic function which calls the proper function depending on the type of input 'qcc' object.

data(pistonrings)
diameter = with(pistonrings, qcc.groups(diameter, sample))
q = qcc(diameter, type="xbar", nsigmas=3, plot=FALSE)
beta = oc.curves.xbar(q)

print(round(beta, digits=4))
##                sample size
## shift (std.dev)    n=5    n=1   n=10   n=15   n=20
##            0    0.9973 0.9973 0.9973 0.9973 0.9973
##            0.05 0.9971 0.9973 0.9970 0.9968 0.9966
##            0.1  0.9966 0.9972 0.9959 0.9952 0.9944
##            0.15 0.9957 0.9970 0.9940 0.9920 0.9900
##            0.2  0.9944 0.9968 0.9909 0.9869 0.9823
##            0.25 0.9925 0.9964 0.9864 0.9789 0.9701
##            0.3  0.9900 0.9960 0.9798 0.9670 0.9514
##            0.35 0.9866 0.9956 0.9708 0.9500 0.9243
##            0.4  0.9823 0.9950 0.9586 0.9266 0.8871
##            0.45 0.9769 0.9943 0.9426 0.8957 0.8383
##            0.5  0.9701 0.9936 0.9220 0.8562 0.7775
##            0.55 0.9616 0.9927 0.8963 0.8078 0.7055
##            0.6  0.9514 0.9916 0.8649 0.7505 0.6243
##            0.65 0.9390 0.9905 0.8275 0.6853 0.5371
##            0.7  0.9243 0.9892 0.7842 0.6137 0.4481
##            0.75 0.9071 0.9877 0.7351 0.5379 0.3616
##            0.8  0.8871 0.9860 0.6809 0.4608 0.2817
##            0.85 0.8642 0.9842 0.6225 0.3851 0.2115
##            0.9  0.8383 0.9821 0.5612 0.3136 0.1527
##            0.95 0.8094 0.9798 0.4983 0.2485 0.1059
##            1    0.7775 0.9772 0.4355 0.1913 0.0705
##            1.05 0.7428 0.9744 0.3743 0.1431 0.0450
##            1.1  0.7055 0.9713 0.3161 0.1038 0.0275
##            1.15 0.6659 0.9678 0.2622 0.0730 0.0161
##            1.2  0.6243 0.9641 0.2134 0.0497 0.0090
##            1.25 0.5812 0.9599 0.1703 0.0328 0.0048
##            1.3  0.5371 0.9554 0.1333 0.0209 0.0024
##            1.35 0.4925 0.9505 0.1022 0.0129 0.0012
##            1.4  0.4481 0.9452 0.0768 0.0077 0.0006
##            1.45 0.4043 0.9394 0.0564 0.0045 0.0002
##            1.5  0.3616 0.9332 0.0406 0.0025 0.0001
##            1.55 0.3206 0.9265 0.0286 0.0013 0.0000
##            1.6  0.2817 0.9192 0.0197 0.0007 0.0000
##            1.65 0.2453 0.9115 0.0133 0.0003 0.0000
##            1.7  0.2115 0.9032 0.0088 0.0002 0.0000
##            1.75 0.1806 0.8943 0.0056 0.0001 0.0000
##            1.8  0.1527 0.8849 0.0036 0.0000 0.0000
##            1.85 0.1278 0.8749 0.0022 0.0000 0.0000
##            1.9  0.1059 0.8643 0.0013 0.0000 0.0000
##            1.95 0.0869 0.8531 0.0008 0.0000 0.0000
##            2    0.0705 0.8413 0.0004 0.0000 0.0000
##            2.05 0.0566 0.8289 0.0002 0.0000 0.0000
##            2.1  0.0450 0.8159 0.0001 0.0000 0.0000
##            2.15 0.0353 0.8023 0.0001 0.0000 0.0000
##            2.2  0.0275 0.7881 0.0000 0.0000 0.0000
##            2.25 0.0211 0.7734 0.0000 0.0000 0.0000
##            2.3  0.0161 0.7580 0.0000 0.0000 0.0000
##            2.35 0.0121 0.7422 0.0000 0.0000 0.0000
##            2.4  0.0090 0.7257 0.0000 0.0000 0.0000
##            2.45 0.0066 0.7088 0.0000 0.0000 0.0000
##            2.5  0.0048 0.6915 0.0000 0.0000 0.0000
##            2.55 0.0034 0.6736 0.0000 0.0000 0.0000
##            2.6  0.0024 0.6554 0.0000 0.0000 0.0000
##            2.65 0.0017 0.6368 0.0000 0.0000 0.0000
##            2.7  0.0012 0.6179 0.0000 0.0000 0.0000
##            2.75 0.0008 0.5987 0.0000 0.0000 0.0000
##            2.8  0.0006 0.5793 0.0000 0.0000 0.0000
##            2.85 0.0004 0.5596 0.0000 0.0000 0.0000
##            2.9  0.0002 0.5398 0.0000 0.0000 0.0000
##            2.95 0.0002 0.5199 0.0000 0.0000 0.0000
##            3    0.0001 0.5000 0.0000 0.0000 0.0000
##            3.05 0.0001 0.4801 0.0000 0.0000 0.0000
##            3.1  0.0000 0.4602 0.0000 0.0000 0.0000
##            3.15 0.0000 0.4404 0.0000 0.0000 0.0000
##            3.2  0.0000 0.4207 0.0000 0.0000 0.0000
##            3.25 0.0000 0.4013 0.0000 0.0000 0.0000
##            3.3  0.0000 0.3821 0.0000 0.0000 0.0000
##            3.35 0.0000 0.3632 0.0000 0.0000 0.0000
##            3.4  0.0000 0.3446 0.0000 0.0000 0.0000
##            3.45 0.0000 0.3264 0.0000 0.0000 0.0000
##            3.5  0.0000 0.3085 0.0000 0.0000 0.0000
##            3.55 0.0000 0.2912 0.0000 0.0000 0.0000
##            3.6  0.0000 0.2743 0.0000 0.0000 0.0000
##            3.65 0.0000 0.2578 0.0000 0.0000 0.0000
##            3.7  0.0000 0.2420 0.0000 0.0000 0.0000
##            3.75 0.0000 0.2266 0.0000 0.0000 0.0000
##            3.8  0.0000 0.2119 0.0000 0.0000 0.0000
##            3.85 0.0000 0.1977 0.0000 0.0000 0.0000
##            3.9  0.0000 0.1841 0.0000 0.0000 0.0000
##            3.95 0.0000 0.1711 0.0000 0.0000 0.0000
##            4    0.0000 0.1587 0.0000 0.0000 0.0000
##            4.05 0.0000 0.1469 0.0000 0.0000 0.0000
##            4.1  0.0000 0.1357 0.0000 0.0000 0.0000
##            4.15 0.0000 0.1251 0.0000 0.0000 0.0000
##            4.2  0.0000 0.1151 0.0000 0.0000 0.0000
##            4.25 0.0000 0.1056 0.0000 0.0000 0.0000
##            4.3  0.0000 0.0968 0.0000 0.0000 0.0000
##            4.35 0.0000 0.0885 0.0000 0.0000 0.0000
##            4.4  0.0000 0.0808 0.0000 0.0000 0.0000
##            4.45 0.0000 0.0735 0.0000 0.0000 0.0000
##            4.5  0.0000 0.0668 0.0000 0.0000 0.0000
##            4.55 0.0000 0.0606 0.0000 0.0000 0.0000
##            4.6  0.0000 0.0548 0.0000 0.0000 0.0000
##            4.65 0.0000 0.0495 0.0000 0.0000 0.0000
##            4.7  0.0000 0.0446 0.0000 0.0000 0.0000
##            4.75 0.0000 0.0401 0.0000 0.0000 0.0000
##            4.8  0.0000 0.0359 0.0000 0.0000 0.0000
##            4.85 0.0000 0.0322 0.0000 0.0000 0.0000
##            4.9  0.0000 0.0287 0.0000 0.0000 0.0000
##            4.95 0.0000 0.0256 0.0000 0.0000 0.0000
##            5    0.0000 0.0228 0.0000 0.0000 0.0000
data(orangejuice)
q = with(orangejuice, qcc(D[trial], sizes=size[trial], type="p", plot=FALSE))
beta = oc.curves(q)

print(round(beta, digits=4))
##      0   0.01   0.02   0.03   0.04   0.05   0.06   0.07   0.08   0.09 
## 0.0000 0.0894 0.2642 0.4447 0.5995 0.7206 0.8100 0.8735 0.9173 0.9468 
##    0.1   0.11   0.12   0.13   0.14   0.15   0.16   0.17   0.18   0.19 
## 0.9662 0.9788 0.9869 0.9920 0.9951 0.9971 0.9983 0.9990 0.9993 0.9995 
##    0.2   0.21   0.22   0.23   0.24   0.25   0.26   0.27   0.28   0.29 
## 0.9995 0.9993 0.9987 0.9978 0.9962 0.9937 0.9900 0.9845 0.9768 0.9662 
##    0.3   0.31   0.32   0.33   0.34   0.35   0.36   0.37   0.38   0.39 
## 0.9522 0.9343 0.9118 0.8844 0.8518 0.8139 0.7711 0.7236 0.6722 0.6176 
##    0.4   0.41   0.42   0.43   0.44   0.45   0.46   0.47   0.48   0.49 
## 0.5610 0.5035 0.4461 0.3901 0.3365 0.2862 0.2398 0.1980 0.1609 0.1287 
##    0.5   0.51   0.52   0.53   0.54   0.55   0.56   0.57   0.58   0.59 
## 0.1013 0.0784 0.0596 0.0446 0.0327 0.0235 0.0166 0.0115 0.0078 0.0052 
##    0.6   0.61   0.62   0.63   0.64   0.65   0.66   0.67   0.68   0.69 
## 0.0034 0.0021 0.0013 0.0008 0.0005 0.0003 0.0002 0.0001 0.0000 0.0000 
##    0.7   0.71   0.72   0.73   0.74   0.75   0.76   0.77   0.78   0.79 
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
##    0.8   0.81   0.82   0.83   0.84   0.85   0.86   0.87   0.88   0.89 
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
##    0.9   0.91   0.92   0.93   0.94   0.95   0.96   0.97   0.98   0.99 
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
##      1 
## 0.0000
data(circuit)
q = with(circuit, qcc(x[trial], sizes=size[trial], type="c", plot=FALSE))
beta = oc.curves(q)

print(round(beta, digits=4))
##      0      1      2      3      4      5      6      7      8      9 
## 0.0000 0.0006 0.0166 0.0839 0.2149 0.3840 0.5543 0.6993 0.8088 0.8843 
##     10     11     12     13     14     15     16     17     18     19 
## 0.9329 0.9625 0.9797 0.9893 0.9945 0.9972 0.9986 0.9991 0.9992 0.9986 
##     20     21     22     23     24     25     26     27     28     29 
## 0.9972 0.9945 0.9895 0.9813 0.9686 0.9502 0.9249 0.8918 0.8505 0.8011 
##     30     31     32     33     34     35     36     37     38     39 
## 0.7444 0.6818 0.6150 0.5461 0.4772 0.4102 0.3470 0.2888 0.2365 0.1907 
##     40     41     42     43     44     45     46     47     48     49 
## 0.1514 0.1184 0.0912 0.0693 0.0519 0.0383 0.0280 0.0201 0.0143 0.0101 
##     50     51     52     53     54     55     56     57     58     59 
## 0.0070 0.0048 0.0033 0.0022 0.0015 0.0010 0.0006 0.0004 0.0003 0.0002 
##     60     61     62     63     64     65 
## 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000

Cusum chart

data(pistonrings)
diameter = with(pistonrings, qcc.groups(diameter, sample))

q1 = cusum(diameter[1:25,], decision.interval = 4, se.shift = 1)

summary(q1)
## 
## Call:
## cusum(data = diameter[1:25, ], decision.interval = 4, se.shift = 1)
## 
## cusum chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   73.99   74.00   74.00   74.00   74.00   74.01 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  74.00118
## Standard deviation:  0.009785039 
## 
## Decision interval (std.err.): 4 
## Shift detection  (std. err.): 1
q2 = cusum(diameter[1:25,], newdata=diameter[26:40,])

summary(q2)
## 
## Call:
## cusum(data = diameter[1:25, ], newdata = diameter[26:40, ])
## 
## cusum chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   73.99   74.00   74.00   74.00   74.00   74.01 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  74.00118
## Standard deviation:  0.009785039 
## 
## Summary of group statistics in diameter[26:40, ]:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 73.99220 74.00290 74.00720 74.00765 74.01270 74.02340 
## 
## Group sample size:  5
## Number of groups:  15 
## 
## Decision interval (std.err.): 5 
## Shift detection  (std. err.): 1
plot(q2, chart.all=FALSE)

EWMA

data(pistonrings)
diameter = with(pistonrings, qcc.groups(diameter, sample))

q1 = ewma(diameter[1:25,], lambda=0.2, nsigmas=3)

summary(q1)
## 
## Call:
## ewma(data = diameter[1:25, ], lambda = 0.2, nsigmas = 3)
## 
## ewma chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 73.99020 73.99820 74.00080 74.00118 74.00420 74.01020 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  74.00118
## Standard deviation:  0.009785039 
## 
## Smoothing parameter: 0.2 
## Control limits:
##          LCL      UCL
## 1   73.99855 74.00380
## 2   73.99781 74.00454
## ...                  
## 25  73.99680 74.00555
q2 =  ewma(diameter[1:25,], lambda=0.2, nsigmas=2.7, 
            newdata=diameter[26:40,]) 

summary(q2)
## 
## Call:
## ewma(data = diameter[1:25, ], lambda = 0.2, nsigmas = 2.7, newdata = diameter[26:40,     ])
## 
## ewma chart for diameter[1:25, ] 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 73.99020 73.99820 74.00080 74.00118 74.00420 74.01020 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  74.00118
## Standard deviation:  0.009785039 
## 
## Summary of group statistics in diameter[26:40, ]:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 73.99220 74.00290 74.00720 74.00765 74.01270 74.02340 
## 
## Group sample size:  5
## Number of groups:  15 
## 
## Smoothing parameter: 0.2 
## Control limits:
##          LCL      UCL
## 1   73.99881 74.00354
## 2   73.99815 74.00420
## ...                  
## 40  73.99724 74.00511
x = c(33.75, 33.05, 34, 33.81, 33.46, 34.02, 33.68, 33.27, 
       33.49, 33.20, 33.62, 33.00, 33.54, 33.12, 33.84)
q3 =  ewma(x, lambda=0.2, nsigmas=2.7)

summary(q3)
## 
## Call:
## ewma(data = x, lambda = 0.2, nsigmas = 2.7)
## 
## ewma chart for x 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 33.00000 33.23500 33.54000 33.52333 33.78000 34.02000 
## 
## Group sample size:  1
## Number of groups:  15
## Center of group statistics:  33.52333
## Standard deviation:  0.4261651 
## 
## Smoothing parameter: 0.2 
## Control limits:
##            LCL      UCL
## [1,]  33.29320 33.75346
## [2,]  33.22862 33.81804
## ...                    
## [15,] 33.14002 33.90664

Process capability analysis

data(pistonrings)
diameter = with(pistonrings, qcc.groups(diameter, sample))

q1 = qcc(diameter[1:25,], type="xbar", nsigmas=3, plot=FALSE)
process.capability(q1, spec.limits=c(73.95,74.05))

## 
## Process Capability Analysis
## 
## Call:
## process.capability(object = q1, spec.limits = c(73.95, 74.05))
## 
## Number of obs = 125          Target = 74
##        Center = 74              LSL = 73.95
##        StdDev = 0.009785        USL = 74.05
## 
## Capability indices:
## 
##       Value   2.5%  97.5%
## Cp    1.703  1.491  1.915
## Cp_l  1.743  1.555  1.932
## Cp_u  1.663  1.483  1.844
## Cp_k  1.663  1.448  1.878
## Cpm   1.691  1.480  1.902
## 
## Exp<LSL 0%    Obs<LSL 0%
## Exp>USL 0%    Obs>USL 0%
process.capability(q1, spec.limits=c(73.95,74.05), target=74.02)

## 
## Process Capability Analysis
## 
## Call:
## process.capability(object = q1, spec.limits = c(73.95, 74.05),     target = 74.02)
## 
## Number of obs = 125          Target = 74.02
##        Center = 74              LSL = 73.95
##        StdDev = 0.009785        USL = 74.05
## 
## Capability indices:
## 
##        Value    2.5%   97.5%
## Cp    1.7033  1.4914  1.9148
## Cp_l  1.7433  1.5548  1.9319
## Cp_u  1.6632  1.4827  1.8437
## Cp_k  1.6632  1.4481  1.8783
## Cpm   0.7856  0.6556  0.9154
## 
## Exp<LSL 0%    Obs<LSL 0%
## Exp>USL 0%    Obs>USL 0%
process.capability(q1, spec.limits=c(73.99,74.01))

## 
## Process Capability Analysis
## 
## Call:
## process.capability(object = q1, spec.limits = c(73.99, 74.01))
## 
## Number of obs = 125          Target = 74
##        Center = 74              LSL = 73.99
##        StdDev = 0.009785        USL = 74.01
## 
## Capability indices:
## 
##        Value    2.5%   97.5%
## Cp    0.3407  0.2983  0.3830
## Cp_l  0.3807  0.3176  0.4439
## Cp_u  0.3006  0.2424  0.3588
## Cp_k  0.3006  0.2312  0.3700
## Cpm   0.3382  0.2960  0.3804
## 
## Exp<LSL 13%   Obs<LSL 12%
## Exp>USL 18%   Obs>USL 16%
process.capability(q1, spec.limits = c(73.99, 74.1))

## 
## Process Capability Analysis
## 
## Call:
## process.capability(object = q1, spec.limits = c(73.99, 74.1))
## 
## Number of obs = 125          Target = 74.05
##        Center = 74              LSL = 73.99
##        StdDev = 0.009785        USL = 74.1
## 
## Capability indices:
## 
##        Value    2.5%   97.5%
## Cp    1.8736  1.6406  2.1063
## Cp_l  0.3807  0.3176  0.4439
## Cp_u  3.3665  3.0115  3.7215
## Cp_k  0.3807  0.3055  0.4559
## Cpm   0.4083  0.3377  0.4788
## 
## Exp<LSL 13%   Obs<LSL 12%
## Exp>USL 0%    Obs>USL 0%

Multivariate Quality Control Charts

Multivariate subgrouped data from Ryan (2000, Table 9.2) with \(p = 2\) variables, \(m = 20\) samples, and \(n = 4\) sample size for each sample:

X1 = matrix(c(72, 56, 55, 44, 97, 83, 47, 88, 57, 26, 46, 
49, 71, 71, 67, 55, 49, 72, 61, 35, 84, 87, 73, 80, 26, 89, 66, 
50, 47, 39, 27, 62, 63, 58, 69, 63, 51, 80, 74, 38, 79, 33, 22, 
54, 48, 91, 53, 84, 41, 52, 63, 78, 82, 69, 70, 72, 55, 61, 62, 
41, 49, 42, 60, 74, 58, 62, 58, 69, 46, 48, 34, 87, 55, 70, 94, 
49, 76, 59, 57, 46), ncol = 4)
X2 = matrix(c(23, 14, 13, 9, 36, 30, 12, 31, 14, 7, 10, 
11, 22, 21, 18, 15, 13, 22, 19, 10, 30, 31, 22, 28, 10, 35, 18, 
11, 10, 11, 8, 20, 16, 19, 19, 16, 14, 28, 20, 11, 28, 8, 6, 
15, 14, 36, 14, 30, 8, 35, 19, 27, 31, 17, 18, 20, 16, 18, 16, 
13, 10, 9, 16, 25, 15, 18, 16, 19, 10, 30, 9, 31, 15, 20, 35, 
12, 26, 17, 14, 16), ncol = 4)
X = list(X1 = X1, X2 = X2) # a list of matrices, one for each variable

q = mqcc(X, type = "T2")

summary(q)
## 
## Call:
## mqcc(data = X, type = "T2")
## 
## T2 chart for X 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.12429  1.32496  2.50272  6.47013  5.34912 63.76042 
## 
## Number of variables:  2
## Number of groups:  20
## Group sample size:  4
## 
## Center: 
##      X1      X2 
## 60.3750 18.4875 
## 
## Covariance matrix:
##          X1        X2
## X1 222.0333 103.11667
## X2 103.1167  56.57917
## |S|:  1929.414 
## 
## Control limits:
##  LCL      UCL
##    0 11.03976
ellipseChart(q)

ellipseChart(q, show.id = TRUE)

q = mqcc(X, type = "T2", pred.limits = TRUE)

Ryan (2000) discussed Xbar-charts for single variables computed adjusting the confidence level of the \(T^2\) chart:

q1 = qcc(X1, type = "xbar", confidence.level = q$confidence.level^(1/2))

summary(q1)
## 
## Call:
## qcc(data = X1, type = "xbar", confidence.level = q$confidence.level^(1/2))
## 
## xbar chart for X1 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  40.000  54.000  61.375  60.375  68.250  81.250 
## 
## Group sample size:  4
## Number of groups:  20
## Center of group statistics:  60.375
## Standard deviation:  14.93443 
## 
## Control limits:
##       LCL      UCL
##  37.97352 82.77648
q2 = qcc(X2, type = "xbar", confidence.level = q$confidence.level^(1/2))

summary(q2)
## 
## Call:
## qcc(data = X2, type = "xbar", confidence.level = q$confidence.level^(1/2))
## 
## xbar chart for X2 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 10.5000 15.3750 19.0000 18.4875 21.5000 29.7500 
## 
## Group sample size:  4
## Number of groups:  20
## Center of group statistics:  18.4875
## Standard deviation:  7.139388 
## 
## Control limits:
##     LCL     UCL
##  7.7785 29.1965

Generate new “in control” data:

Xnew = list(X1 = matrix(NA, 10, 4), X2 =  matrix(NA, 10, 4))
for(i in 1:4)
   { x = MASS::mvrnorm(10, mu = q$center, Sigma = q$cov)
     Xnew$X1[,i] = x[,1]
     Xnew$X2[,i] = x[,2] 
   }
qq = mqcc(X, type = "T2", newdata = Xnew, pred.limits = TRUE)

summary(qq)
## 
## Call:
## mqcc(data = X, type = "T2", pred.limits = TRUE, newdata = Xnew)
## 
## T2 chart for X 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.12429  1.32496  2.50272  6.47013  5.34912 63.76042 
## 
## Number of variables:  2
## Number of groups:  20
## Group sample size:  4
## 
## Center: 
##      X1      X2 
## 60.3750 18.4875 
## 
## Covariance matrix:
##          X1        X2
## X1 222.0333 103.11667
## X2 103.1167  56.57917
## |S|:  1929.414 
## 
## Summary of group statistics in Xnew:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.356345 1.148565 1.602482 2.347734 3.318625 5.342685 
## 
## Number of groups:  10
## Group sample size:  4
## 
## Control limits:
##  LCL      UCL
##    0 11.03976
## 
## Prediction limits:
##  LPL      UPL
##    0 12.20184
ellipseChart(qq, show.id = TRUE)

Generate new “out of control” data:

Xnew = list(X1 = matrix(NA, 10, 4), X2 =  matrix(NA, 10, 4))
for(i in 1:4)
   { x = MASS::mvrnorm(10, mu = 1.2*q$center, Sigma = q$cov)
     Xnew$X1[,i] = x[,1]
     Xnew$X2[,i] = x[,2] 
   }
qq = mqcc(X, type = "T2", newdata = Xnew, pred.limits = TRUE)

summary(qq)
## 
## Call:
## mqcc(data = X, type = "T2", pred.limits = TRUE, newdata = Xnew)
## 
## T2 chart for X 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.12429  1.32496  2.50272  6.47013  5.34912 63.76042 
## 
## Number of variables:  2
## Number of groups:  20
## Group sample size:  4
## 
## Center: 
##      X1      X2 
## 60.3750 18.4875 
## 
## Covariance matrix:
##          X1        X2
## X1 222.0333 103.11667
## X2 103.1167  56.57917
## |S|:  1929.414 
## 
## Summary of group statistics in Xnew:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  3.045001  4.935299  5.498651  6.758231  6.989239 16.962860 
## 
## Number of groups:  10
## Group sample size:  4
## 
## Control limits:
##  LCL      UCL
##    0 11.03976
## 
## Prediction limits:
##  LPL      UPL
##    0 12.20184
ellipseChart(qq, show.id = TRUE)

Individual observations data:

data(boiler)

q1 = mqcc(boiler, type = "T2.single", confidence.level = 0.999)

summary(q1)
## 
## Call:
## mqcc(data = boiler, type = "T2.single", confidence.level = 0.999)
## 
## T2.single chart for boiler 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  1.316342  5.305689  7.074224  7.680000  9.775744 17.575293 
## 
## Number of variables:  8
## Number of groups:  25
## Group sample size:  1
## 
## Center: 
##     t1     t2     t3     t4     t5     t6     t7     t8 
## 525.00 513.56 538.92 521.68 503.80 512.44 478.72 477.24 
## 
## Covariance matrix:
##            t1        t2        t3         t4         t5         t6
## t1 54.0000000 0.9583333 20.583333 31.2916667 20.3333333 -2.2916667
## t2  0.9583333 4.8400000  2.963333  2.6866667  0.3250000  3.0766667
## t3 20.5833333 2.9633333 22.993333 10.0566667  4.9833333  2.0366667
## t4 31.2916667 2.6866667 10.056667 22.3100000 13.5583333 -2.2700000
## t5 20.3333333 0.3250000  4.983333 13.5583333 11.4166667 -1.6583333
## t6 -2.2916667 3.0766667  2.036667 -2.2700000 -1.6583333  4.5066667
## t7 20.3750000 0.7050000  6.476667 12.6983333 10.6500000 -0.7466667
## t8  0.2083333 3.4433333  2.770000  0.1633333 -0.2833333  3.7233333
##            t7         t8
## t1 20.3750000  0.2083333
## t2  0.7050000  3.4433333
## t3  6.4766667  2.7700000
## t4 12.6983333  0.1633333
## t5 10.6500000 -0.2833333
## t6 -0.7466667  3.7233333
## t7 11.6266667  0.5283333
## t8  0.5283333  3.8566667
## |S|:  8313.241 
## 
## Control limits:
##  LCL      UCL
##    0 17.41705

Generate new “in control” data:

boilerNew = MASS::mvrnorm(10, mu = q1$center, Sigma = q1$cov)
q2 = mqcc(boiler, type = "T2.single", confidence.level = 0.999, 
          newdata = boilerNew, pred.limits = TRUE)

summary(q2)
## 
## Call:
## mqcc(data = boiler, type = "T2.single", pred.limits = TRUE, newdata = boilerNew,     confidence.level = 0.999)
## 
## T2.single chart for boiler 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  1.316342  5.305689  7.074224  7.680000  9.775744 17.575293 
## 
## Number of variables:  8
## Number of groups:  25
## Group sample size:  1
## 
## Center: 
##     t1     t2     t3     t4     t5     t6     t7     t8 
## 525.00 513.56 538.92 521.68 503.80 512.44 478.72 477.24 
## 
## Covariance matrix:
##            t1        t2        t3         t4         t5         t6
## t1 54.0000000 0.9583333 20.583333 31.2916667 20.3333333 -2.2916667
## t2  0.9583333 4.8400000  2.963333  2.6866667  0.3250000  3.0766667
## t3 20.5833333 2.9633333 22.993333 10.0566667  4.9833333  2.0366667
## t4 31.2916667 2.6866667 10.056667 22.3100000 13.5583333 -2.2700000
## t5 20.3333333 0.3250000  4.983333 13.5583333 11.4166667 -1.6583333
## t6 -2.2916667 3.0766667  2.036667 -2.2700000 -1.6583333  4.5066667
## t7 20.3750000 0.7050000  6.476667 12.6983333 10.6500000 -0.7466667
## t8  0.2083333 3.4433333  2.770000  0.1633333 -0.2833333  3.7233333
##            t7         t8
## t1 20.3750000  0.2083333
## t2  0.7050000  3.4433333
## t3  6.4766667  2.7700000
## t4 12.6983333  0.1633333
## t5 10.6500000 -0.2833333
## t6 -0.7466667  3.7233333
## t7 11.6266667  0.5283333
## t8  0.5283333  3.8566667
## |S|:  8313.241 
## 
## Summary of group statistics in boilerNew:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  2.387716  6.055542  8.214904  8.262216 11.238031 13.329123 
## 
## Number of groups:  10
## Group sample size:  1
## 
## Control limits:
##  LCL      UCL
##    0 17.41705
## 
## Prediction limits:
##  LPL      UPL
##    0 70.02943

Generate new “out of control” data:

boilerNew = MASS::mvrnorm(10, mu = 1.01*q1$center, Sigma = q1$cov)
q3 = mqcc(boiler, type = "T2.single", confidence.level = 0.999, 
          newdata = boilerNew, pred.limits = TRUE)

summary(q3)
## 
## Call:
## mqcc(data = boiler, type = "T2.single", pred.limits = TRUE, newdata = boilerNew,     confidence.level = 0.999)
## 
## T2.single chart for boiler 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  1.316342  5.305689  7.074224  7.680000  9.775744 17.575293 
## 
## Number of variables:  8
## Number of groups:  25
## Group sample size:  1
## 
## Center: 
##     t1     t2     t3     t4     t5     t6     t7     t8 
## 525.00 513.56 538.92 521.68 503.80 512.44 478.72 477.24 
## 
## Covariance matrix:
##            t1        t2        t3         t4         t5         t6
## t1 54.0000000 0.9583333 20.583333 31.2916667 20.3333333 -2.2916667
## t2  0.9583333 4.8400000  2.963333  2.6866667  0.3250000  3.0766667
## t3 20.5833333 2.9633333 22.993333 10.0566667  4.9833333  2.0366667
## t4 31.2916667 2.6866667 10.056667 22.3100000 13.5583333 -2.2700000
## t5 20.3333333 0.3250000  4.983333 13.5583333 11.4166667 -1.6583333
## t6 -2.2916667 3.0766667  2.036667 -2.2700000 -1.6583333  4.5066667
## t7 20.3750000 0.7050000  6.476667 12.6983333 10.6500000 -0.7466667
## t8  0.2083333 3.4433333  2.770000  0.1633333 -0.2833333  3.7233333
##            t7         t8
## t1 20.3750000  0.2083333
## t2  0.7050000  3.4433333
## t3  6.4766667  2.7700000
## t4 12.6983333  0.1633333
## t5 10.6500000 -0.2833333
## t6 -0.7466667  3.7233333
## t7 11.6266667  0.5283333
## t8  0.5283333  3.8566667
## |S|:  8313.241 
## 
## Summary of group statistics in boilerNew:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  6.66978 13.74057 18.72251 18.54601 20.76794 37.98010 
## 
## Number of groups:  10
## Group sample size:  1
## 
## Control limits:
##  LCL      UCL
##    0 17.41705
## 
## Prediction limits:
##  LPL      UPL
##    0 70.02943

Recompute by providing “robust” estimates for the means and the covariance matrix:

rob = MASS::cov.rob(boiler)
q4 = mqcc(boiler, type = "T2.single", center = rob$center, cov = rob$cov)

summary(q4)
## 
## Call:
## mqcc(data = boiler, type = "T2.single", center = rob$center,     cov = rob$cov)
## 
## T2.single chart for boiler 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##   1.13929   5.35641   7.52501  14.21985  10.68735 117.73932 
## 
## Number of variables:  8
## Number of groups:  25
## Group sample size:  1
## 
## Center: 
##       t1       t2       t3       t4       t5       t6       t7       t8 
## 525.4348 513.4348 539.9130 521.6087 503.8261 512.4783 478.8696 477.2609 
## 
## Covariance matrix:
##            t1        t2        t3         t4         t5         t6
## t1 41.0750988 2.9387352 16.221344 26.8596838 16.7154150 -2.7173913
## t2  2.9387352 4.9841897  4.903162  3.4051383  0.8063241  3.4189723
## t3 16.2213439 4.9031621 12.173913 11.6007905  4.8932806  1.7252964
## t4 26.8596838 3.4051383 11.600791 20.4308300 12.1561265 -2.4407115
## t5 16.7154150 0.8063241  4.893281 12.1561265 10.6047431 -1.8221344
## t6 -2.7173913 3.4189723  1.725296 -2.4407115 -1.8221344  4.8972332
## t7 15.4683794 1.4683794  4.897233 11.0375494  9.5217391 -0.8893281
## t8  0.1086957 3.7905138  2.750988  0.1976285 -0.3162055  4.0513834
##            t7         t8
## t1 15.4683794  0.1086957
## t2  1.4683794  3.7905138
## t3  4.8972332  2.7509881
## t4 11.0375494  0.1976285
## t5  9.5217391 -0.3162055
## t6 -0.8893281  4.0513834
## t7 10.1185771  0.5355731
## t8  0.5355731  4.2015810
## |S|:  1071.809 
## 
## Control limits:
##  LCL      UCL
##    0 14.26225

Pareto chart

defect = c(80, 27, 66, 94, 33)
names(defect) = c("price code", "schedule date", "supplier code", "contact num.", "part num.")
pareto.chart(defect, ylab = "Error frequency")

##                
## Pareto chart analysis for defect
##                 Frequency Cum.Freq. Percentage Cum.Percent.
##   contact num.   94.00000  94.00000   31.33333     31.33333
##   price code     80.00000 174.00000   26.66667     58.00000
##   supplier code  66.00000 240.00000   22.00000     80.00000
##   part num.      33.00000 273.00000   11.00000     91.00000
##   schedule date  27.00000 300.00000    9.00000    100.00000

Cause and effect diagram

cause.and.effect(cause = list(Measurements = c("Micrometers", "Microscopes", "Inspectors"),
                              Materials    = c("Alloys", "Lubricants", "Suppliers"),
                              Personnel    = c("Shifts", "Supervisors", "Training", "Operators"),
                              Environment  = c("Condensation", "Moisture"),
                              Methods      = c("Brake", "Engager", "Angle"),
                              Machines     = c("Speed", "Lathes", "Bits", "Sockets")),
                 effect = "Surface Flaws")

Process variation examples

In the following simulated data are used to describe some models for process variation. For further details see Wetherill, G.B. and Brown, D.W. (1991) Statistical Process Control, New York, Chapman and Hall, Chapter 3.

Simple random variation

\(x_{ij} = \mu + \sigma_W \epsilon_{ij}\)

mu = 100
sigma_W = 10
epsilon = rnorm(500)
x = matrix(mu + sigma_W*epsilon, ncol=10, byrow=TRUE)
q = qcc(x, type="xbar")

q = qcc(x, type="R")

q = qcc(x, type="S")

Between and within sample extra variation

\(x_{ij} = \mu + \sigma_B u_i + \sigma_W \epsilon_{ij}\)

mu = 100
sigma_W = 10
sigma_B = 5
epsilon = rnorm(500)
u = as.vector(sapply(rnorm(50), rep, 10))
x = mu + sigma_B*u + sigma_W*epsilon
x = matrix(x, ncol=10, byrow=TRUE)
q = qcc(x, type="xbar")

q = qcc(x, type="R")

q = qcc(x, type="S")

Autocorrelation

\(x_{ij} = \mu + W_i + \sigma_W \epsilon_{ij}\)
where \(W_i = \rho W_{i-1} + \sigma_B u_i = \sigma_B u_i + \rho \sigma_B u_{i-1} + \rho^2 \sigma_B u_{i-2} + \ldots\),
and \(W_0 = 0\).

mu = 100
rho = 0.8
sigma_W = 10
sigma_B = 5
epsilon = rnorm(500)
u = rnorm(500)
W = rep(0,100)
for(i in 2:length(W))
    W[i] = rho*W[i-1] + sigma_B*u[i]
x = mu + sigma_B*u + sigma_W*epsilon
x = matrix(x, ncol=10, byrow=TRUE)
q = qcc(x, type="xbar")

q = qcc(x, type="R")

q = qcc(x, type="S")

Recurring cycles

Assume we have 3 working turns of 8 hours each for each working day, so \(8 \times 3 = 24\) points in time, and at each point we sample 5 units.

\(x_{ij} = \mu + W_i + \sigma_W \epsilon_{ij}\)
where \(W_i\) (\(i=1,\ldots,8\)) is the cycle.

mu = 100
sigma_W = 10
epsilon = rnorm(120, sd=0.3)
W = c(-4, 0, 1, 2, 4, 2, 0, -2) # assumed workers cycle
W = rep(rep(W, rep(5,8)), 3)
x = mu + W + sigma_W*epsilon
x = matrix(x, ncol=5, byrow=TRUE)
q = qcc(x, type="xbar")

q = qcc(x, type="R")

q = qcc(x, type="S")

Mixture

\(x_{ij} = \mu_1 p + \mu_2 (1-p) + \sigma_W \epsilon_{ij}\)
where \(p = \Pr(\text{Process #1})\).

mu1 = 90
mu2 = 110
sigma_W = 10
epsilon = rnorm(500)
p = rbinom(50, 1, 0.5)
mu = mu1*p + mu2*(1-p)
x = rep(mu, rep(10, length(mu))) + sigma_W*epsilon
x = matrix(x, ncol=10, byrow=TRUE)
q = qcc(x, type="xbar")

q = qcc(x, type="R")

q = qcc(x, type="S")

Sudden jumps

\(x_{ij} = \mu_i + \sigma_W \epsilon_{ij}\)
where \(\mu_i\) is the mean of the process for state \(i\) (\(i=1,\ldots,k)\).

mu = rep(c(95,110,100,90), c(20,35,25,20))
sigma_W = 10
epsilon = rnorm(500)
x = rep(mu, rep(5, length(mu))) + sigma_W*epsilon
x = matrix(x, ncol=10, byrow=TRUE)
q = qcc(x, type="xbar")

q = qcc(x, type="R")

q = qcc(x, type="S")