Process Control was one of the major keys to the to the Japanese boom in the 70’s and 80’s. By the 90’s the U.S. had become believers as well. However, by 2005 many of the nation’s companies had taken process control too far. In fact Forbes Magazine just recently put out another article questioning over emphasis on process control.
As good practitioners we understand that there are times to get our process an in control and there are times to be innovative and we are careful not to let one objective step on the other. We will be using the qcc package in R. In 2004 there was a good article on the package. Akritas has provided all the example functions we will use.
After loading the two function below (click the code button) into R, run this line of code pdr = plot_process(process_data(step_shift=-1)). The plot_process function plots the xbar control chart step-by-step as if you were on the factory floor. The process_data function generates random process data. When the step_shift variable is negative a random step in the process is selected and the prescribed changes in the process mean and/or sd are implemented. Your goal is to use the chart and the process rules to decide when the process is out of control.
step_shift that prints to screen. Report how different your process was from the actual in the file name (e.g. Rplot_-4.png or Rplot_8.png) of the process control chart that you upload.The above example was based on a sample size of 5 per time step. Let’s see how well you can do with a xbar.one type. This means that only one sample is observed at each step. Use pdr = plot_process_one(process_data(step_shift=-1)) and attempt the same process as above.
#install.packages("qcc")
library(qcc)
process_data = function(mean=5,sd=.2,shift_mean=.1,shift_sd=0,step_shift=15,num_steps=25,num_obs=5,seed_number=NA){
if (is.na(seed_number)==FALSE) set.seed(seed_number)
if (step_shift==0){
start_values = NULL
end_values = t(replicate(num_steps,rnorm(num_obs,mean=mean,sd=sd)))
}
if (step_shift>0){
start_values = t(replicate(step_shift,rnorm(num_obs,mean=mean,sd=sd)))
end_values = t(replicate(num_steps-step_shift,rnorm(num_obs,mean=mean+shift_mean,sd=sd+shift_sd)))
}
if (step_shift<0){
step_shift = sample(abs(num_steps-10),1)
start_values = t(replicate(step_shift,rnorm(num_obs,mean=mean,sd=sd)))
end_values = t(replicate(num_steps-step_shift,rnorm(num_obs,mean=mean+shift_mean,sd=sd+shift_sd)))
}
list(data=rbind(start_values,end_values),input=c("mean"=mean,"sd"=sd,"shift_mean"=shift_mean,
"shift_sd"=shift_sd,"step_shift"=step_shift,
"num_steps"=num_steps,"num_obs"=num_obs,"seed_number"=seed_number))
}
plot_process = function(x,mean=5,sd=2){
for (i in 2:nrow(x$data)){
qcc(x$data[1:i,],type="xbar",center=mean,sd=sd)$call
cat ("Press [enter] to continue")
line <- readline()
}
print(x$input)
return(x)
}
plot_process_one = function(x,mean=5,sd=2){
for (i in 2:nrow(x$data)){
qcc(x$data[1:i,1],type="xbar.one",center=mean,sd=sd)$call
cat ("Press [enter] to continue")
line <- readline()
}
print(x$input)
return(x)
}
$$\hat{\sigma}_2=\frac{1}{B_n}\frac{1}{k}\sum_{i=1}^kr_i=\frac{1}{B_n}\bar{r}$$
$$\hat{\sigma}_1=\frac{1}{A_n}\frac{1}{k}\sum_{i=1}^ks_i=\frac{1}{A_n}\bar{s}$$
Here is a table for the constants An and Bn up to n = 20
library(DT)
library(qcc)
library(ggplot2)
Bn = c(1.128,1.693,2.059,2.326,2.534,2.704,2.847,2.97,3.078,3.173,3.258,3.336,3.407,3.472,3.532,3.588,3.64,3.689,3.735)
An = c(0.7979,0.8862,0.9213,0.94,0.9515,0.9594,0.965,0.9693,0.9727,0.9754,0.9776,0.9794,0.981,0.9823,0.9835,0.9845,0.9854,0.9862,0.9869)
n = 2:20datatable(data.frame(n=n,An=An,Bn=Bn),rownames=FALSE)pd = process_data(mean=20,num_steps=30,step_shift=15,shift_mean=-2,shift_sd=0,sd=3)$data
qcc(pd[1:15,],type="xbar",newdata=pd[16:nrow(pd),])$call
pdr = plot_process(process_data(step_shift=-1,seed_number=5))
pd_sd = process_data(mean=20,num_steps=30,step_shift=0,shift_mean=0,shift_sd=0,sd=3)$data
(qcc(pd_sd,type="xbar",std.dev="UWAVE-SD")$limits[2] - qcc(pd_sd,type="xbar",std.dev="UWAVE-SD")$center)*(sqrt(ncol(pd_sd))/3)
(qcc(pd_sd,type="xbar",std.dev="UWAVE-R")$limits[2] - qcc(pd_sd,type="xbar",std.dev="UWAVE-R")$center)*(sqrt(ncol(pd_sd))/3)
sigma2 = sd.R(pd_sd,std.dev="UWAVE-R")
sigma1 = sd.S(pd_sd,std.dev="UWAVE-SD")
raw_sigma = mean(apply(pd_sd,1,sd))pd = process_data(mean=3.555,num_steps=21,step_shift=9,shift_mean=.042,shift_sd=0,sd=.075)$data
qcc(pd[1:8,],type="xbar",newdata=pd[9:nrow(pd),])$call
pdr = plot_process(process_data(step_shift=-1,seed_number=5))nistHH = read.csv("http://www.nist.gov/sites/default/files/documents/mml/acmd/structural_materials/HH-105-48877.csv",stringsAsFactors = FALSE,skip = 2)
data.dir = "http://media.pearsoncmg.com/cmg/pmmg_mml_shared/mathstatsresources/Akritas"
shaft = read.table(file.path(data.dir,"CcShaftDM.txt"),header=T)
data(pistonrings)
pistons = qcc.groups(pistonrings$diameter,pistonrings$sample)
qcc(shaft,type="xbar")
qcc(pistons[1:25,],type="xbar",newdata=pistons[26:40,])
data("orangejuice")
orangejuice$d <- orangejuice$D/orangejuice$size
boxplot(d ~ trial,data=orangejuice)
mark <- ifelse(orangejuice$trial, 1, 2)
plot(orangejuice$sample, orangejuice$d, type="b", col=mark, pch=mark)
qcc(subset(nistHH,n==5)[,4:8],type="xbar")
qcc(subset(nistHH,n==5)[,4:8],type="R")
qcc(subset(nistHH,n==5)[,4:8],type="S")