qcr: Quality Control Review

#-----------------------------------------------------------------------------#
#                                                                             #
#                  QUALITY CONTROL STATISTICS IN R                            #
#                                                                             #
#  An R package for statistical in-line quality control.                      #
#                                                                             #
#  Written by: Miguel A. Flores Sanchez                                       #
#              Professor of the Mathematics Department                        #
#              Escuela Politecnica Nacional, Ecuador                          #
#              miguel.flores@epn.edu.ec                                       #
#                                                                             #
#-----------------------------------------------------------------------------#
#-----------------------------------------------------------------------------#
# Main function to create a 'qcs.ca' object
#-----------------------------------------------------------------------------#
##' Capability Analysis
##' 
##' Calculates the process capability indices Cp, Cpk, Cpl, Cpu, Cpm, Cpmk for a 
##' qcs object and normal distribution. Also, this function calculates confidence
##' limits for \eqn{C_p}{C_p} using the method described by Chou et al. (1990). 
##' Approximate confidence limits for \eqn{C_{pl}}{C_pl}, \eqn{C_{pu}}{C_pu} and
##' \eqn{C_{pk}}{C_pk} are computed using the method in Bissell (1990). 
##' Confidence limits for \eqn{C_{pm}}{C_pm} are based on the method of Boyles (1991); 
##' this method is approximate and it assumes the target is midway between the 
##' specification limits.
##' Moreover, calculates the process capability indices CNp, CNpk, CNpm, CNpmk for a qcs object. 
##' A histogram with a density curve is displayed along with the specification limits, a 
##' Quantile-Quantile Plot for the specified distribution and contour graph is plotted 
##' for estimate the index Cpm.
##' @aliases qcs.ca 
##' @param object qcs object of type \code{"qcs.xbar"} or \code{"qcs.one"}.
##' @param limits A vector specifying the lower and upper specification limits.
##' @param target A value specifying the target of the process. 
##' If it is \code{NULL}, the target is set at the middle value between specification limits.
##' @param std.dev A value specifying the within-group standard deviation.
##' @param nsigmas A numeric value specifying the number of sigmas to use.
##' @param confidence A numeric value between 0 and 1 specifying the probabilities 
##' for computing the quantiles.
##' This values is used only when \code{object} values is provided. 
##' By default \code{confidence}=0.9973.
##' @param plot Logical value indicating whether graph should be plotted.
##' @param main Title of the plot.
##' @param ... Arguments to be passed to or from methods.
##' @export
##' @references 
##' Montgomery, D.C. (1991) \emph{Introduction to Statistical Quality Control}, 2nd
##' ed, New York, John Wiley & Sons. \cr
##' Tong, L.I. and Chen, J.P. (1998), \emph{Lower con???dence limits of process capability 
##' indices for nonnormal process distributions.} International Journal of Quality & Reliability Management, 
##' Vol. 15 No. 8/9, pp. 907-19.\cr
##' Vannman, K (1995) \emph{A Unified Approach to Capability Indices}. Statitica Sinica,5,805-820.\cr
##' Vannman, K. (2001). \emph{A Graphical Method to Control Process Capability}. Frontiers in Statistical Quality Control, 
##' No 6, Editors: H-J Lenz and P-TH Wilrich. Physica-Verlag, Heidelberg, 290-311.\cr
##' Hubele and Vannman (2004). \emph{The E???ect of Pooled and Un-pooled Variance Estimators on Cpm When Using Subsamples}.
##' Journal Quality Technology, 36, 207-222.\cr
##' @examples
##' library(qcr)
##' data(pistonrings) 
##' xbar <- qcs.xbar(pistonrings[1:125,],plot = TRUE)
##' LSL=73.99; USL=74.01
##' limits = c(lsl = 73.99, usl = 74.01)
##' qcs.ca(xbar, limits = limits)
qcs.ca <- function (object, 
                      limits = c(lsl = -3, usl = 3),
                      target = NULL, std.dev = NULL,nsigmas = 3,
                      confidence = 0.9973, plot = TRUE, main = NULL, ...)
 {

   if (missing(object)){
     if (!inherits(object, "qcs"))
       stop("an object of class 'qcs' is required")
     if (!(object$type == "xbar" | object$type == "one"))
       stop("Process Capability Analysis only available for charts type
            \"qcs.xbar\" and \"qcs.one\" charts")
   }

   x<-object$statistics
   n = sum(object$sizes)
   center = object$center
   std.dev = object$std.dev
   data.name = object$data.name
   parList = list(...)

   if (length(limits)!=2)
     stop("specification limits must be two")
   lsl <- limits[1]
   usl <- limits[2]
   if (lsl>= usl)
     stop("lsl >= usl")
   if (!(is.numeric(usl) & is.finite(lsl)))
     lsl <- NA
   if (!(is.numeric(usl) & is.finite(lsl)))
    usl <- NA
   if (is.na(lsl) & is.na(usl))
     stop("invalid specification limits")
   
   if (is.null(target)) target <- mean(limits, na.rm = TRUE)
   
   if (is.na(lsl)) {
       if (target > usl)
         warning("target value larger than one-sided specification limit...")
     }
   
   if (is.na(usl)) {
       if (target < lsl)
         warning("target value smaller than one-sided specification limit...")
     }
   
   if (!is.na(lsl) & !is.na(usl)) {
       if (target < lsl || target > usl)
         warning("target value is not within specification limits...")
     }
   
   m <- (lsl+usl)/2
   
   if (confidence < 0 | confidence > 1)
     stop("conf.level must be a value between 0 and 1")

   alpha <- 1 - confidence
   
   if (!is.null(lsl) && !is.null(usl)){
     cp <- (usl - lsl)/(2 * nsigmas * std.dev)
     cp.limits <- cp * sqrt(qchisq(c(alpha/2, 1 - alpha/2), n -
                                    1)/(n - 1))
     
     cpm <- cp/sqrt(1 + ((center - target)/std.dev)^2)
     
     df <- n * (1 + ((center - target)/std.dev)^2)/(1 + 2 * ((center -
                                                                target)/std.dev)^2)
     cpm.limits <- cpm * sqrt(qchisq(c(alpha/2, 1 - alpha/2),
                                     df)/df)
     
     cpmk <- qcs.cp(object,parameters = c(1,1),limits = c(lsl,usl),
                    contour = FALSE)[1]
     
     cnp <- qcs.cpn(object,parameters = c(0,0),limits = c(lsl,usl))[1]
     cnpk <- qcs.cpn(object,parameters = c(1,0),limits = c(lsl,usl))[1]
     cnpm <- qcs.cpn(object,parameters = c(0,1),limits = c(lsl,usl))[1]
     cnpmk <- qcs.cpn(object,parameters = c(1,1),limits = c(lsl,usl))[1]
     
     
     
   }else{
     cp <- NA
     cpm <- NA
     cpmk <- NA
     cnp <- NA
     cnpk <- NA
     cnpm <- NA
     cnpmk <- NA
    }
   
   if (!is.null(usl)) {
     cpu <- (usl - center)/(nsigmas * std.dev)
     cpu.limits <- cpu * (1 + c(-1, 1) * qnorm(confidence) *
                              sqrt(1/(9 * n * cpu^2) + 1/(2 * (n - 1))))
   }
   
   if (!is.null(lsl)) {
     cpl <- (center - lsl)/(nsigmas * std.dev)
     cpl.limits <- cpl * (1 + c(-1, 1) * qnorm(confidence) *
                              sqrt(1/(9 * n * cpl^2) + 1/(2 * (n - 1))))
   }
   
   cpk = min(cpu, cpl)
   cpk.limits <- cpk * (1 + c(-1, 1) * qnorm(1 - alpha/2) *
                            sqrt(1/(9 * n * cpk^2) + 1/(2 * (n - 1))))
   


   
   names(cp.limits) <- names(cpk.limits) <- names(cpm.limits) <- 
     c(paste(round(100 *alpha/2, 1), "%", sep = ""), 
       paste(round(100 * (1 - alpha/2), 1), "%", sep = ""))
   
   if (is.na(lsl))
     exp.LSL <- NA
   else {
     exp.LSL <- pnorm((lsl - center)/std.dev) * 100
     if (exp.LSL < 0.01)
       exp.LSL <- 0
   }
   
   if (is.na(usl))
     exp.USL <- NA
   else {
     exp.USL <- (1 - pnorm((usl - center)/std.dev)) * 100
     if (exp.USL < 0.01)
       exp.USL <- 0
   }
   
   obs.LSL <- sum(x < lsl)/n * 100
   obs.USL <- sum(x > usl)/n * 100

   tab <- cbind(c(cp, cpl, cpu, cpk, cpm), rbind(cp.limits,
                                                    cpl.limits, cpu.limits, cpk.limits, cpm.limits))
   rownames(tab) <- c("Cp", "Cp_l", "Cp_u", "Cp_k", "Cpm")
   colnames(tab) <- c("Value", names(cp.limits))
   
   tabn <- cbind(c(cnp, cnpk, cnpm,cnpmk))
   rownames(tabn) <- c("CNp", "CNpK", "CNpm", "CNpmk")
   colnames(tabn) <- c("Value")
   
   if (plot == TRUE) {

     par.orig <- par(c("mar", "oma", "mfrow"))

     if (is.null(parList[["col"]]))
       parList$col = "lightblue"
     if (is.null(parList[["border"]]))
       parList$border = 1
     if (is.null(parList[["lwd"]]))
       parList$lwd = 1
     if (is.null(parList[["cex.axis"]]))
       parList$cex.axis = 1.5

     lineWidth = 1
     lineCol = "red"
     lineType = "solid"
     specCol = "red3"
     specWidth = 1
     cex.text = 1.5
     cex.val = 1
     cex.col = "darkgray"
     bounds.lty = 3
     bounds.col = "red"
     
     xlim <- range(x, usl,lsl, target, na.rm = TRUE)
     xlim <- xlim + diff(xlim) * c(-0.1, 0.1)
     xx <- seq(min(xlim), max(xlim), length = 250)
     dx <- dnorm(xx, center, std.dev)
     if (class(object)[1]=="qcs.xbar"){
        std.dev2 <- sd(object$qcd[,1])
        dx2<- dnorm(xx, center, std.dev2) 
     }
     h <- hist(x[, 1], plot = FALSE)
     ylim <- range(h$density, dx)
     ylim <- ylim + diff(ylim) * c(0, 0.05)
     xlim <- range(x[, 1], usl, lsl)
     xlim <- xlim + diff(xlim) * c(-0.2, 0.2)


     par(bg="#CCCCCC",mar = c(0, 0, 0, 0) + 0.1)

     par(oma = c(2, 4, 7, 4) + 0.1)
     layout(matrix(c(1, 1, 2, 3, 1, 1, 4, 5, 1, 1, 6, 7), 
                   nrow = 3, byrow = TRUE))

          do.call(hist, c(list(x[, 1], freq = FALSE, xlim = xlim,
                          ylim = ylim, main = "",col = "#EEEEEE")))

     tempList = parList
     tempList$col = "#EEEEEE"
     tempList$border = NULL

     abline(v = 0)
     abline(h = 0)
     lines(xx, dx, lwd =1 , col = "red", lty = 1)
     if (class(object)[1]=="qcs.xbar") 
       lines(xx, dx2, lwd = lineWidth, col = "blue", lty = lineType)
     abline(v = usl, col = specCol, lwd = specWidth, lty = 5)
     abline(v = lsl, col = specCol, lwd = specWidth, lty = 5)
     abline(v = target, col = specCol, lwd = specWidth, lty = 5)

     do.call(grid,list(col = "#EEEEEE"))
     do.call(box,tempList)

     if (!is.null(lsl))
       axis(side = 3, at = lsl, labels = paste("LSL =",
                                               format(lsl, digits = 3)), col = specCol)
     if (!is.null(usl))
       axis(side = 3, at = usl, labels = paste("USL =",
                                               format(usl, digits = 3)), col = specCol)
     if (!is.null(lsl) && !is.null(usl))
       axis(side = 3, at = c(lsl, usl), labels = c(paste("LSL =",
                                                         format(lsl, digits = 3)), paste("USL =", format(usl,
                                                                                                         digits = 3))), col = specCol)
     if (!is.null(target))
       text(target, max(ylim), "TARGET", pos = 1, col = "black",
            cex = cex.text)

     if (is.null(main))
         main = paste("Process Capability for",data.name)
     
     if (class(object)[1]=="qcs.xbar") 
       legend("topright",legend = c("ST","LT"), 
              col=c(lineCol,"blue"),lty = c(1,2),bty="n")

#################
     par(mar = c(0, 0, 0, 0) + 0.1)
     title(main = main, outer = TRUE)
     plot(0:5, 0:5, type = "n", axes = FALSE, xlab = "", ylab = "",
          main = "")
     
     rect(par("usr")[1],
          par("usr")[3],
          par("usr")[2],
          par("usr")[4],
          col  =  "white")
     
     
     dis="normal"
     adTestStats = .myADTest(x, dis)
     A = adTestStats$statistic
     p = adTestStats$p.value

     text(3.2, 5, "Process Data", pos = 2,
          cex = cex.val, col = "blue")
     
     text(4, 4.5, paste("Sample n =",n), pos = 2,
          cex = cex.val)
     
     text(4, 4, paste("A =",format(A, digits = 3)),
          pos = 2, cex = cex.val)
     
     if (!is.null(adTestStats$smaller) && adTestStats$smaller)
       text(4, 3.5, paste("p <", format(p, digits = 3)),
            pos = 2, cex = cex.val)
     if (!is.null(adTestStats$smaller) && !adTestStats$smaller)
       text(4, 3.5, paste("p >=", format(p,
                                         digits = 3)), pos = 2, cex = cex.val)
     if (is.null(adTestStats$smaller))
       text(4, 3.5, paste("p =", format(p, digits = 3)),
            pos = 2, cex = cex.val)
     
     text(4, 3, paste("LSL =",format(lsl, digits = 3)), 
          pos = 2,cex = cex.val)

     text(4, 2.5, paste("Target =",format(target, digits = 3)), 
          pos = 2,cex = cex.val)
     
     text(4, 2, paste("USL =",format(usl, digits = 3)), 
          pos = 2,cex = cex.val)
     
     text(4, 1.5, paste("std.dev (ST) =",format(std.dev, digits = 3)), pos = 2,
          cex = cex.val)
     
     text(4, 1, paste("mean =",format(center, digits = 3)), 
          pos = 2,cex = cex.val)
     
    if (class(object)[1]=="qcs.xbar") 
     text(4, 0.5, paste("std.dev (LT) =",format(std.dev2, digits = 3)), pos = 2,
          cex = cex.val)

     
     

     qqPlot(x[, 1], y = dis, ylab = "", main = "",
            axes = FALSE, bounds.lty = bounds.lty,
            bounds.col = bounds.col)
     box()
     do.call(grid,list(col = "#EEEEEE"))
     
 #################
     par(mar = c(0, 0, 0, 0) + 0.1)
     title(main = main, outer = TRUE)
     plot(0:5, 0:5, type = "n", axes = FALSE, xlab = "", ylab = "",
          main = "")
     
     rect(par("usr")[1],
          par("usr")[3],
          par("usr")[2],
          par("usr")[4],
          col  =  "white")
     
     text(4.7, 5, " Parametric Capability Process (ST)", pos = 2,
          cex = cex.val, col = "blue")
 
     
     text(4, 4.2, paste("Cp =",format(cp, digits = 3)), pos = 2,
          cex = cex.val)
     
     text(4, 3.4, paste("Cpu =",format(cpu, digits = 3)),
          pos = 2, cex = cex.val)
     
     text(4, 2.6, paste("Cpl =",format(cpl, digits = 3)),
          pos = 2, cex = cex.val)
     
     text(4, 1.8, paste("Cpk =",format(cpk, digits = 3)), 
          pos = 2,cex = cex.val)
     
     text(4, 1, paste("Cpm =",format(cpm, digits = 3)), 
          pos = 2,cex = cex.val)
     
     text(4, 0.2, paste("Cpmk =",format(cpmk, digits = 3)), 
          pos = 2,cex = cex.val)
     
     plot(0:5, 0:5, type = "n", axes = FALSE, xlab = "", ylab = "",
          main = "")
     rect(par("usr")[1],
          par("usr")[3],
          par("usr")[2],
          par("usr")[4],
          col  =  "white")
     
     text(5, 5, " NonParametric Capability Process (ST)", pos = 2,
          cex = cex.val, col = "blue")

     
     text(4, 4, paste("CNp =",format(cnp, digits = 3)), pos = 2,
          cex = cex.val)
     
     text(4, 3, paste("CNpk =",format(cnpk, digits = 3)),
          pos = 2, cex = cex.val)
     
     text(4, 2, paste("CNpm =",format(cnpm, digits = 3)), 
          pos = 2,cex = cex.val)
     
     text(4, 1, paste("CNpmk =",format(cnpmk, digits = 3)), 
          pos = 2,cex = cex.val)
     
     
     
#################     
     
     par(mar = c(0, 0, 0, 0) + 0.1)
     title(main = main, outer = TRUE)
     plot(0:5, 0:5, type = "n", axes = FALSE, xlab = "", ylab = "",
          main = "")
     
     rect(par("usr")[1],
          par("usr")[3],
          par("usr")[2],
          par("usr")[4],
          col  =  "white")
     
     text(4, 5, " Performance - ST [%]", pos = 2,
          cex = cex.val, col = "blue")
     
     
     text(4, 4.2, paste("Exp<LSL =",format(exp.LSL, digits = 3)), pos = 2,
          cex = cex.val)
     
     text(4, 3.4, paste("Exp>USL =",format(exp.USL, digits = 3)),
          pos = 2, cex = cex.val)
     
     text(4, 2.6, paste("Exp Total =",format(exp.USL+exp.LSL, digits = 3)),
          pos = 2, cex = cex.val)
     
     text(4, 1.8, paste("Obs<LSL =",format(obs.LSL, digits = 3)), 
          pos = 2,cex = cex.val)
     
     text(4, 1, paste("Obs>USL =",format(obs.USL, digits = 3)), 
          pos = 2,cex = cex.val)
     
     text(4, 0.2, paste("Obs Total =",format(obs.USL+obs.LSL, digits = 3)), 
          pos = 2,cex = cex.val)
     
     par(mar = c(0, 0, 0, 0) + 0.2)
     hat.cpm <- qcs.hat.cpm(object, limits = limits,contour = FALSE)
     delta.t<- hat.cpm$delta.t
     gamma.t<-hat.cpm$gamma.t
     delta.e<-hat.cpm$delta.e
     gamma.e<-hat.cpm$gamma.e
     point.delta<-hat.cpm$delta
     point.gamma<-hat.cpm$gamma
     ylim<-c(0,point.gamma+0.2)
     xlim<-c(-max(delta.t)-0.2,max(delta.t)+0.2)
     plot(delta.t, gamma.t, type = "n", axes = FALSE, xlab = "", ylab = "",
          main = "",ylim=ylim,xlim=xlim)

     do.call(grid,list(col = "#EEEEEE"))
     abline(v = 0)
     lines(delta.t, gamma.t, col = "black", lwd = 2,lty = 1)
     lines(delta.e, gamma.e, col = "blue", lwd = 2,lty = 2)
     points(x = point.delta,y = point.gamma,col="red",pch = 21, bg = "red", 
            lwd=2)
     legend("topleft", c("Theory", "Empirical"), 
            lwd = c(2,2), bty="n",lty = c(1,2), col =c("black","blue"), cex =0.7)
     box()
     
     
     
     
     
#################     
     on.exit(par(par.orig))
   }

   digits = 4  
   cat("\nProcess Capability Analysis\n")
   cat("\nCall:\n", deparse(match.call()), "\n\n", sep = "")
   cat(paste(formatC("Number of obs = ", width = 16), formatC(n, 
                                                              width = 12, flag = "-"), formatC("Target = ", width = 10), 
             formatC(signif(target, digits = digits), 
                     flag = "-"), "\n", sep = ""))
   cat(paste(formatC("Center = ", width = 16), formatC(signif(center,digits = digits), width = 12, flag = "-"),
             formatC("LSL = ", width = 10),formatC(signif(lsl,digits = digits), flag = "-")), "\n", sep = "")
   
   cat(paste(formatC("StdDev = ", width = 16), formatC(signif(std.dev,digits = digits), width = 12, flag = "-"), 
             formatC("USL = ", width = 10), formatC(signif(usl, digits = digits), flag = "-")), "\n", sep = "")
   
   cat("\nParemetric Capability indices:\n\n")
   print(tab, digits = 4, na.print = "", print.gap = 2)
   cat("\n")
   
   
   cat("\nNon parametric Capability indices:\n\n")
   print(tabn, digits = 4, na.print = "", print.gap = 2)
   cat("\n")
   
   
   cat("\nPPM:\n\n")
   cat(paste(formatC("Exp<LSL", width = 16), formatC(exp.LSL* 1e+06, width = 12, flag = "-"), 
             formatC("Obs<LSL", width = 10), formatC(obs.LSL* 1e+06, flag = "-")), "\n", sep = "")
   
   cat(paste(formatC("Exp>USL", width = 16), formatC(exp.USL* 1e+06, width = 12, flag = "-"), 
             formatC("Obs>USL", width = 10), formatC(obs.USL* 1e+06, flag = "-")), "\n", sep = "")
   
   cat(paste(formatC("Exp Total", width = 16), formatC((exp.USL+exp.LSL)* 1e+06, width = 12, flag = "-"), 
             formatC("Obs Total", width = 10), formatC((obs.USL+obs.LSL)* 1e+06, flag = "-")), "\n", sep = "")
   
   cat("\nTest:\n\n")
   print(adTestStats)



   result <- list(cp = cp, cpk = cpk, cpl = cpl, cpu = cpu, cpm = cpm, cpmk = cpmk, 
                  cnp = cnp, cnpk = cnpk,cnpm = cnpm, cnpmk = cnpmk,
                  exp.LSL = exp.LSL* 1e+06, exp.USL = exp.USL* 1e+06, exp.T = (exp.LSL+exp.USL)* 1e+06, 
                  usl = usl, lsl = lsl, target = target,
                  obs.LSL = obs.LSL* 1e+06, obs.USL = obs.USL* 1e+06, obs.T = (obs.LSL+obs.USL)* 1e+06)

   invisible(result)
 }