R/qcc.R

Defines functions qcc print.qcc summary.qcc plot.qcc .qcc.c4 stats.xbar sd.xbar limits.xbar stats.S sd.S limits.S stats.R sd.R limits.R stats.xbar.one sd.xbar.one limits.xbar.one stats.p sd.p limits.p stats.np sd.np limits.np stats.c sd.c limits.c stats.u sd.u limits.u shewhart.rules beyond.limits violating.runs oc.curves oc.curves.xbar oc.curves.p oc.curves.c oc.curves.R oc.curves.S qcc.groups qcc.overdispersion.test

Documented in beyond.limits limits.c limits.np limits.p limits.R limits.S limits.u limits.xbar limits.xbar.one oc.curves oc.curves.c oc.curves.p oc.curves.R oc.curves.S oc.curves.xbar plot.qcc print.qcc qcc qcc.groups qcc.overdispersion.test sd.c sd.np sd.p sd.R sd.S sd.u sd.xbar sd.xbar.one shewhart.rules stats.c stats.np stats.p stats.R stats.S stats.u stats.xbar stats.xbar.one summary.qcc violating.runs

#-----------------------------------------------------------------------------#
#                                                                             #
#                     QUALITY CONTROL CHARTS IN R                             #
#                                                                             #
#  An R package for statistical in-line quality control.                      #
#                                                                             #
#  Written by: Luca Scrucca                                                   #
#              Department of Statistics                                       #
#              University of Perugia, ITALY                                   #
#              luca@stat.unipg.it                                             #
#                                                                             #
#-----------------------------------------------------------------------------#

#
#  Main function to create a 'qcc' object
#

qcc <- function(data, type = c("xbar", "R", "S", "xbar.one", "p", "np", "c", "u", "g"), sizes, center, std.dev, limits, data.name, labels, newdata, newsizes, newdata.name, newlabels, nsigmas = 3, confidence.level, rules = shewhart.rules, plot = TRUE, ...)
{
  call <- match.call()
  
  if (missing(data))
     stop("'data' argument is not specified")
  
  if(identical(type, eval(formals(qcc)$type)))
    { type <- as.character(type)[1]
      warning("chart 'type' not specified, assuming \"", type, "\"",
              immediate. = TRUE) }
  if(!exists(paste("stats.", type, sep = ""), mode="function") |
     !exists(paste("sd.", type, sep = ""), mode="function") |
     !exists(paste("limits.", type, sep = ""), mode="function"))
    stop(paste("invalid", type, "control chart. See help(qcc) "))

  if (missing(data.name)) 
     data.name <- deparse(substitute(data))
  data <- data.matrix(data)
  if (missing(sizes)) 
     { if (any(type==c("p", "np", "u")))
          stop(paste("sample 'sizes' must be given for a", type, "Chart"))
       else
          sizes <- apply(data, 1, function(x) sum(!is.na(x)))  }
  else
     { if (length(sizes)==1)
          sizes <- rep(sizes, nrow(data))
       else if (length(sizes) != nrow(data))
                stop("sizes length doesn't match with data") }

  if (missing(labels))
     { if (is.null(rownames(data))) labels <- 1:nrow(data)
       else                         labels <- rownames(data) }

  stats <- paste("stats.", type, sep = "")
  if (!exists(stats, mode="function"))
     stop(paste("function", stats, "is not defined"))
  stats <- do.call(stats, list(data, sizes))
  statistics <- stats$statistics
  if (missing(center)) center <- stats$center

  sd <- paste("sd.", type, sep = "")
  if (!exists(sd, mode="function"))
     stop(paste("function", sd, "is not defined!"))
  missing.std.dev <- missing(std.dev)
  if (missing.std.dev)
     { std.dev <- NULL
       std.dev <- switch(type, 
                         "xbar" = { if(any(sizes > 25)) "RMSDF"
                                    else                "UWAVE-R" },
                         "xbar.one" = "MR",
                         "R" = "UWAVE-R",
                         "S" = "UWAVE-SD",
                         NULL)
       std.dev <- do.call(sd, list(data, sizes, std.dev)) }
  else 
     { if (is.character(std.dev))
          { std.dev <- do.call(sd, list(data, sizes, std.dev)) }
       else
          { if (!is.numeric(std.dev))
               stop("if provided the argument 'std.dev' must be a method available or a numerical value. See help(qcc).")  }
     }

  names(statistics) <-  rownames(data) <-  labels
  names(dimnames(data)) <- list("Group", "Samples")

  object <- list(call = call, type = type, 
                 data.name = data.name, data = data, 
                 statistics = statistics, sizes = sizes, 
                 center = center, std.dev = std.dev)
  # check for new data provided and update object
  if (!missing(newdata))
     {   if (missing(newdata.name))
			{newdata.name <- deparse(substitute(newdata))}
       newdata <- data.matrix(newdata)
       if (missing(newsizes))
          { if (any(type==c("p", "np", "u")))
               stop(paste("sample sizes must be given for a", type, "Chart"))
            else
               newsizes <- apply(newdata, 1, function(x) sum(!is.na(x))) }
       else
          { if (length(newsizes)==1)
               newsizes <- rep(newsizes, nrow(newdata))
            else if (length(newsizes) != nrow(newdata))
                     stop("newsizes length doesn't match with newdata") }
       stats <- paste("stats.", type, sep = "")
       if (!exists(stats, mode="function"))
          stop(paste("function", stats, "is not defined"))
       newstats <- do.call(stats, list(newdata, newsizes))$statistics
       if (missing(newlabels))
          { if (is.null(rownames(newdata)))
               { start <- length(statistics)
                 newlabels <- seq(start+1, start+length(newstats)) }
            else
               { newlabels <- rownames(newdata) }
          }
       names(newstats) <- newlabels
       object$newstats <- newstats
       object$newdata  <- newdata
       object$newsizes <- newsizes
       object$newdata.name <- newdata.name
       statistics <- c(statistics, newstats)
       sizes <- c(sizes, newsizes)
     }

  conf <- nsigmas
  if (!missing(confidence.level))
     conf <- confidence.level
  if (conf >= 1)
     { object$nsigmas <- conf }
  else
     if (conf > 0 & conf < 1)
        { object$confidence.level <- conf } 

  # get control limits
  if (missing(limits))
     { limits <- paste("limits.", type, sep = "")
       if (!exists(limits, mode="function"))
          stop(paste("function", limits, "is not defined"))
       limits <- do.call(limits, list(center = center, std.dev = std.dev,
                                      sizes = sizes, conf = conf)) 
     }
  else 
     { if (!missing.std.dev)
          warning("'std.dev' is not used when limits is given")
       if (!is.numeric(limits))
          stop("'limits' must be a vector of length 2 or a 2-columns matrix")
       limits <- matrix(limits, ncol = 2)
       dimnames(limits) <- list(rep("",nrow(limits)), c("LCL ", "UCL"))
     }

  lcl <- limits[,1]
  ucl <- limits[,2]
  object$limits <- limits
  if (is.function(rules)) violations <- rules(object)
  else                    violations <- NULL
  object$violations <- violations

  class(object) <- "qcc"
  if(plot) plot(object, ...) 
  return(object)
}

print.qcc <- function(x, ...) str(x,1)

summary.qcc <- function(object, digits =  getOption("digits"), ...)
{
  #object <- x   # Argh.  Really want to use 'object' anyway
  cat("\nCall:\n",deparse(object$call),"\n\n",sep="")
  data.name <- object$data.name
  type <- object$type
  cat(paste(type, "chart for", data.name, "\n"))
  statistics <- object$statistics
  cat("\nSummary of group statistics:\n")
  print(summary(statistics), digits = digits, ...)
  sizes <- object$sizes
  if(length(unique(sizes))==1)
     sizes <- sizes[1]
  if(length(sizes) == 1)
     cat("\nGroup sample size: ", format(sizes))
  else {
         cat("\nSummary of group sample sizes: ")
         tab <- table(sizes)
         print(matrix(c(as.numeric(names(tab)), tab), 
                      ncol = length(tab), byrow = TRUE, 
                      dimnames = list(c("  sizes", "  counts"),
                      character(length(tab)))), 
               digits = digits, ...)
        }
  cat("\nNumber of groups: ", length(statistics))
  
  center <- object$center
  if(length(center) == 1)
    { cat("\nCenter of group statistics: ", format(center, digits = digits)) }
  else
    { out <- paste(format(center, digits = digits))
      out <- out[which(cumsum(nchar(out)+1) < getOption("width")-40)]      
      out <- paste0(paste(out, collapse = " "), " ...")
      cat("\nCenter of group statistics: ", out, sep = "")
    }
  
  sd <- object$std.dev
  if(length(sd) == 1)
    { cat("\nStandard deviation: ", format(sd, digits = digits), "\n") }
  else
    { out <- paste(format(sd, digits = digits))
      out <- out[which(cumsum(nchar(out)+1) < getOption("width")-40)]
      out <- paste0(paste(out, collapse = " "), " ...")
      cat("\nStandard deviation: ", out, "\n", sep = "")
    }

  newdata.name <- object$newdata.name
  newstats <- object$newstats
  if (!is.null(newstats)) 
     { cat(paste("\nSummary of group statistics in ", 
                 newdata.name, ":\n", sep = ""))
       print(summary(newstats), digits = digits, ...)
       newsizes <- object$newsizes
       if (length(unique(newsizes)) == 1)
          newsizes <- newsizes[1]
       if (length(newsizes) == 1)
             cat("\nGroup sample size: ", format(newsizes))
       else 
           { cat("\nSummary of group sample sizes:\n")
             new.tab <- table(newsizes)
             print(matrix(c(as.numeric(names(new.tab)), new.tab),
                          ncol = length(new.tab), byrow = TRUE, 
                          dimnames = list(c("  sizes", "  counts"),
                                          character(length(new.tab)))), 
                   digits = digits, ...)
           }
       cat("\nNumber of groups: ", length(newstats), "\n")
     }

  limits <- object$limits
  if (!is.null(limits)) 
     { cat("\nControl limits:\n")
       .printShortMatrix(limits, digits = digits, ...) }

  invisible()
}


plot.qcc <- function(x, add.stats = TRUE, chart.all = TRUE, 
                     label.limits = c("LCL ", "UCL"),
                     title, xlab, ylab, ylim, axes.las = 0,
                     digits =  getOption("digits"),
                     restore.par = TRUE, ...) 
{
  object <- x  # Argh.  Really want to use 'object' anyway
  if ((missing(object)) | (!inherits(object, "qcc")))
     stop("an object of class `qcc' is required")

  # collect info from object
  type <- object$type
  std.dev <- object$std.dev
  data.name <- object$data.name
  center <- object$center
  stats <- object$statistics
  limits <- object$limits 
  lcl <- limits[,1]
  ucl <- limits[,2]
  newstats <- object$newstats
  newdata.name <- object$newdata.name
  violations <- object$violations
  if(chart.all) 
    { statistics <- c(stats, newstats)
      indices <- 1:length(statistics) }
  else
     { if(is.null(newstats))
         { statistics <- stats
           indices <- 1:length(statistics) }
       else
         { statistics <- newstats 
           indices <- seq(length(stats)+1, length(stats)+length(newstats)) }
     }
   
  if (missing(title))
     { if (is.null(newstats))
            main.title <- paste(type, "Chart\nfor", data.name)
       else if (chart.all)
                 main.title <- paste(type, "Chart\nfor", data.name, 
                                     "and", newdata.name)
            else main.title <- paste(type, "Chart\nfor", newdata.name) 
     }
  else main.title <- paste(title)

  oldpar <- par(no.readonly = TRUE)
  if(restore.par) on.exit(par(oldpar))
  mar <- pmax(oldpar$mar, c(4.1,4.1,3.1,2.1))
  par(bg  = qcc.options("bg.margin"), 
      cex = oldpar$cex * qcc.options("cex"),
      mar = if(add.stats) pmax(mar, c(7.6,0,0,0)) else mar)

  # plot Shewhart chart
  plot(indices, statistics, type="n",
       ylim = if(!missing(ylim)) ylim 
              else range(statistics, limits, center),
       ylab = if(missing(ylab)) "Group summary statistics" else ylab,
       xlab = if(missing(xlab)) "Group" else xlab, 
       axes = FALSE)
  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], 
       col = qcc.options("bg.figure"))
  axis(1, at = indices, las = axes.las,
       labels = if(is.null(names(statistics))) 
                   as.character(indices) else names(statistics))
  axis(2, las = axes.las)
  box()
  top.line <- par("mar")[3]-length(capture.output(cat(main.title)))
  top.line <- top.line - if(chart.all & (!is.null(newstats))) 0.1 else 0.5
  mtext(main.title, side = 3, line = top.line,
        font = par("font.main"), 
        cex  = qcc.options("cex"), 
        col  = par("col.main"))

  lines(indices, statistics, type = "b", pch=20) 

  if(length(center) == 1)
       abline(h = center)
  else lines(indices, center[indices], type="s")

  if(length(lcl) == 1) 
    { abline(h = lcl, lty = 2)
      abline(h = ucl, lty = 2) }
  else 
    { lines(indices, lcl[indices], type="s", lty = 2)
      lines(indices, ucl[indices], type="s", lty = 2) }
  mtext(label.limits, side = 4, at = c(rev(lcl)[1], rev(ucl)[1]), 
        las = 1, line = 0.1, col = gray(0.3), cex = par("cex"))
  mtext("CL", side = 4, at = rev(center)[1], 
        las = 1, line = 0.1, col = gray(0.3), cex = par("cex"))

  if(is.null(qcc.options("violating.runs")))
     stop(".qcc.options$violating.runs undefined. See help(qcc.options).")
  if(length(violations$violating.runs))
    { v <- violations$violating.runs
      if(!chart.all & !is.null(newstats))
        { v <- v - length(stats) 
          v <- v[v>0] }
      points(indices[v], statistics[v], 
             col = qcc.options("violating.runs")$col, 
             pch = qcc.options("violating.runs")$pch) 
    }

  if(is.null(qcc.options("beyond.limits")))
     stop(".qcc.options$beyond.limits undefined. See help(qcc.options).")
  if(length(violations$beyond.limits))
    { v <- violations$beyond.limits
      if(!chart.all & !is.null(newstats))
        { v <- v - length(stats) 
          v <- v[v>0] }
      points(indices[v], statistics[v], 
             col = qcc.options("beyond.limits")$col, 
             pch = qcc.options("beyond.limits")$pch) 
    }

  if(chart.all & (!is.null(newstats)))
    { len.obj.stats <- length(object$statistics)
      len.new.stats <- length(statistics) - len.obj.stats
      abline(v = len.obj.stats + 0.5, lty = 3)
      mtext(# paste("Calibration data in", data.name),
            "Calibration data", cex = par("cex")*0.8,
            at = len.obj.stats/2, line = 0, adj = 0.5)
      mtext(# paste("New data in", object$newdata.name),  
            "New data", cex = par("cex")*0.8, 
            at = len.obj.stats + len.new.stats/2, line = 0, adj = 0.5)
    }

  if(add.stats) 
    { 
      # computes the x margins of the figure region
      plt <- par()$plt; usr <- par()$usr
      px <- diff(usr[1:2])/diff(plt[1:2])
      xfig <- c(usr[1]-px*plt[1], usr[2]+px*(1-plt[2]))
      at.col <- xfig[1] + diff(xfig[1:2])*c(0.10, 0.40, 0.65)
      top.line <- 4.5
      # write info at bottom
      mtext(paste("Number of groups = ", length(statistics), sep = ""), 
            side = 1, line = top.line, adj = 0, at = at.col[1],
            font = qcc.options("font.stats"),
            cex = par("cex")*qcc.options("cex.stats"))
      center <- object$center
      if(length(center) == 1)
        { mtext(paste("Center = ", signif(center[1], digits), sep = ""),
                side = 1, line = top.line+1, adj = 0, at = at.col[1],
                font = qcc.options("font.stats"),
                cex = par("cex")*qcc.options("cex.stats"))
        }
      else 
        { mtext("Center is variable",
                side = 1, line = top.line+1, adj = 0, at = at.col[1],
                font = qcc.options("font.stats"),
                cex = par("cex")*qcc.options("cex.stats"))
        }
      
      if(length(std.dev) == 1)
        { mtext(paste("StdDev = ", signif(std.dev, digits), sep = ""),
                side = 1, line = top.line+2, adj = 0, at = at.col[1],
                font = qcc.options("font.stats"),
                cex = par("cex")*qcc.options("cex.stats"))
        }      
      else
        { mtext("StdDev is variable",
                side = 1, line = top.line+2, adj = 0, at = at.col[1],
                font = qcc.options("font.stats"),
                cex = par("cex")*qcc.options("cex.stats"))
        }
      
      if(length(unique(lcl)) == 1)
        { mtext(paste("LCL = ", signif(lcl[1], digits), sep = ""), 
                side = 1, line = top.line+1, adj = 0, at = at.col[2],
                font = qcc.options("font.stats"),
                cex = par("cex")*qcc.options("cex.stats"))
        }
      else 
        { mtext("LCL is variable", 
                side = 1, line = top.line+1, adj = 0, at = at.col[2],
                font = qcc.options("font.stats"),
                cex = par("cex")*qcc.options("cex.stats"))
        }
      
      if(length(unique(ucl)) == 1)
         { mtext(paste("UCL = ", signif(ucl[1], digits), sep = ""),
                 side = 1, line = top.line+2, adj = 0, at = at.col[2],
                 font = qcc.options("font.stats"),
                 cex = par("cex")*qcc.options("cex.stats")) 
         }
       else 
         { mtext("UCL is variable", 
                 side = 1, line = top.line+2, adj = 0, at = at.col[2],
                 font = qcc.options("font.stats"),
                 cex = par("cex")*qcc.options("cex.stats"))
         }
      
       if(!is.null(violations))
         { mtext(paste("Number beyond limits =",
                       length(unique(violations$beyond.limits))), 
                 side = 1, line = top.line+1, adj = 0, at = at.col[3],
                 font = qcc.options("font.stats"),
                 cex = par("cex")*qcc.options("cex.stats"))
           mtext(paste("Number violating runs =",
                       length(unique(violations$violating.runs))), 
                 side = 1, line = top.line+2, adj = 0, at = at.col[3],
                 font = qcc.options("font.stats"),
                 cex = par("cex")*qcc.options("cex.stats"))
         }
     }

  invisible()
}

#
#  Functions used to compute Shewhart charts statistics
#

.qcc.c4 <- function(n)
{ sqrt(2/(n - 1)) * exp(lgamma(n/2) - lgamma((n - 1)/2)) }
        
# xbar

stats.xbar <- function(data, sizes)
{
  if (missing(sizes))
     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
  statistics <- apply(data, 1, mean, na.rm=TRUE)
  center <- sum(sizes * statistics)/sum(sizes)
  list(statistics = statistics, center = center)
}

sd.xbar <- function(data, sizes, std.dev = c("UWAVE-R", "UWAVE-SD", "MVLUE-R", "MVLUE-SD", "RMSDF"))
{
  if (!is.numeric(std.dev))
     std.dev <- match.arg(std.dev)
  if (missing(sizes))
     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
  if (any(sizes == 1))
     stop("group sizes must be larger than one")
  c4 <- .qcc.c4
  if (is.numeric(std.dev))
     { sd <- std.dev }
  else
     { switch(std.dev, 
              "UWAVE-R" = {  R <- apply(data, 1, function(x) 
                                                 diff(range(x, na.rm = TRUE)))
                             d2 <- qcc.options("exp.R.unscaled")[sizes]
                             sd <- sum(R/d2)/length(sizes) }, 
              "UWAVE-SD" = { S <- apply(data, 1, sd, na.rm = TRUE)
                             sd <- sum(S/c4(sizes))/length(sizes) },
              "MVLUE-R"  = { R <- apply(data, 1, function(x) 
                                                 diff(range(x, na.rm = TRUE)))
                             d2 <- qcc.options("exp.R.unscaled")[sizes]
                             d3 <- qcc.options("se.R.unscaled")[sizes]
                             w  <- (d2/d3)^2
                             sd <- sum(R/d2*w)/sum(w) }, 
              "MVLUE-SD" = { S <- apply(data, 1, sd, na.rm = TRUE)
                             w  <- c4(sizes)^2/(1-c4(sizes)^2)
                             sd <- sum(S/c4(sizes)*w)/sum(w) },
              "RMSDF" =    { S <- apply(data, 1, sd, na.rm = TRUE)
                             w  <- sizes-1
                             sd <- sqrt(sum(S^2*w)/sum(w))/c4(sum(w)+1) }
              )
     }
#  if (missing(std.dev))
#     var.within <- apply(data, 1, var, na.rm=TRUE)
#  else 
#     var.within <- std.dev^2
#  var.df <- sum(sizes - 1)
#  if (equal.sd) 
#     { std.dev <- sqrt(sum((sizes - 1) * var.within)/var.df) / c4(var.df + 1) }
#  else 
#     { c <- c4(sizes)/(1 - c4(sizes)^2)
#       std.dev <- sum(c * sqrt(var.within))/sum(c * c4(sizes)) }
  return(sd)
}

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

# S chart

stats.S <- function(data, sizes)
{
  if (missing(sizes))
     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
  if(ncol(data)==1) 
    { statistics <- as.vector(data) }
  else 
    { statistics <- sqrt(apply(data, 1, var, na.rm=TRUE)) }
  if (length(sizes == 1))
     sizes <- rep(sizes, length(statistics))
  center <- sum(sizes * statistics)/sum(sizes)
  list(statistics = statistics, center = center)
}

sd.S <- function(data, sizes, std.dev = c("UWAVE-SD", "MVLUE-SD", "RMSDF"))
{
  if (!is.numeric(std.dev))
     std.dev <- match.arg(std.dev)
  sd.xbar(data, sizes, std.dev)
}

limits.S <- function(center, std.dev, sizes, conf)
{
  if (length(unique(sizes))==1)
     sizes <- sizes[1]
  c4 <- .qcc.c4
  se.stats <- std.dev * sqrt(1 - c4(sizes)^2)
  if (conf >= 1) 
     { lcl <- pmax(0, center - conf * se.stats)
       ucl <- center + conf * se.stats
     }
  else 
     { if (conf > 0 & conf < 1) 
          { ucl <- std.dev * sqrt(qchisq(1 - (1 - conf)/2, sizes - 1)/
                                  (sizes - 1))
            lcl <- std.dev * sqrt(qchisq((1 - conf)/2, sizes - 1)/
                                  (sizes - 1))
          }
          else stop("invalid conf argument. See help.")
     }
  limits <- matrix(c(lcl, ucl), ncol = 2)
  rownames(limits) <- rep("", length = nrow(limits))
  colnames(limits) <- c("LCL", "UCL")
  limits
}

# R Chart 

stats.R <- function(data, sizes)
{
  if (missing(sizes))
     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
  if(ncol(data)==1) 
    { statistics <- as.vector(data) }
  else 
    { statistics <- apply(data, 1, function(x) diff(range(x, na.rm=TRUE))) }
  if (length(sizes == 1))
     sizes <- rep(sizes, length(statistics))
  center <- sum(sizes * statistics)/sum(sizes)
  list(statistics = statistics, center = center)
}

sd.R <- function(data, sizes, std.dev = c("UWAVE-R", "MVLUE-R"))
{
  if (!is.numeric(std.dev))
     std.dev <- match.arg(std.dev)
  sd.xbar(data, sizes, std.dev)
}

limits.R <- function(center, std.dev, sizes, conf)
{
  if (length(unique(sizes))==1)
     sizes <- sizes[1]
  se.R.unscaled <- qcc.options("se.R.unscaled")
  Rtab <- length(se.R.unscaled)
  if (conf >= 1) 
     { if (any(sizes > Rtab))
          stop(paste("group size must be less than", 
                      Rtab + 1, "when giving nsigmas"))
       se.R <- se.R.unscaled[sizes] * std.dev
       lcl <- pmax(0, center - conf * se.R)
       ucl <- center + conf * se.R
     }
  else 
     { if (conf > 0 && conf < 1) 
          { ucl <- qtukey(1 - (1 - conf)/2, sizes, 1e100) * std.dev
            lcl <- qtukey((1 - conf)/2, sizes, 1e100) * std.dev
          }
       else stop("invalid conf argument. See help.")
     }
  limits <- matrix(c(lcl, ucl), ncol = 2)
  rownames(limits) <- rep("", length = nrow(limits))
  colnames(limits) <- c("LCL", "UCL")
  return(limits)
}

# xbar Chart for one-at-time data

stats.xbar.one <- function(data, sizes)
{
  statistics <- as.vector(data)
  center <- mean(statistics)
  list(statistics = statistics, center = center)
}

sd.xbar.one <- function(data, sizes, std.dev = c("MR", "SD"), k = 2)
{
  data <- as.vector(data)
  n <- length(data)
  if(!is.numeric(std.dev)) 
     std.dev <- match.arg(std.dev)
  c4 <- .qcc.c4
  if(is.numeric(std.dev)) 
    { sd <- std.dev }
  else
    { switch(std.dev, 
             "MR" = { d2 <- qcc.options("exp.R.unscaled")
                      if(is.null(d2))
                         stop(".qcc.options$exp.R.unscaled is null")
                      d <- 0
                      for(j in k:n)
                          d <- d+abs(diff(range(data[c(j:(j-k+1))])))
                      sd <- (d/(n-k+1))/d2[k] },
             "SD" = { sd <- sd(data)/c4(n) },
             sd <- NULL)
    }
  return(sd)
}


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


# p Chart

stats.p <- function(data, sizes)
{
  data <- as.vector(data)
  sizes <- as.vector(sizes)
  pbar <- sum(data)/sum(sizes)
  list(statistics = data/sizes, center = pbar)
}

sd.p <- function(data, sizes, ...)
{
  data <- as.vector(data)
  sizes <- as.vector(sizes)
  pbar <- sum(data)/sum(sizes)
  std.dev <- sqrt(pbar * (1 - pbar))
  return(std.dev)
}

limits.p <- function(center, std.dev, sizes, conf)
{ 
  limits.np(center * sizes, std.dev, sizes, conf) / sizes
}

# np Chart

stats.np <- function(data, sizes)
{
  data <- as.vector(data)
  sizes <- as.vector(sizes)
  pbar <- sum(data)/sum(sizes)
  center <- sizes * pbar
  if (length(unique(center)) == 1)
     center <- center[1]
  list(statistics = data, center = center)
}

sd.np <- function(data, sizes, ...)
{
  data <- as.vector(data)
  sizes <- as.vector(sizes)
  pbar <- sum(data)/sum(sizes)
  std.dev <- sqrt(sizes * pbar * (1 - pbar))
  if (length(unique(std.dev)) == 1)
     std.dev <- std.dev[1]
  return(std.dev)
}

limits.np <- function(center, std.dev, sizes, conf)
{ 
  sizes <- as.vector(sizes)
  if (length(unique(sizes)) == 1)
     sizes <- sizes[1]
  pbar <- mean(center / sizes)
  if (conf >= 1)
     { tol <- conf * sqrt(pbar * (1 - pbar) * sizes)
       lcl <- pmax(center - tol, 0)
       ucl <- pmin(center + tol, sizes)
     }
  else
     { if (conf > 0 & conf < 1)
          { lcl <- qbinom((1 - conf)/2, sizes, pbar)
            ucl <- qbinom((1 - conf)/2, sizes, pbar, lower.tail = FALSE)
          }
       else stop("invalid conf argument. See help.")
     }
  limits <- matrix(c(lcl, ucl), ncol = 2)
  rownames(limits) <- rep("", length = nrow(limits))
  colnames(limits) <- c("LCL", "UCL")
  return(limits)
}

# c Chart

stats.c <- function(data, sizes)
{
  data <- as.vector(data)
  sizes <- as.vector(sizes)
  if (length(unique(sizes)) != 1)
     stop("all sizes must be be equal for a c chart")
  statistics <- data
  center <- mean(statistics)
  list(statistics = statistics, center = center)
}

sd.c <- function(data, sizes, ...)
{
  data <- as.vector(data)
  std.dev <- sqrt(mean(data))
  return(std.dev)
}

limits.c <- function(center, std.dev, sizes, conf)
{
  if (conf >= 1) 
     { lcl <- center - conf * sqrt(center)
       lcl[lcl < 0] <- 0
       ucl <- center + conf * sqrt(center)
     }
  else 
     { if (conf > 0 & conf < 1) 
          { ucl <- qpois(1 - (1 - conf)/2, center)
            lcl <- qpois((1 - conf)/2, center)
          }
       else stop("invalid conf argument. See help.")
     }
  limits <- matrix(c(lcl, ucl), ncol = 2)
  rownames(limits) <- rep("", length = nrow(limits))
  colnames(limits) <- c("LCL", "UCL")
  return(limits)
}

# u Chart

stats.u <- function(data, sizes)
{
  data <- as.vector(data)
  sizes <- as.vector(sizes)
  statistics <- data/sizes
  center <- sum(sizes * statistics)/sum(sizes)
  list(statistics = statistics, center = center)
}

sd.u <- function(data, sizes, ...)
{
  data <- as.vector(data)
  sizes <- as.vector(sizes)
  std.dev <- sqrt(sum(data)/sum(sizes))
  return(std.dev)
}

limits.u <- function(center, std.dev, sizes, conf)
{
  sizes <- as.vector(sizes)
  if (length(unique(sizes))==1)
     sizes <- sizes[1]
  limits.c(center * sizes, std.dev, sizes, conf) / sizes
}

#
# Functions used to signal points out of control 
#

shewhart.rules <- function(object, limits = object$limits, run.length = qcc.options("run.length"))
{
# Return a list of cases beyond limits and violating runs
  bl <- beyond.limits(object, limits = limits)
  vr <- violating.runs(object, run.length = run.length)
  list(beyond.limits = bl, violating.runs = vr)
}

beyond.limits <- function(object, limits = object$limits)
{
# Return cases beyond limits
  statistics <- c(object$statistics, object$newstats) 
  lcl <- limits[,1]
  ucl <- limits[,2]
  index.above.ucl <- seq(along = statistics)[statistics > ucl]
  index.below.lcl <- seq(along = statistics)[statistics < lcl]
  return(c(index.above.ucl, index.below.lcl))
}

violating.runs <- function(object, run.length = qcc.options("run.length"))
{
# Return indices of points violating runs
  if(run.length == 0)
    return(numeric())
  center <- object$center
  statistics <- c(object$statistics, object$newstats)
  cl <- object$limits
  diffs <- statistics - center
  diffs[diffs > 0] <- 1
  diffs[diffs < 0] <- -1
  runs <- rle(diffs)
  vruns <- rep(runs$lengths >= run.length, runs$lengths)
  vruns.above <- (vruns & (diffs > 0))
  vruns.below <- (vruns & (diffs < 0))
  rvruns.above <- rle(vruns.above)
  rvruns.below <- rle(vruns.below)
  vbeg.above <- cumsum(rvruns.above$lengths)[rvruns.above$values] -
                      (rvruns.above$lengths - run.length)[rvruns.above$values]
  vend.above <- cumsum(rvruns.above$lengths)[rvruns.above$values]
  vbeg.below <- cumsum(rvruns.below$lengths)[rvruns.below$values] -
                      (rvruns.below$lengths - run.length)[rvruns.below$values]
  vend.below <- cumsum(rvruns.below$lengths)[rvruns.below$values]
  violators <- numeric()
  if (length(vbeg.above)) 
     { for (i in 1:length(vbeg.above))
           violators <- c(violators, vbeg.above[i]:vend.above[i]) }
  if (length(vbeg.below)) 
     { for (i in 1:length(vbeg.below))
           violators <- c(violators, vbeg.below[i]:vend.below[i]) }
  return(violators)
}

#-------------------------------------------------------------------#
#                                                                   #
#          Operating Characteristic Function                        #
#                                                                   #
#-------------------------------------------------------------------#

oc.curves <- function(object, ...)
{
# Draws the operating characteristic curves for the qcc object 

  if ((missing(object)) | (!inherits(object, "qcc")))
     stop("an object of class 'qcc' is required")

  size <- unique(object$sizes)
  if (length(size)>1)
     stop("Operating characteristic curves available only for equal sample sizes!")

  beta <- switch(object$type,
                 xbar = oc.curves.xbar(object, ...),
                 R    = oc.curves.R(object, ...),
                 S    = oc.curves.S(object, ...),
                 np   =,
                 p    = oc.curves.p(object, ...),
                 u    =,
                 c    = oc.curves.c(object, ...))
  if (is.null(beta))
     stop("Operating characteristic curves not available for this type of chart.")

  invisible(beta)
}

oc.curves.xbar <- function(object, n, c = seq(0, 5, length=101), nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
{
# Draw the operating-characteristic curves for the xbar-chart with nsigmas
# limits. The values on the vertical axis give the probability of not detecting
# a shift of c*sigma in the mean on the first sample following the shift.

  if (!(object$type=="xbar"))
     stop("not a `qcc' object of type \"xbar\".")

  size <- unique(object$sizes)
  if (length(size) > 1)
     stop("Operating characteristic curves available only for equal sample sizes!")
  if (missing(n))
     n <- unique(c(size, c(1,5,10,15,20)))
  if (is.null(nsigmas))
     nsigmas <- qnorm(1 - (1 - object$confidence.level) / 2)

  beta <- matrix(as.double(NA), length(n), length(c))
  for (i in 1:length(n))
      beta[i,] <- pnorm(nsigmas-c*sqrt(n[i])) - pnorm(-nsigmas-c*sqrt(n[i]))
  rownames(beta) <- paste("n=",n,sep="")
  colnames(beta) <- c

  oldpar <- par(no.readonly = TRUE)
  if(restore.par) on.exit(par(oldpar))
  par(bg  = qcc.options("bg.margin"), 
      cex = oldpar$cex * qcc.options("cex"),
      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))

  plot(c, beta[1,], type="n",
       ylim = c(0,1), xlim = c(0,max(c)),
       xlab = "Process shift (std.dev)",
       ylab = "Prob. type II error ")
  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], 
       col = qcc.options("bg.figure"))
  box()
  mtext(paste("OC curves for", object$type, "Chart"), 
        side = 3, line = par("mar")[3]/3,
        font = par("font.main"), 
        cex  = qcc.options("cex"), 
        col  = par("col.main"))
  for(i in 1:length(n))
     lines(c, beta[i,], type = "l", lty=i)
  beta <- t(beta)
  names(dimnames(beta)) <- c("shift (std.dev)", "sample size")

  if (identify)
     { cs <- rep(c,length(n))
       betas <- as.vector(beta)
       labels <- paste("c=", formatC(cs, 2, flag="-"), 
                       ": beta=", formatC(betas, 4, flag="-"), 
                       ", ARL=", formatC(1/(1-betas), 2, flag="-"), sep="")
       i <- identify(cs, betas, labels, pos=4, offset=0.2)
       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
     }
  else
     { legend(max(c), 1, legend = paste("n =", n), 
              bg = qcc.options("bg.figure"),
              lty = 1:length(n), xjust = 1, yjust = 1)
     }
  invisible(beta)
}

oc.curves.p <- function(object, nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
{
  if (!(object$type=="p" | object$type=="np"))
     stop("not a `qcc' object of type \"p\" or \"np\".")

  size <- unique(object$sizes)
  if (length(size) > 1)
     stop("Operating characteristic curves available only for equal sample sizes!")

  if (is.null(object$limits))
     stop("the `qcc' object does not have control limits!")
  limits <- object$limits
  p <- seq(0, 1, length=101)

  if (object$type=="p") 
     { UCL <- min(floor(size*limits[,2]), size)
       LCL <- max(floor(size*limits[,1]), 0) }
  else
     { UCL <- min(floor(limits[,2]), size)
       LCL <- max(floor(limits[,1]), 0) }
  beta <- pbinom(UCL, size, p) - pbinom(LCL-1, size, p)
  names(beta) <- p

  oldpar <- par(no.readonly = TRUE)
  if(restore.par) on.exit(par(oldpar))
  par(bg  = qcc.options("bg.margin"), 
      cex = oldpar$cex * qcc.options("cex"),
      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))

  plot(p, beta, type = "n", 
       ylim = c(0,1), xlim = c(0,1),
       xlab = expression(p), 
       ylab = "Prob. type II error ")
  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], 
       col = qcc.options("bg.figure"))
  box()
  mtext(paste("OC curves for", object$type, "Chart"), 
        side = 3, line = par("mar")[3]/3,
        font = par("font.main"), 
        cex  = qcc.options("cex"), 
        col  = par("col.main"))

  lines(p, beta)
  lines(rep(p[which.max(beta)], 2), c(0, max(beta)), lty = 2)

  warning("Some computed values for the type II error have been rounded due to the discreteness of the binomial distribution. Thus, some ARL values might be meaningless.")

  if (identify)
     { labels <- paste("p=", formatC(p, 2, flag="-"), 
                       ": beta=", formatC(beta, 4, flag="-"), 
                       ", ARL=", formatC(1/(1-beta), 2, flag="-"), sep="")
       i <- identify(p, beta, labels, pos=4, offset=0.2)
       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
     }
  invisible(beta)  
}

oc.curves.c <- function(object, nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
{
  type <- object$type
  if (!(object$type=="c" | object$type=="u"))
     stop("not a `qcc' object of type \"c\" or \"u\".")

  size <- unique(object$sizes)
  if (length(size) > 1)
     stop("Operating characteristic curves available only for equal sample size!")

  if (is.null(object$limits))
     stop("the `qcc' object does not have control limits!")
  limits <- object$limits
  CL  <- object$center
  std.dev <- object$std.dev
  if (object$type=="c") 
     { max.lambda <- ceiling(CL+10*std.dev)
       UCL <- floor(limits[1,2])
       LCL <- floor(limits[1,1])
     }
  else
     { max.lambda <- ceiling(CL*size+10*std.dev)[1]
       UCL <- floor(size*limits[1,2])
       LCL <- floor(size*limits[1,1])
     }
  lambda <- seq(0, max.lambda)
  beta <- ppois(UCL, lambda) - ppois(LCL-1, lambda)
  names(beta) <- lambda

  oldpar <- par(no.readonly = TRUE)
  if(restore.par) on.exit(par(oldpar))
  par(bg  = qcc.options("bg.margin"), 
      cex = oldpar$cex * qcc.options("cex"),
      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))

  plot(lambda, beta, type = "n", 
       ylim = c(0,1), xlim = range(lambda),
       xlab = "Mean", 
       ylab = "Prob. type II error ")
  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], 
       col = qcc.options("bg.figure"))
  box()
  mtext(paste("OC curves for", object$type, "Chart"), 
        side = 3, line = par("mar")[3]/3,
        font = par("font.main"), 
        cex  = qcc.options("cex"), 
        col  = par("col.main"))

  lines(lambda, beta)     
  lines(rep(lambda[which.max(beta)], 2), c(0, max(beta)), lty = 2)

  warning("Some computed values for the type II error have been rounded due to the discreteness of the Poisson distribution. Thus, some ARL values might be meaningless.")

  if (identify)
     { labels <- paste("lambda=", formatC(lambda, 0, flag="-"), 
                       ": beta=", formatC(beta, 4, flag="-"), 
                       ", ARL=", formatC(1/(1-beta), 2, flag="-"), sep="")
       i <- identify(lambda, beta, labels, pos=4, offset=0.2)
       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
     }
  invisible(beta)
}

oc.curves.R <-
function(object, n, c = seq(1, 6, length=101), nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
{
# Draw the operating-characteristic curves for the R-chart with nsigmas
# limits. The values on the vertical axis give the probability of not detecting
# a change from sigma to c*sigma on the first sample following the change.

  if (!(object$type=="R"))
     stop("not a `qcc' object of type \"R\".")

  size <- unique(object$sizes)
  if (length(size) > 1)
     stop("Operating characteristic curves available only for equal sample sizes!")
  if (missing(n))
     n <- unique(c(size, c(2,5,10,15,20)))
  if (is.null(nsigmas))
    { tail.prob <- (1 - object$confidence.level) / 2
      beta.fun1 <- function(c, n, p)
      {
        lcl <- qtukey(p, n, Inf)
        ucl <- qtukey(p, n, Inf, lower.tail = FALSE)
        ptukey(ucl / c, n, Inf) - ptukey(lcl / c, n, Inf)
      }
      beta <- outer(c, n, beta.fun1, tail.prob)
    }
  else
    { exp.R.unscaled <- qcc.options("exp.R.unscaled")
      se.R.unscaled <- qcc.options("se.R.unscaled")
      Rtab <- min(length(exp.R.unscaled), length(se.R.unscaled))
      if (any(n > Rtab))
          stop(paste("group size must be less than",
                      Rtab + 1, "when giving nsigmas"))
      beta.fun2 <- function(c, n, conf)
      {
        d2 <- exp.R.unscaled[n]
        d3 <- se.R.unscaled[n]
        lcl <- pmax(0, d2 - conf * d3)
        ucl <- d2 + conf * d3
        ptukey(ucl / c, n, Inf) - ptukey(lcl / c, n, Inf)
      }
      beta <- outer(c, n, beta.fun2, nsigmas)
    }

  colnames(beta) <- paste("n=",n,sep="")
  rownames(beta) <- c

  oldpar <- par(no.readonly = TRUE)
  if(restore.par) on.exit(par(oldpar))
  par(bg  = qcc.options("bg.margin"), 
      cex = oldpar$cex * qcc.options("cex"),
      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))

  plot(c, beta[,1], type="n",
       ylim = c(0,1), xlim = c(1,max(c)),
       xlab = "Process scale multiplier",
       ylab = "Prob. type II error ")
  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
       col = qcc.options("bg.figure"))
  box()
  mtext(paste("OC curves for", object$type, "Chart"), 
        side = 3, line = par("mar")[3]/3,
        font = par("font.main"), 
        cex  = qcc.options("cex"), 
        col  = par("col.main"))
  
  matlines(c, beta, lty = 1:length(n), col = 1)

  names(dimnames(beta)) <- c("scale multiplier", "sample size")

  if (identify)
     { cs <- rep(c,length(n))
       betas <- as.vector(beta)
       labels <- paste("c=", formatC(cs, 2, flag="-"),
                       ": beta=", formatC(betas, 4, flag="-"),
                       ", ARL=", formatC(1/(1-betas), 2, flag="-"), sep="")
       i <- identify(cs, betas, labels, pos=4, offset=0.2)
       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
     }
  else
     { legend(max(c), 1, legend = paste("n =", n),
              bg = qcc.options("bg.figure"),
              lty = 1:length(n), xjust = 1, yjust = 1)
     }
  invisible(beta)
}

oc.curves.S <- function(object, n, c = seq(1, 6, length=101), nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
{
# Draw the operating-characteristic curves for the S-chart with nsigmas
# limits. The values on the vertical axis give the probability of not detecting
# a change from sigma to c*sigma on the first sample following the change.

  if (!(object$type=="S"))
     stop("not a `qcc' object of type \"S\".")

  size <- unique(object$sizes)
  if (length(size) > 1)
     stop("Operating characteristic curves available only for equal sample sizes!")
  if (missing(n))
     n <- unique(c(size, c(2,5,10,15,20)))
  if (is.null(nsigmas))
    { tail.prob <- (1 - object$confidence.level) / 2
      beta.fun1 <- function(c, n, p)
      {
        ucl <- sqrt(qchisq(1 - p, n - 1) / (n - 1))
        lcl <- sqrt(qchisq(p, n - 1) / (n - 1))
        pchisq((n - 1) * (ucl / c)^2, n - 1) - pchisq((n - 1)* (lcl / c)^2, n - 1)
      }
      beta <- outer(c, n, beta.fun1, tail.prob)
    }
  else
    { c4 <- .qcc.c4
      beta.fun2 <- function(c, n)
      {
        center <- c4(n)
        tol <- sqrt(1 - c4(n)^2)
        lcl <- pmax(0, center - nsigmas * tol)
        ucl <- center + nsigmas * tol
        pchisq((n - 1) * (ucl / c)^2, n - 1) - pchisq((n - 1) * (lcl / c)^2, n - 1)
      }
      beta <- outer(c, n, beta.fun2)
    }

  colnames(beta) <- paste("n=",n,sep="")
  rownames(beta) <- c

  oldpar <- par(no.readonly = TRUE)
  if(restore.par) on.exit(par(oldpar))
  par(bg  = qcc.options("bg.margin"), 
      cex = oldpar$cex * qcc.options("cex"),
      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))

  plot(c, beta[,1], type="n",
       ylim = c(0,1), xlim = c(1,max(c)),
       xlab = "Process scale multiplier",
       ylab = "Prob. type II error ")
  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
       col = qcc.options("bg.figure"))
  box()
  mtext(paste("OC curves for", object$type, "Chart"), 
        side = 3, line = par("mar")[3]/3,
        font = par("font.main"), 
        cex  = qcc.options("cex"), 
        col  = par("col.main"))
  matlines(c, beta, lty = 1:length(n), col = 1)

  names(dimnames(beta)) <- c("scale multiplier", "sample size")

  if (identify)
     { cs <- rep(c,length(n))
       betas <- as.vector(beta)
       labels <- paste("c=", formatC(cs, 2, flag="-"),
                       ": beta=", formatC(betas, 4, flag="-"),
                       ", ARL=", formatC(1/(1-betas), 2, flag="-"), sep="")
       i <- identify(cs, betas, labels, pos=4, offset=0.2)
       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
     }
  else
     { legend(max(c), 1, legend = paste("n =", n),
              bg = qcc.options("bg.figure"),
              lty = 1:length(n), xjust = 1, yjust = 1)
     }
  invisible(beta)
}

#-------------------------------------------------------------------#
#                                                                   #
# Miscellaneous functions                                           #
#                                                                   #
#-------------------------------------------------------------------#

qcc.groups <- function(data, sample)
{
  if(length(data)!=length(sample))
    stop("data and sample must be vectors of equal length")
  x <- lapply(split(data, sample), as.vector)
  lx <- sapply(x, length)
  for(i in which(lx != max(lx)))
      x[[i]] <- c(x[[i]], rep(NA, max(lx)-lx[i]))
  x <- t(sapply(x, as.vector))
  return(x)
}

qcc.overdispersion.test <- function(x, size, 
                            type=ifelse(missing(size), "poisson", "binomial"))
{
  type <- match.arg(type, c("poisson", "binomial"))
  if (type=="binomial" & missing(size))
     stop("binomial data require argument \"size\"")
  if (!missing(size))
     if (length(x) != length(size))   
        stop("arguments \"x\" and \"size\" must be vector of same length")

  n <- length(x)
  obs.var <- var(x)
  if (type=="binomial")
     { p <- sum(x)/sum(size)
       theor.var <- mean(size)*p*(1-p) }
  else if (type=="poisson")
          { theor.var <- mean(x) }
       else
          stop("invalid \"type\" argument. See help.")

  D <- (obs.var * (n-1)) / theor.var
  p.value <- 1-pchisq(D, n-1)

  out <- matrix(c(obs.var/theor.var, D, signif(p.value,5)), 1, 3)
  rownames(out) <- paste(type, "data")
  colnames(out) <- c("Obs.Var/Theor.Var", "Statistic", "p-value") 
  names(dimnames(out)) <- c(paste("Overdispersion test"), "")
  return(out)
}

Try the qcc package in your browser

Any scripts or data that you put into this service are public.

qcc documentation built on May 2, 2019, 9:15 a.m.