Rohit GitHub
Statistical Quality Control (SQC) using “R”
1. 99.73% observations to be between +(-)3 Sigma / Six Sigma .
2. Control Charts - UCL , Central Line and LCL , leads to outlier or Special cause detection.
3. Using package - “qcc” to create a Control Chart within “R”
4. Multivariate Control Charts - Control Ellipsoid - http://arxiv.org/ftp/arxiv/papers/0901/0901.2880.pdf
library(MSQC)
library("MSQC", lib.loc="~/R/win-library/3.1")
data(package="MSQC")
library(qcc)
citation("qcc")To cite qcc in publications use:
Scrucca, L. (2004). qcc: an R package for quality control charting and statistical process control. R News 4/1, 11-17.
A BibTeX entry for LaTeX users is
@Article{, title = {qcc: an R package for quality control charting and statistical process control}, author = {Luca Scrucca}, journal = {R News}, year = {2004}, pages = {11–17}, volume = {4/1}, url = {http://CRAN.R-project.org/doc/Rnews/}, }
#Scrucca, L. (2004). qcc: an R package for quality control
#charting and statistical process control. R News 4/1, 11-17.
set.seed(20)
x <-round(rnorm(120,20,2),2)
length <-matrix(x, ncol=4, byrow=TRUE)
par(mfrow=c(1,2))
qcc(length, type="xbar", std.dev="RMSDF"); qcc(length, type="R")
List of 11 $ call : language qcc(data = length, type = “xbar”, std.dev = “RMSDF”) $ type : chr “xbar” $ data.name : chr “length” $ data : num [1:30, 1:4] 22.3 19.1 19.1 18.7 21.9 … ..- attr(, “dimnames”)=List of 2 $ statistics: Named num [1:30] 20.5 18.2 19.4 19.4 20.7 … ..- attr(, “names”)= chr [1:30] “1” “2” “3” “4” … $ sizes : int [1:30] 4 4 4 4 4 4 4 4 4 4 … $ center : num 20 $ std.dev : num 2.03 $ nsigmas : num 3 $ limits : num [1, 1:2] 17 23.1 ..- attr(, “dimnames”)=List of 2 $ violations:List of 2 - attr(, “class”)= chr “qcc” 
List of 11 $ call : language qcc(data = length, type = “R”) $ type : chr “R” $ data.name : chr “length” $ data : num [1:30, 1:4] 22.3 19.1 19.1 18.7 21.9 … ..- attr(, “dimnames”)=List of 2 $ statistics: Named num [1:30] 6.24 6.92 1.07 5.69 2.11 … ..- attr(, “names”)= chr [1:30] “1” “2” “3” “4” … $ sizes : int [1:30] 4 4 4 4 4 4 4 4 4 4 … $ center : num 4.07 $ std.dev : num 1.98 $ nsigmas : num 3 $ limits : num [1, 1:2] 0 9.28 ..- attr(, “dimnames”)=List of 2 $ violations:List of 2 - attr(, “class”)= chr “qcc”
#
cap<-qcc(length, type="xbar", nsigmas=3, plot=T)
process.capability(cap, spec.limits=c(14,26))
Process Capability Analysis
Call: process.capability(object = cap, spec.limits = c(14, 26))
Number of obs = 120 Target = 20 Center = 20.01 LSL = 14 StdDev = 1.976 USL = 26
Capability indices:
Value 2.5% 97.5%Cp 1.012 0.8837 1.141 Cp_l 1.014 0.8952 1.134 Cp_u 1.010 0.8913 1.129 Cp_k 1.010 0.8685 1.152 Cpm 1.012 0.8842 1.140
Exp
#
mu <-c(0,0) # The Mean vector - "mu" of a Bivariate Normal Distribution
sigma <- matrix(c(10,3,3,6),2,2) # Sigma the Covariance Matrix
# The Square Covariance Matrix structure needs to be seen and understood ...
sigma [,1] [,2][1,] 10 3 [2,] 3 6
#
sigma[1,1][1] 10
# Which is the Observation "10" t location - Row - 1 and Column -1, within the Square Matrix "sigma"
sigma[2,2][1] 6
# Which is the Observation "6" at location - Row - 2 and Column -2 , within the Square Matrix "sigma"
#
rho <- sigma[1,2]/(sqrt(sigma[1,1]*sigma[2,2]))
rho[1] 0.3872983
# rho=3/(sqrt(10*6))
# dummy<-3/sqrt(60) - same as value of "rho"
# https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cor.html
#
var1<-seq(-12,12,.7)
var1[1] -12.0 -11.3 -10.6 -9.9 -9.2 -8.5 -7.8 -7.1 -6.4 -5.7 -5.0 [12] -4.3 -3.6 -2.9 -2.2 -1.5 -0.8 -0.1 0.6 1.3 2.0 2.7 [23] 3.4 4.1 4.8 5.5 6.2 6.9 7.6 8.3 9.0 9.7 10.4 [34] 11.1 11.8
#
var2<-var1
knitr::kable(var2)| -12.0 |
| -11.3 |
| -10.6 |
| -9.9 |
| -9.2 |
| -8.5 |
| -7.8 |
| -7.1 |
| -6.4 |
| -5.7 |
| -5.0 |
| -4.3 |
| -3.6 |
| -2.9 |
| -2.2 |
| -1.5 |
| -0.8 |
| -0.1 |
| 0.6 |
| 1.3 |
| 2.0 |
| 2.7 |
| 3.4 |
| 4.1 |
| 4.8 |
| 5.5 |
| 6.2 |
| 6.9 |
| 7.6 |
| 8.3 |
| 9.0 |
| 9.7 |
| 10.4 |
| 11.1 |
| 11.8 |
#
f<-matrix(0, length(var1), length(var1))
for(i in 1:length(var1))
{
for(j in 1:length(var1))
{
f[i,j]<-1/(2*pi*sqrt(sigma[1,1]*sigma[2,2]*(1-rho^2)))*exp(-1/(2*(1-rho^2)) * ((var1[i]-mu[1])^2/sigma[1,1] + (var2[j] - mu[2])^2/sigma[2,2]-2 *rho*((var1[i] - mu[1])*(var2[j]-mu[2]))/(sqrt(sigma[1,1])*sqrt(sigma[2,2]))))}}
#
str(f)num [1:35, 1:35] 1.65e-08 2.62e-08 3.94e-08 5.60e-08 7.50e-08 …
#
persp(var1,var2,f,xlab="Var-1",ylab="Var-2",zlab="f(var1,var2)",
theta=30,phi=30,col = "red",shade = 0.05,r=50,d=0.02,expand=0.9,ltheta=90,lphi=180, nticks=4,box=F)
#
jet.colors <- colorRampPalette( c("pink", "red") )
#
nbcol <- 100
color <- jet.colors(nbcol)
persp(var1,var2,f,xlab="Var-1",ylab="Var-2",zlab="f(var1,var2)",
theta=30,phi=30,col = color,shade = 0.05,r=50,d=0.02,expand=0.9,ltheta=90,lphi=180, nticks=4,box=T)
# Further Shading - using Code from help -persp {graphics}
z <- outer(var1, var2, function(a, b) a*b^2)
nrz <- nrow(z)
ncz <- ncol(z)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
persp(var1,var2,f,xlab="Var-1",ylab="Var-2",zlab="f(var1,var2)",
theta=30,phi=30,col = color[facetcol],shade = 0.05,r=50,d=0.02,expand=0.9,ltheta=90,lphi=180, nticks=4,box=T)
# In the Above seen Perspective Plot - persp {graphics}
# f == the matrix containing values to be plotted.# R Six Sigma 3D Scatter Plots / Graphs - using data set "dowel1" and dowel2"
# from package / library(MSQC)
# data(dowel1)
# View(dowel1)
# ?dowel1
library(scatterplot3d)
library("scatterplot3d", lib.loc="~/R/win-library/3.1")
# Sample 3D Scatter Plot
attach(mtcars)
View(mtcars)
scatterplot3d(x=mtcars$wt,
y=mtcars$disp,
z=mtcars$mpg,angle=10,color=rgb(.5, .5, .5, .5),pch=20,grid=T,highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue") # as we have used - highlight.3d=TRUE , the - color=rgb(.5, .5, .5, .5) , is ignored but still retained in the code. ## Warning in scatterplot3d(x = mtcars$wt, y = mtcars$disp, z = mtcars$mpg, :
## color is ignored when highlight.3d = TRUE
scatterplot3d(x=mtcars$wt,
y=mtcars$disp,
z=mtcars$mpg,angle=30,color=rgb(.5, .5, .5, .5),pch=20,grid=T,highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue")## Warning in scatterplot3d(x = mtcars$wt, y = mtcars$disp, z = mtcars$mpg, :
## color is ignored when highlight.3d = TRUE
#
scatterplot3d(x=mtcars$wt,
y=mtcars$disp,
z=mtcars$mpg,angle=50,color=rgb(.5, .5, .5, .5),pch=20,grid=T,highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue")## Warning in scatterplot3d(x = mtcars$wt, y = mtcars$disp, z = mtcars$mpg, :
## color is ignored when highlight.3d = TRUE
#
# Now plotting with own data - S=Qlty SCORE, L=Location of Employee , D= Department of Employee
#
dQlt <- read.csv("C:/STAT/SQC_Stat_Qlty_Ctrl/dQlt.csv")
summary(dQlt) S L D Min. :-10.000 Min. :1.000 Min. :1.000
1st Qu.: -1.750 1st Qu.:2.000 1st Qu.:1.000
Median : 8.000 Median :4.000 Median :2.000
Mean : 8.604 Mean :3.604 Mean :2.412
3rd Qu.: 18.000 3rd Qu.:5.000 3rd Qu.:3.000
Max. : 30.000 Max. :6.000 Max. :4.000
#
#
s3D_dQlt<-scatterplot3d(x=dQlt$S,
y=dQlt$D,
z=dQlt$L,angle=30,color=rgb(.5, .5, .5, .5),pch=20,grid=T,highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue")## Warning in scatterplot3d(x = dQlt$S, y = dQlt$D, z = dQlt$L, angle = 30, :
## color is ignored when highlight.3d = TRUE# As seen above - Random Scatter as No Relation Amongst variables ...
s3D_dQlt<-scatterplot3d(x=dQlt$S,
y=dQlt$D,
z=dQlt$L,angle=30,color=rgb(.5, .5, .5, .5),pch=20,grid=T,highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue")## Warning in scatterplot3d(x = dQlt$S, y = dQlt$D, z = dQlt$L, angle = 30, :
## color is ignored when highlight.3d = TRUE
#
# ,angle=30, scale.y=0.45, type="n",
# pch=16, color=rgb(0, 0, 0, .5),
# x.ticklabs=pretty(LONG, 3),
# grid=T, zlim=c(-20, 0))
# par(lab=c(3, 3, 0))
# s3d <- with(dowel1,scatterplot3d(diameter,
# angle=30, scale.y=0.45, type="n",
# pch=16, color=rgb(0, 0, 0, .5),
# x.ticklabs=pretty(LONG, 3),
# grid=FALSE, zlim=c(-20, 0)))
# Code below this point is for Testin need not be Run
# Experiment -1
zz <- seq(-10, 10,length.out=1000)
# zz
xx<-cos(zz)
yy<-sin(zz)
scatterplot3d(xx, yy, zz, highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue", main="Scatter-plot3d-1 [ Data Length =1000 Obs.] ", cex.symbols=2,pch=20)
# Experiment -2
zz2 <- seq(-10, 10,length.out=1000)
# zz
xx2<-1-(zz2)
yy2<-10-(zz2)
scatterplot3d(xx2, yy2, zz2, highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue", main="Scatter-plot3d-1 [ Data Length =1000 Obs.] ", cex.symbols=2,pch=20)
# R Code Ends here ....
## Contact… - email - Github - LinkedIn - Twitter
## License - GPL - No Copyright (http://www.gnu.org/licenses/licenses.en.html#GPL)