Quality Control-2
The ASTA team
OC curves
Usual setup:
- \(m\) samples with sample size \(n\).
- Process mean \(\mu\) and standard deviation \(\sigma\).
- Sample means have standard deviation \(\frac{\sigma}{\sqrt{n}}\).
- By default we use the 3*sigma rule.
Assume that the mean is shifted by \(c\times\sigma\).
What is the probability of NOT getting an immediate alarm?
This is also called a type II error.
library(qcc)
data(pistonrings)
diam <- pistonrings$diameter
phaseI <- matrix(diam[1:125],25,byrow=TRUE)
phaseII <- matrix(diam[126:200],15,byrow=TRUE)
xbcc <- qcc(phaseI, std.dev = "UWAVE-SD", type = "xbar", plot = FALSE,
newdata = phaseII, title = "qq chart: pistonrings")
OC curves - example
For the actual sample size (n=5) and in case of a process shift c=2, the error probability is around 7%.
If we increase sample size to n=10 then we are almost sure to immediate detection of a shift by \(2\sigma\).
CUSUM chart
The CUSUM chart is generally better than the xbar chart for detecting small shifts in the mean of a process.
Consider the standardized residuals \[z_i=\frac{\sqrt{n}(x_i-\hat{\mu}_0)}{\hat{\sigma}}\] where
- \(\hat{\mu}_0\) is an estimate of the process mean.
- \(\hat{\sigma}\) is an estimate of the process standard deviation.
For an in-control proces the cumulative sum (CUSUM) of residuals should vary around \(zero\).
The CUSUM chart is developed to see, if there is a drift away from zero. To that end we define 2 processes controlling for downward(D) respectively upward(U) drift:
\begin{align*}
D(i)&=\max\{0,D(i-1)-k-z_i\}\\
U(i)&=\max\{0,U(i-1)+z_i-k\}
\end{align*}
where \(k\) is a positive number.
Interpretation of CUSUM chart
Interpretation of
\begin{align*}
D(i)&=\max\{0,D(i-1)-k-z_i\}\\
U(i)&=\max\{0,U(i-1)+z_i-k\}
\end{align*}
- If \(z_i< -k\), i.e. \(z_i\) is more that \(k\) below zero, then \(D\) is increased. If this happens a number of times in a row, then \(D\) grows “big”.
- If \(z_i>k\), i.e. \(z_i\) is more that \(k\) above zero, then \(U\) is increased. If this happens a number of times in a row, then \(U\) grows “big”.
The process is considered out of control if \(D\) or \(U\) exceeds a limit \(h\).
In qcc \(k\) and \(h\) is specified by the arguments:
se.shift is \(k\), which has a default value of 1.decision.interval is \(h\), which has a default value of 5.
CUSUM chart example
h <- cusum(phaseI, newdata = phaseII, title = "CUSUM chart: pistonrings")
The chart includes a plot of
- U controlling for positive drift, which is clearly happening in phase II.
- D controlling for negative drift.
EWMA chart
The Exponentially Weighted Moving Average (EWMA) is a statistic for monitoring the process, which averages the data in a way that gives most weight to recent data.
The EWMA is formally defined by \[M_t=\lambda x_t+(1-\lambda)M_{t-1},\ \ t=1,2,\ldots,T\] where
- \(M_0\) is the mean of some historical data.
- \(x_t\) is the measurement at time \(t\).
- \(T\) is the length of the sampling period.
- \(0 < \lambda \leq 1\) is a smoothing parameter, where \(\lambda=1\) corresponds to “no memory”.
The influence of \(x_t\) on \(M_{t+s}\) is of the order \((1-\lambda)^s\), i.e. exponentially decreasing, which explains the term “Exponentially Weighted”.
EWMA chart
Estimated variance for the EWMA process
- \(s_M^2=\frac{\lambda}{2-\lambda}s^2\)
- \(s\) is the standard deviation from historical data.
\[
\begin{aligned}
\text{UCL:} \quad &M_0+k s_M\\
\text{CL :} \quad &M_0\\
\text{LCL:} \quad &M_0-k s_M
\end{aligned}
\]
Conventional choise of \(k\) is 3.
EWMA chart example
h <- ewma(phaseI, newdata = phaseII, title = "EWMA chart: pistonrings")
- Observed values are “plus” and predicted values are “dots”.
- Default value of smoothing is \(0.2\). May be set by the argument
lambda.
- Default value of \(k\) is \(3\). May be set by the argument
nsigmas.
Multivariate charts
- In some situations , it is important to watch the covariation between two or more variables.
- We shall limit our considerations to two variables \(y\) and \(x\).
- We assume a linear regression equation for \(y\) depending on \(x\).
- We have data to estimate the line and the expected deviation from the line.
- We have estimates of mean and standard deviation for \(x\).
An outlier would deviate from the line and/or the x-mean, so we calculate
- \(Zy_x\) denoting the standardized residual from the regression line.
- \(Zx\) denoting the standardized residual from the expected mean of x.
We expect that both \(Zy_x\) and \(Zx\) should be within the limits \(\pm2\). In order to get an overall teststatistic, we calculate \[T^2=\frac{m-1}{m-2}Zy_x^2+Zx^2\] where \(m\) is sample size.
Multivariate charts
Alternatively, one might interchange the role of \(x\) and \(y\). But that does not matter. Actually it holds that \[T^2=\frac{m-1}{m-2}Zy_x^2+Zx^2=\frac{m-1}{m-2}Zx_y^2+Zy^2\]
- If the process is in control, then \(T^2\) has a chi-square distribution with 2 degrees of freedom.
- It is critical to the process if \(T^2\) is large, i.e. we only need an upper limit.
load(url("https://asta.math.aau.dk/datasets?file=T2Example.RData"))
## [,1] [,2] [,3] [,4]
## [1,] 72 84 79 49
## [2,] 56 87 33 42
## [3,] 55 73 22 60
Multivariate chart example
## [,1] [,2] [,3] [,4]
## [1,] 23 30 28 10
## [2,] 14 31 8 9
## [3,] 13 22 6 16
X is a listX$X1 is a matrix, where the rows are samples of variable X1- Actual sample size is 4 and the number of samples is 20.
Similarly X$X2 has samples for variable X2
Multivariate chart example
h <- mqcc(X, type = "T2")
When the parameters are estimated from phase I samples, then the reference distribution is not chi-square, but rather a scaled F-distribution.
Multivariate chart example
ellipseChart(h, show.id = TRUE)
Observation 10 deviates a lot from the regression line.
The acceptance area is an ellipse.
Acceptance sampling
Set-up:
- Production proces where we produce lots of size \(N\).
- From the lot we take a sample of size \(n\).
- If the number of defective items exceeds the number \(c\), then the lot is rejected.
We term this a \((N,n,c)\) sampling plan. Possible reasons for acceptance sampling:
- Testing is destructive
- The cost of 100% inspection is very high
- 100% inspection takes too long
Sampling distributions
- Let \(X\) denote the number of defective items in the sample.
- Let \(p\) denote the fraction of defective items in the lot.
The correct sampling distribution of \(X\) is the so called hypergeometric distribution with parameters \((N,n,p)\).
If \(N>>n\), then the sampling distribution of \(X\) is well approximated by the simpler binomial distribution with parameters \((n,p)\).
If \(N>>n\) and \(p\) is smal, then the sampling distribution of \(X\) is well approximated by the much simpler poisson distribution with parameter \(np\).
OC curve of a sampling plan
For a given sampling plan \((N,n,c)\) the probability of accepting the lot depends on
- the fraction \(p\) of defective items in the lot
- the assumed sampling distribution
We can use the function OC2c in the package AcceptanceSampling to determine these probabilities.
Sampling plan: \((1000,100,2)\)
library(AcceptanceSampling)
OCbin <- OC2c(100, 2) #default binomial
OCpoi <- OC2c(100, 2, type = "poisson")
OChyp <- OC2c(100, 2, type = "hypergeom", N = 1000)
OC curve of a sampling plan
par(mfrow=c(2,2)) #division of plot window
plot(1:10, type = "n", axes = FALSE, xlab = "", ylab = "")
text(4, 5, "(N,n,c)\n(1000,100,2)", cex = 2)
xl <- c(0, 0.1)
plot(OChyp, xlim = xl, main = "hyper")
plot(OCbin, xlim = xl, main = "binom")
plot(OCpoi, xlim = xl, main = "poiss")
Find sampling plan
Suppose that \(N\) is fixed. If we specify 2 points on the OC curne, then this determines \((n,c)\).
PRP: Producer Risk Point with coordinates (p1,q1): p1 is a fraction of defectives, that the producer finds acceptable, e.g. p1=0.01. The corresponding probability q1 of accept should then be high, e.g. q1=0.95.
CRP: Consumer Risk Point with coordinates (p2,q2): p2\(>\)p1 is a fraction of defectives, that the consumer finds unacceptable, e.g. p2=0.05. The corresponding probability q2 of accept should then be low,e.g. q2=0.01.
find.plan(PRP = c(.01,.95), CRP = c(.05,.01))[1:2]
## $n
## [1] 259
##
## $c
## [1] 5
Default assumption is binomial sampling.
Find sampling plan
plan <- find.plan(c(.01,.95), c(.05,.01), type = "hyp", N = 200)[1:2]
plan
## $n
## [1] 121
##
## $c
## [1] 2
OChyp <- OC2c(plan$n, plan$c, type = "hyp", N = 200, pd = c(.01,.05))
attr(OChyp, "paccept")
## [1] 1.000000000 0.009609746
We cannot have an exact match of the required values \((0.95,0.01)\) since \(n\) and \(c\) must be integers.
Double sampling
Let \(0\leq c_1 < r_1\) be integers.
Furthermore, \(c_2\) is an integer such that \(c_1<c_2\).
- \(x_1\): number of defectives in an initial sample of size \(n_1\).
- If \(x_1\leq c_1\) accept the lot.
- If \(r_1\leq x_1\) reject the lot.
- If \(c_1< x_1<r_1\): Take a second sample of size \(n_2\) and let \(x_2\) be the number of defectives.
- If \(x_1+x_2\leq c_2\) accept the lot. Otherwise reject.
This is known as a double sampling plan.
OC curve of a double sampling plan
Determining the OC curve of a double sampling plan requires input of \(n=c(n_1,n_2)\), \(c=c(c_1,c_2)\) and \(r=c(r_1,r_2)\), where \(r_2=c_2+1\).
x <- OC2c(c(125,125), c(1,4), c(4,5), pd = seq(0,0.1,0.001))
x
## Acceptance Sampling Plan (binomial)
##
## Sample 1 Sample 2
## Sample size(s) 125 125
## Acc. Number(s) 1 4
## Rej. Number(s) 4 5
OC curve of a double sampling plan
