1#-----------------------------------------------------------------------------#
2#                                                                             #
3#                     QUALITY CONTROL CHARTS IN R                             #
4#                                                                             #
5#  An R package for statistical in-line quality control.                      #
6#                                                                             #
7#  Written by: Luca Scrucca                                                   #
8#              Department of Statistics                                       #
9#              University of Perugia, ITALY                                   #
10#              luca@stat.unipg.it                                             #
11#                                                                             #
12#-----------------------------------------------------------------------------#
13
14#
15#  Main function to create a 'qcc' object
16#
17
18qcc <- 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, ...)
19{
20  call <- match.call()
21
22  if (missing(data))
23     stop("'data' argument is not specified")
24
25  if(identical(type, eval(formals(qcc)$type)))
26    { type <- as.character(type)[1]
27      warning("chart 'type' not specified, assuming \"", type, "\"",
28              immediate. = TRUE) }
29  if(!exists(paste("stats.", type, sep = ""), mode="function") |
30     !exists(paste("sd.", type, sep = ""), mode="function") |
31     !exists(paste("limits.", type, sep = ""), mode="function"))
32    stop(paste("invalid", type, "control chart. See help(qcc) "))
33
34  if (missing(data.name))
35     data.name <- deparse(substitute(data))
36  data <- data.matrix(data)
37  if (missing(sizes))
38     { if (any(type==c("p", "np", "u")))
39          stop(paste("sample 'sizes' must be given for a", type, "Chart"))
40       else
41          sizes <- apply(data, 1, function(x) sum(!is.na(x)))  }
42  else
43     { if (length(sizes)==1)
44          sizes <- rep(sizes, nrow(data))
45       else if (length(sizes) != nrow(data))
46                stop("sizes length doesn't match with data") }
47
48  if (missing(labels))
49     { if (is.null(rownames(data))) labels <- 1:nrow(data)
50       else                         labels <- rownames(data) }
51
52  stats <- paste("stats.", type, sep = "")
53  if (!exists(stats, mode="function"))
54     stop(paste("function", stats, "is not defined"))
55  stats <- do.call(stats, list(data, sizes))
56  statistics <- stats$statistics
57  if (missing(center)) center <- stats$center
58
59  sd <- paste("sd.", type, sep = "")
60  if (!exists(sd, mode="function"))
61     stop(paste("function", sd, "is not defined!"))
62  missing.std.dev <- missing(std.dev)
63  if (missing.std.dev)
64     { std.dev <- NULL
65       std.dev <- switch(type,
66                         "xbar" = { if(any(sizes > 25)) "RMSDF"
67                                    else                "UWAVE-R" },
68                         "xbar.one" = "MR",
69                         "R" = "UWAVE-R",
70                         "S" = "UWAVE-SD",
71                         NULL)
72       std.dev <- do.call(sd, list(data, sizes, std.dev)) }
73  else
74     { if (is.character(std.dev))
75          { std.dev <- do.call(sd, list(data, sizes, std.dev)) }
76       else
77          { if (!is.numeric(std.dev))
78               stop("if provided the argument 'std.dev' must be a method available or a numerical value. See help(qcc).")  }
79     }
80
81  names(statistics) <-  rownames(data) <-  labels
82  names(dimnames(data)) <- list("Group", "Samples")
83
84  object <- list(call = call, type = type,
85                 data.name = data.name, data = data,
86                 statistics = statistics, sizes = sizes,
87                 center = center, std.dev = std.dev)
88  # check for new data provided and update object
89  if (!missing(newdata))
90     {   if (missing(newdata.name))
91			{newdata.name <- deparse(substitute(newdata))}
92       newdata <- data.matrix(newdata)
93       if (missing(newsizes))
94          { if (any(type==c("p", "np", "u")))
95               stop(paste("sample sizes must be given for a", type, "Chart"))
96            else
97               newsizes <- apply(newdata, 1, function(x) sum(!is.na(x))) }
98       else
99          { if (length(newsizes)==1)
100               newsizes <- rep(newsizes, nrow(newdata))
101            else if (length(newsizes) != nrow(newdata))
102                     stop("newsizes length doesn't match with newdata") }
103       stats <- paste("stats.", type, sep = "")
104       if (!exists(stats, mode="function"))
105          stop(paste("function", stats, "is not defined"))
106       newstats <- do.call(stats, list(newdata, newsizes))$statistics
107       if (missing(newlabels))
108          { if (is.null(rownames(newdata)))
109               { start <- length(statistics)
110                 newlabels <- seq(start+1, start+length(newstats)) }
111            else
112               { newlabels <- rownames(newdata) }
113          }
114       names(newstats) <- newlabels
115       object$newstats <- newstats
116       object$newdata  <- newdata
117       object$newsizes <- newsizes
118       object$newdata.name <- newdata.name
119       statistics <- c(statistics, newstats)
120       sizes <- c(sizes, newsizes)
121     }
122
123  conf <- nsigmas
124  if (!missing(confidence.level))
125     conf <- confidence.level
126  if (conf >= 1)
127     { object$nsigmas <- conf }
128  else
129     if (conf > 0 & conf < 1)
130        { object$confidence.level <- conf }
131
132  # get control limits
133  if (missing(limits))
134     { limits <- paste("limits.", type, sep = "")
135       if (!exists(limits, mode="function"))
136          stop(paste("function", limits, "is not defined"))
137       limits <- do.call(limits, list(center = center, std.dev = std.dev,
138                                      sizes = sizes, conf = conf))
139     }
140  else
141     { if (!missing.std.dev)
142          warning("'std.dev' is not used when limits is given")
143       if (!is.numeric(limits))
144          stop("'limits' must be a vector of length 2 or a 2-columns matrix")
145       limits <- matrix(limits, ncol = 2)
146       dimnames(limits) <- list(rep("",nrow(limits)), c("LCL ", "UCL"))
147     }
148
149  lcl <- limits[,1]
150  ucl <- limits[,2]
151  object$limits <- limits
152  if (is.function(rules)) violations <- rules(object)
153  else                    violations <- NULL
154  object$violations <- violations
155
156  class(object) <- "qcc"
157  if(plot) plot(object, ...)
158  return(object)
159}
160
161print.qcc <- function(x, ...) str(x,1)
162
163summary.qcc <- function(object, digits =  getOption("digits"), ...)
164{
165  #object <- x   # Argh.  Really want to use 'object' anyway
166  cat("\nCall:\n",deparse(object$call),"\n\n",sep="")
167  data.name <- object$data.name
168  type <- object$type
169  cat(paste(type, "chart for", data.name, "\n"))
170  statistics <- object$statistics
171  cat("\nSummary of group statistics:\n")
172  print(summary(statistics), digits = digits, ...)
173  sizes <- object$sizes
174  if(length(unique(sizes))==1)
175     sizes <- sizes[1]
176  if(length(sizes) == 1)
177     cat("\nGroup sample size: ", format(sizes))
178  else {
179         cat("\nSummary of group sample sizes: ")
180         tab <- table(sizes)
181         print(matrix(c(as.numeric(names(tab)), tab),
182                      ncol = length(tab), byrow = TRUE,
183                      dimnames = list(c("  sizes", "  counts"),
184                      character(length(tab)))),
185               digits = digits, ...)
186        }
187  cat("\nNumber of groups: ", length(statistics))
188
189  center <- object$center
190  if(length(center) == 1)
191    { cat("\nCenter of group statistics: ", format(center, digits = digits)) }
192  else
193    { out <- paste(format(center, digits = digits))
194      out <- out[which(cumsum(nchar(out)+1) < getOption("width")-40)]
195      out <- paste0(paste(out, collapse = " "), " ...")
196      cat("\nCenter of group statistics: ", out, sep = "")
197    }
198
199  sd <- object$std.dev
200  if(length(sd) == 1)
201    { cat("\nStandard deviation: ", format(sd, digits = digits), "\n") }
202  else
203    { out <- paste(format(sd, digits = digits))
204      out <- out[which(cumsum(nchar(out)+1) < getOption("width")-40)]
205      out <- paste0(paste(out, collapse = " "), " ...")
206      cat("\nStandard deviation: ", out, "\n", sep = "")
207    }
208
209  newdata.name <- object$newdata.name
210  newstats <- object$newstats
211  if (!is.null(newstats))
212     { cat(paste("\nSummary of group statistics in ",
213                 newdata.name, ":\n", sep = ""))
214       print(summary(newstats), digits = digits, ...)
215       newsizes <- object$newsizes
216       if (length(unique(newsizes)) == 1)
217          newsizes <- newsizes[1]
218       if (length(newsizes) == 1)
219             cat("\nGroup sample size: ", format(newsizes))
220       else
221           { cat("\nSummary of group sample sizes:\n")
222             new.tab <- table(newsizes)
223             print(matrix(c(as.numeric(names(new.tab)), new.tab),
224                          ncol = length(new.tab), byrow = TRUE,
225                          dimnames = list(c("  sizes", "  counts"),
226                                          character(length(new.tab)))),
227                   digits = digits, ...)
228           }
229       cat("\nNumber of groups: ", length(newstats), "\n")
230     }
231
232  limits <- object$limits
233  if (!is.null(limits))
234     { cat("\nControl limits:\n")
235       .printShortMatrix(limits, digits = digits, ...) }
236
237  invisible()
238}
239
240
241plot.qcc <- function(x, add.stats = TRUE, chart.all = TRUE,
242                     label.limits = c("LCL ", "UCL"),
243                     title, xlab, ylab, ylim, axes.las = 0,
244                     digits =  getOption("digits"),
245                     restore.par = TRUE, ...)
246{
247  object <- x  # Argh.  Really want to use 'object' anyway
248  if ((missing(object)) | (!inherits(object, "qcc")))
249     stop("an object of class `qcc' is required")
250
251  # collect info from object
252  type <- object$type
253  std.dev <- object$std.dev
254  data.name <- object$data.name
255  center <- object$center
256  stats <- object$statistics
257  limits <- object$limits
258  lcl <- limits[,1]
259  ucl <- limits[,2]
260  newstats <- object$newstats
261  newdata.name <- object$newdata.name
262  violations <- object$violations
263  if(chart.all)
264    { statistics <- c(stats, newstats)
265      indices <- 1:length(statistics) }
266  else
267     { if(is.null(newstats))
268         { statistics <- stats
269           indices <- 1:length(statistics) }
270       else
271         { statistics <- newstats
272           indices <- seq(length(stats)+1, length(stats)+length(newstats)) }
273     }
274
275  if (missing(title))
276     { if (is.null(newstats))
277            main.title <- paste(type, "Chart\nfor", data.name)
278       else if (chart.all)
279                 main.title <- paste(type, "Chart\nfor", data.name,
280                                     "and", newdata.name)
281            else main.title <- paste(type, "Chart\nfor", newdata.name)
282     }
283  else main.title <- paste(title)
284
285  oldpar <- par(no.readonly = TRUE)
286  if(restore.par) on.exit(par(oldpar))
287  mar <- pmax(oldpar$mar, c(4.1,4.1,3.1,2.1))
288  par(bg  = qcc.options("bg.margin"),
289      cex = oldpar$cex * qcc.options("cex"),
290      mar = if(add.stats) pmax(mar, c(7.6,0,0,0)) else mar)
291
292  # plot Shewhart chart
293  plot(indices, statistics, type="n",
294       ylim = if(!missing(ylim)) ylim
295              else range(statistics, limits, center),
296       ylab = if(missing(ylab)) "Group summary statistics" else ylab,
297       xlab = if(missing(xlab)) "Group" else xlab,
298       axes = FALSE)
299  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
300       col = qcc.options("bg.figure"))
301  axis(1, at = indices, las = axes.las,
302       labels = if(is.null(names(statistics)))
303                   as.character(indices) else names(statistics))
304  axis(2, las = axes.las)
305  box()
306  top.line <- par("mar")[3]-length(capture.output(cat(main.title)))
307  top.line <- top.line - if(chart.all & (!is.null(newstats))) 0.1 else 0.5
308  mtext(main.title, side = 3, line = top.line,
309        font = par("font.main"),
310        cex  = qcc.options("cex"),
311        col  = par("col.main"))
312
313  lines(indices, statistics, type = "b", pch=20)
314
315  if(length(center) == 1)
316       abline(h = center)
317  else lines(indices, center[indices], type="s")
318
319  if(length(lcl) == 1)
320    { abline(h = lcl, lty = 2)
321      abline(h = ucl, lty = 2) }
322  else
323    { lines(indices, lcl[indices], type="s", lty = 2)
324      lines(indices, ucl[indices], type="s", lty = 2) }
325  mtext(label.limits, side = 4, at = c(rev(lcl)[1], rev(ucl)[1]),
326        las = 1, line = 0.1, col = gray(0.3), cex = par("cex"))
327  mtext("CL", side = 4, at = rev(center)[1],
328        las = 1, line = 0.1, col = gray(0.3), cex = par("cex"))
329
330  if(is.null(qcc.options("violating.runs")))
331     stop(".qcc.options$violating.runs undefined. See help(qcc.options).")
332  if(length(violations$violating.runs))
333    { v <- violations$violating.runs
334      if(!chart.all & !is.null(newstats))
335        { v <- v - length(stats)
336          v <- v[v>0] }
337      points(indices[v], statistics[v],
338             col = qcc.options("violating.runs")$col,
339             pch = qcc.options("violating.runs")$pch)
340    }
341
342  if(is.null(qcc.options("beyond.limits")))
343     stop(".qcc.options$beyond.limits undefined. See help(qcc.options).")
344  if(length(violations$beyond.limits))
345    { v <- violations$beyond.limits
346      if(!chart.all & !is.null(newstats))
347        { v <- v - length(stats)
348          v <- v[v>0] }
349      points(indices[v], statistics[v],
350             col = qcc.options("beyond.limits")$col,
351             pch = qcc.options("beyond.limits")$pch)
352    }
353
354  if(chart.all & (!is.null(newstats)))
355    { len.obj.stats <- length(object$statistics)
356      len.new.stats <- length(statistics) - len.obj.stats
357      abline(v = len.obj.stats + 0.5, lty = 3)
358      mtext(# paste("Calibration data in", data.name),
359            "Calibration data", cex = par("cex")*0.8,
360            at = len.obj.stats/2, line = 0, adj = 0.5)
361      mtext(# paste("New data in", object$newdata.name),
362            "New data", cex = par("cex")*0.8,
363            at = len.obj.stats + len.new.stats/2, line = 0, adj = 0.5)
364    }
365
366  if(add.stats)
367    {
368      # computes the x margins of the figure region
369      plt <- par()$plt; usr <- par()$usr
370      px <- diff(usr[1:2])/diff(plt[1:2])
371      xfig <- c(usr[1]-px*plt[1], usr[2]+px*(1-plt[2]))
372      at.col <- xfig[1] + diff(xfig[1:2])*c(0.10, 0.40, 0.65)
373      top.line <- 4.5
374      # write info at bottom
375      mtext(paste("Number of groups = ", length(statistics), sep = ""),
376            side = 1, line = top.line, adj = 0, at = at.col[1],
377            font = qcc.options("font.stats"),
378            cex = par("cex")*qcc.options("cex.stats"))
379      center <- object$center
380      if(length(center) == 1)
381        { mtext(paste("Center = ", signif(center[1], digits), sep = ""),
382                side = 1, line = top.line+1, adj = 0, at = at.col[1],
383                font = qcc.options("font.stats"),
384                cex = par("cex")*qcc.options("cex.stats"))
385        }
386      else
387        { mtext("Center is variable",
388                side = 1, line = top.line+1, adj = 0, at = at.col[1],
389                font = qcc.options("font.stats"),
390                cex = par("cex")*qcc.options("cex.stats"))
391        }
392
393      if(length(std.dev) == 1)
394        { mtext(paste("StdDev = ", signif(std.dev, digits), sep = ""),
395                side = 1, line = top.line+2, adj = 0, at = at.col[1],
396                font = qcc.options("font.stats"),
397                cex = par("cex")*qcc.options("cex.stats"))
398        }
399      else
400        { mtext("StdDev is variable",
401                side = 1, line = top.line+2, adj = 0, at = at.col[1],
402                font = qcc.options("font.stats"),
403                cex = par("cex")*qcc.options("cex.stats"))
404        }
405
406      if(length(unique(lcl)) == 1)
407        { mtext(paste("LCL = ", signif(lcl[1], digits), sep = ""),
408                side = 1, line = top.line+1, adj = 0, at = at.col[2],
409                font = qcc.options("font.stats"),
410                cex = par("cex")*qcc.options("cex.stats"))
411        }
412      else
413        { mtext("LCL is variable",
414                side = 1, line = top.line+1, adj = 0, at = at.col[2],
415                font = qcc.options("font.stats"),
416                cex = par("cex")*qcc.options("cex.stats"))
417        }
418
419      if(length(unique(ucl)) == 1)
420         { mtext(paste("UCL = ", signif(ucl[1], digits), sep = ""),
421                 side = 1, line = top.line+2, adj = 0, at = at.col[2],
422                 font = qcc.options("font.stats"),
423                 cex = par("cex")*qcc.options("cex.stats"))
424         }
425       else
426         { mtext("UCL is variable",
427                 side = 1, line = top.line+2, adj = 0, at = at.col[2],
428                 font = qcc.options("font.stats"),
429                 cex = par("cex")*qcc.options("cex.stats"))
430         }
431
432       if(!is.null(violations))
433         { mtext(paste("Number beyond limits =",
434                       length(unique(violations$beyond.limits))),
435                 side = 1, line = top.line+1, adj = 0, at = at.col[3],
436                 font = qcc.options("font.stats"),
437                 cex = par("cex")*qcc.options("cex.stats"))
438           mtext(paste("Number violating runs =",
439                       length(unique(violations$violating.runs))),
440                 side = 1, line = top.line+2, adj = 0, at = at.col[3],
441                 font = qcc.options("font.stats"),
442                 cex = par("cex")*qcc.options("cex.stats"))
443         }
444     }
445
446  invisible()
447}
448
449#
450#  Functions used to compute Shewhart charts statistics
451#
452
453.qcc.c4 <- function(n)
454{ sqrt(2/(n - 1)) * exp(lgamma(n/2) - lgamma((n - 1)/2)) }
455
456# xbar
457
458stats.xbar <- function(data, sizes)
459{
460  if (missing(sizes))
461     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
462  statistics <- apply(data, 1, mean, na.rm=TRUE)
463  center <- sum(sizes * statistics)/sum(sizes)
464  list(statistics = statistics, center = center)
465}
466
467sd.xbar <- function(data, sizes, std.dev = c("UWAVE-R", "UWAVE-SD", "MVLUE-R", "MVLUE-SD", "RMSDF"))
468{
469  if (!is.numeric(std.dev))
470     std.dev <- match.arg(std.dev)
471  if (missing(sizes))
472     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
473  if (any(sizes == 1))
474     stop("group sizes must be larger than one")
475  c4 <- .qcc.c4
476  if (is.numeric(std.dev))
477     { sd <- std.dev }
478  else
479     { switch(std.dev,
480              "UWAVE-R" = {  R <- apply(data, 1, function(x)
481                                                 diff(range(x, na.rm = TRUE)))
482                             d2 <- qcc.options("exp.R.unscaled")[sizes]
483                             sd <- sum(R/d2)/length(sizes) },
484              "UWAVE-SD" = { S <- apply(data, 1, sd, na.rm = TRUE)
485                             sd <- sum(S/c4(sizes))/length(sizes) },
486              "MVLUE-R"  = { R <- apply(data, 1, function(x)
487                                                 diff(range(x, na.rm = TRUE)))
488                             d2 <- qcc.options("exp.R.unscaled")[sizes]
489                             d3 <- qcc.options("se.R.unscaled")[sizes]
490                             w  <- (d2/d3)^2
491                             sd <- sum(R/d2*w)/sum(w) },
492              "MVLUE-SD" = { S <- apply(data, 1, sd, na.rm = TRUE)
493                             w  <- c4(sizes)^2/(1-c4(sizes)^2)
494                             sd <- sum(S/c4(sizes)*w)/sum(w) },
495              "RMSDF" =    { S <- apply(data, 1, sd, na.rm = TRUE)
496                             w  <- sizes-1
497                             sd <- sqrt(sum(S^2*w)/sum(w))/c4(sum(w)+1) }
498              )
499     }
500#  if (missing(std.dev))
501#     var.within <- apply(data, 1, var, na.rm=TRUE)
502#  else
503#     var.within <- std.dev^2
504#  var.df <- sum(sizes - 1)
505#  if (equal.sd)
506#     { std.dev <- sqrt(sum((sizes - 1) * var.within)/var.df) / c4(var.df + 1) }
507#  else
508#     { c <- c4(sizes)/(1 - c4(sizes)^2)
509#       std.dev <- sum(c * sqrt(var.within))/sum(c * c4(sizes)) }
510  return(sd)
511}
512
513limits.xbar <- function(center, std.dev, sizes, conf)
514{
515  if (length(unique(sizes))==1)
516     sizes <- sizes[1]
517  se.stats <- std.dev/sqrt(sizes)
518  if (conf >= 1)
519     { lcl <- center - conf * se.stats
520       ucl <- center + conf * se.stats
521     }
522  else
523     { if (conf > 0 & conf < 1)
524          { nsigmas <- qnorm(1 - (1 - conf)/2)
525            lcl <- center - nsigmas * se.stats
526            ucl <- center + nsigmas * se.stats
527          }
528       else stop("invalid 'conf' argument. See help.")
529     }
530  limits <- matrix(c(lcl, ucl), ncol = 2)
531  rownames(limits) <- rep("", length = nrow(limits))
532  colnames(limits) <- c("LCL", "UCL")
533  return(limits)
534}
535
536# S chart
537
538stats.S <- function(data, sizes)
539{
540  if (missing(sizes))
541     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
542  if(ncol(data)==1)
543    { statistics <- as.vector(data) }
544  else
545    { statistics <- sqrt(apply(data, 1, var, na.rm=TRUE)) }
546  if (length(sizes == 1))
547     sizes <- rep(sizes, length(statistics))
548  center <- sum(sizes * statistics)/sum(sizes)
549  list(statistics = statistics, center = center)
550}
551
552sd.S <- function(data, sizes, std.dev = c("UWAVE-SD", "MVLUE-SD", "RMSDF"))
553{
554  if (!is.numeric(std.dev))
555     std.dev <- match.arg(std.dev)
556  sd.xbar(data, sizes, std.dev)
557}
558
559limits.S <- function(center, std.dev, sizes, conf)
560{
561  if (length(unique(sizes))==1)
562     sizes <- sizes[1]
563  c4 <- .qcc.c4
564  se.stats <- std.dev * sqrt(1 - c4(sizes)^2)
565  if (conf >= 1)
566     { lcl <- pmax(0, center - conf * se.stats)
567       ucl <- center + conf * se.stats
568     }
569  else
570     { if (conf > 0 & conf < 1)
571          { ucl <- std.dev * sqrt(qchisq(1 - (1 - conf)/2, sizes - 1)/
572                                  (sizes - 1))
573            lcl <- std.dev * sqrt(qchisq((1 - conf)/2, sizes - 1)/
574                                  (sizes - 1))
575          }
576          else stop("invalid conf argument. See help.")
577     }
578  limits <- matrix(c(lcl, ucl), ncol = 2)
579  rownames(limits) <- rep("", length = nrow(limits))
580  colnames(limits) <- c("LCL", "UCL")
581  limits
582}
583
584# R Chart
585
586stats.R <- function(data, sizes)
587{
588  if (missing(sizes))
589     sizes <- apply(data, 1, function(x) sum(!is.na(x)))
590  if(ncol(data)==1)
591    { statistics <- as.vector(data) }
592  else
593    { statistics <- apply(data, 1, function(x) diff(range(x, na.rm=TRUE))) }
594  if (length(sizes == 1))
595     sizes <- rep(sizes, length(statistics))
596  center <- sum(sizes * statistics)/sum(sizes)
597  list(statistics = statistics, center = center)
598}
599
600sd.R <- function(data, sizes, std.dev = c("UWAVE-R", "MVLUE-R"))
601{
602  if (!is.numeric(std.dev))
603     std.dev <- match.arg(std.dev)
604  sd.xbar(data, sizes, std.dev)
605}
606
607limits.R <- function(center, std.dev, sizes, conf)
608{
609  if (length(unique(sizes))==1)
610     sizes <- sizes[1]
611  se.R.unscaled <- qcc.options("se.R.unscaled")
612  Rtab <- length(se.R.unscaled)
613  if (conf >= 1)
614     { if (any(sizes > Rtab))
615          stop(paste("group size must be less than",
616                      Rtab + 1, "when giving nsigmas"))
617       se.R <- se.R.unscaled[sizes] * std.dev
618       lcl <- pmax(0, center - conf * se.R)
619       ucl <- center + conf * se.R
620     }
621  else
622     { if (conf > 0 && conf < 1)
623          { ucl <- qtukey(1 - (1 - conf)/2, sizes, 1e100) * std.dev
624            lcl <- qtukey((1 - conf)/2, sizes, 1e100) * std.dev
625          }
626       else stop("invalid conf argument. See help.")
627     }
628  limits <- matrix(c(lcl, ucl), ncol = 2)
629  rownames(limits) <- rep("", length = nrow(limits))
630  colnames(limits) <- c("LCL", "UCL")
631  return(limits)
632}
633
634# xbar Chart for one-at-time data
635
636stats.xbar.one <- function(data, sizes)
637{
638  statistics <- as.vector(data)
639  center <- mean(statistics)
640  list(statistics = statistics, center = center)
641}
642
643sd.xbar.one <- function(data, sizes, std.dev = c("MR", "SD"), k = 2)
644{
645  data <- as.vector(data)
646  n <- length(data)
647  if(!is.numeric(std.dev))
648     std.dev <- match.arg(std.dev)
649  c4 <- .qcc.c4
650  if(is.numeric(std.dev))
651    { sd <- std.dev }
652  else
653    { switch(std.dev,
654             "MR" = { d2 <- qcc.options("exp.R.unscaled")
655                      if(is.null(d2))
656                         stop(".qcc.options$exp.R.unscaled is null")
657                      d <- 0
658                      for(j in k:n)
659                          d <- d+abs(diff(range(data[c(j:(j-k+1))])))
660                      sd <- (d/(n-k+1))/d2[k] },
661             "SD" = { sd <- sd(data)/c4(n) },
662             sd <- NULL)
663    }
664  return(sd)
665}
666
667
668limits.xbar.one <- function(center, std.dev, sizes, conf)
669{
670  se.stats <- std.dev
671  if (conf >= 1)
672     { lcl <- center - conf * se.stats
673       ucl <- center + conf * se.stats
674     }
675  else
676     { if (conf > 0 & conf < 1)
677          { nsigmas <- qnorm(1 - (1 - conf)/2)
678            lcl <- center - nsigmas * se.stats
679            ucl <- center + nsigmas * se.stats
680          }
681       else stop("invalid conf argument. See help.")
682     }
683  limits <- matrix(c(lcl, ucl), ncol = 2)
684  rownames(limits) <- rep("", length = nrow(limits))
685  colnames(limits) <- c("LCL", "UCL")
686  return(limits)
687}
688
689
690# p Chart
691
692stats.p <- function(data, sizes)
693{
694  data <- as.vector(data)
695  sizes <- as.vector(sizes)
696  pbar <- sum(data)/sum(sizes)
697  list(statistics = data/sizes, center = pbar)
698}
699
700sd.p <- function(data, sizes, ...)
701{
702  data <- as.vector(data)
703  sizes <- as.vector(sizes)
704  pbar <- sum(data)/sum(sizes)
705  std.dev <- sqrt(pbar * (1 - pbar))
706  return(std.dev)
707}
708
709limits.p <- function(center, std.dev, sizes, conf)
710{
711  limits.np(center * sizes, std.dev, sizes, conf) / sizes
712}
713
714# np Chart
715
716stats.np <- function(data, sizes)
717{
718  data <- as.vector(data)
719  sizes <- as.vector(sizes)
720  pbar <- sum(data)/sum(sizes)
721  center <- sizes * pbar
722  if (length(unique(center)) == 1)
723     center <- center[1]
724  list(statistics = data, center = center)
725}
726
727sd.np <- function(data, sizes, ...)
728{
729  data <- as.vector(data)
730  sizes <- as.vector(sizes)
731  pbar <- sum(data)/sum(sizes)
732  std.dev <- sqrt(sizes * pbar * (1 - pbar))
733  if (length(unique(std.dev)) == 1)
734     std.dev <- std.dev[1]
735  return(std.dev)
736}
737
738limits.np <- function(center, std.dev, sizes, conf)
739{
740  sizes <- as.vector(sizes)
741  if (length(unique(sizes)) == 1)
742     sizes <- sizes[1]
743  pbar <- mean(center / sizes)
744  if (conf >= 1)
745     { tol <- conf * sqrt(pbar * (1 - pbar) * sizes)
746       lcl <- pmax(center - tol, 0)
747       ucl <- pmin(center + tol, sizes)
748     }
749  else
750     { if (conf > 0 & conf < 1)
751          { lcl <- qbinom((1 - conf)/2, sizes, pbar)
752            ucl <- qbinom((1 - conf)/2, sizes, pbar, lower.tail = FALSE)
753          }
754       else stop("invalid conf argument. See help.")
755     }
756  limits <- matrix(c(lcl, ucl), ncol = 2)
757  rownames(limits) <- rep("", length = nrow(limits))
758  colnames(limits) <- c("LCL", "UCL")
759  return(limits)
760}
761
762# c Chart
763
764stats.c <- function(data, sizes)
765{
766  data <- as.vector(data)
767  sizes <- as.vector(sizes)
768  if (length(unique(sizes)) != 1)
769     stop("all sizes must be be equal for a c chart")
770  statistics <- data
771  center <- mean(statistics)
772  list(statistics = statistics, center = center)
773}
774
775sd.c <- function(data, sizes, ...)
776{
777  data <- as.vector(data)
778  std.dev <- sqrt(mean(data))
779  return(std.dev)
780}
781
782limits.c <- function(center, std.dev, sizes, conf)
783{
784  if (conf >= 1)
785     { lcl <- center - conf * sqrt(center)
786       lcl[lcl < 0] <- 0
787       ucl <- center + conf * sqrt(center)
788     }
789  else
790     { if (conf > 0 & conf < 1)
791          { ucl <- qpois(1 - (1 - conf)/2, center)
792            lcl <- qpois((1 - conf)/2, center)
793          }
794       else stop("invalid conf argument. See help.")
795     }
796  limits <- matrix(c(lcl, ucl), ncol = 2)
797  rownames(limits) <- rep("", length = nrow(limits))
798  colnames(limits) <- c("LCL", "UCL")
799  return(limits)
800}
801
802# u Chart
803
804stats.u <- function(data, sizes)
805{
806  data <- as.vector(data)
807  sizes <- as.vector(sizes)
808  statistics <- data/sizes
809  center <- sum(sizes * statistics)/sum(sizes)
810  list(statistics = statistics, center = center)
811}
812
813sd.u <- function(data, sizes, ...)
814{
815  data <- as.vector(data)
816  sizes <- as.vector(sizes)
817  std.dev <- sqrt(sum(data)/sum(sizes))
818  return(std.dev)
819}
820
821limits.u <- function(center, std.dev, sizes, conf)
822{
823  sizes <- as.vector(sizes)
824  if (length(unique(sizes))==1)
825     sizes <- sizes[1]
826  limits.c(center * sizes, std.dev, sizes, conf) / sizes
827}
828
829#
830# Functions used to signal points out of control
831#
832
833shewhart.rules <- function(object, limits = object$limits, run.length = qcc.options("run.length"))
834{
835# Return a list of cases beyond limits and violating runs
836  bl <- beyond.limits(object, limits = limits)
837  vr <- violating.runs(object, run.length = run.length)
838  list(beyond.limits = bl, violating.runs = vr)
839}
840
841beyond.limits <- function(object, limits = object$limits)
842{
843# Return cases beyond limits
844  statistics <- c(object$statistics, object$newstats)
845  lcl <- limits[,1]
846  ucl <- limits[,2]
847  index.above.ucl <- seq(along = statistics)[statistics > ucl]
848  index.below.lcl <- seq(along = statistics)[statistics < lcl]
849  return(c(index.above.ucl, index.below.lcl))
850}
851
852violating.runs <- function(object, run.length = qcc.options("run.length"))
853{
854# Return indices of points violating runs
855  if(run.length == 0)
856    return(numeric())
857  center <- object$center
858  statistics <- c(object$statistics, object$newstats)
859  cl <- object$limits
860  diffs <- statistics - center
861  diffs[diffs > 0] <- 1
862  diffs[diffs < 0] <- -1
863  runs <- rle(diffs)
864  vruns <- rep(runs$lengths >= run.length, runs$lengths)
865  vruns.above <- (vruns & (diffs > 0))
866  vruns.below <- (vruns & (diffs < 0))
867  rvruns.above <- rle(vruns.above)
868  rvruns.below <- rle(vruns.below)
869  vbeg.above <- cumsum(rvruns.above$lengths)[rvruns.above$values] -
870                      (rvruns.above$lengths - run.length)[rvruns.above$values]
871  vend.above <- cumsum(rvruns.above$lengths)[rvruns.above$values]
872  vbeg.below <- cumsum(rvruns.below$lengths)[rvruns.below$values] -
873                      (rvruns.below$lengths - run.length)[rvruns.below$values]
874  vend.below <- cumsum(rvruns.below$lengths)[rvruns.below$values]
875  violators <- numeric()
876  if (length(vbeg.above))
877     { for (i in 1:length(vbeg.above))
878           violators <- c(violators, vbeg.above[i]:vend.above[i]) }
879  if (length(vbeg.below))
880     { for (i in 1:length(vbeg.below))
881           violators <- c(violators, vbeg.below[i]:vend.below[i]) }
882  return(violators)
883}
884
885#-------------------------------------------------------------------#
886#                                                                   #
887#          Operating Characteristic Function                        #
888#                                                                   #
889#-------------------------------------------------------------------#
890
891oc.curves <- function(object, ...)
892{
893# Draws the operating characteristic curves for the qcc object
894
895  if ((missing(object)) | (!inherits(object, "qcc")))
896     stop("an object of class 'qcc' is required")
897
898  size <- unique(object$sizes)
899  if (length(size)>1)
900     stop("Operating characteristic curves available only for equal sample sizes!")
901
902  beta <- switch(object$type,
903                 xbar = oc.curves.xbar(object, ...),
904                 R    = oc.curves.R(object, ...),
905                 S    = oc.curves.S(object, ...),
906                 np   =,
907                 p    = oc.curves.p(object, ...),
908                 u    =,
909                 c    = oc.curves.c(object, ...))
910  if (is.null(beta))
911     stop("Operating characteristic curves not available for this type of chart.")
912
913  invisible(beta)
914}
915
916oc.curves.xbar <- function(object, n, c = seq(0, 5, length=101), nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
917{
918# Draw the operating-characteristic curves for the xbar-chart with nsigmas
919# limits. The values on the vertical axis give the probability of not detecting
920# a shift of c*sigma in the mean on the first sample following the shift.
921
922  if (!(object$type=="xbar"))
923     stop("not a `qcc' object of type \"xbar\".")
924
925  size <- unique(object$sizes)
926  if (length(size) > 1)
927     stop("Operating characteristic curves available only for equal sample sizes!")
928  if (missing(n))
929     n <- unique(c(size, c(1,5,10,15,20)))
930  if (is.null(nsigmas))
931     nsigmas <- qnorm(1 - (1 - object$confidence.level) / 2)
932
933  beta <- matrix(as.double(NA), length(n), length(c))
934  for (i in 1:length(n))
935      beta[i,] <- pnorm(nsigmas-c*sqrt(n[i])) - pnorm(-nsigmas-c*sqrt(n[i]))
936  rownames(beta) <- paste("n=",n,sep="")
937  colnames(beta) <- c
938
939  oldpar <- par(no.readonly = TRUE)
940  if(restore.par) on.exit(par(oldpar))
941  par(bg  = qcc.options("bg.margin"),
942      cex = oldpar$cex * qcc.options("cex"),
943      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))
944
945  plot(c, beta[1,], type="n",
946       ylim = c(0,1), xlim = c(0,max(c)),
947       xlab = "Process shift (std.dev)",
948       ylab = "Prob. type II error ")
949  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
950       col = qcc.options("bg.figure"))
951  box()
952  mtext(paste("OC curves for", object$type, "Chart"),
953        side = 3, line = par("mar")[3]/3,
954        font = par("font.main"),
955        cex  = qcc.options("cex"),
956        col  = par("col.main"))
957  for(i in 1:length(n))
958     lines(c, beta[i,], type = "l", lty=i)
959  beta <- t(beta)
960  names(dimnames(beta)) <- c("shift (std.dev)", "sample size")
961
962  if (identify)
963     { cs <- rep(c,length(n))
964       betas <- as.vector(beta)
965       labels <- paste("c=", formatC(cs, 2, flag="-"),
966                       ": beta=", formatC(betas, 4, flag="-"),
967                       ", ARL=", formatC(1/(1-betas), 2, flag="-"), sep="")
968       i <- identify(cs, betas, labels, pos=4, offset=0.2)
969       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
970     }
971  else
972     { legend(max(c), 1, legend = paste("n =", n),
973              bg = qcc.options("bg.figure"),
974              lty = 1:length(n), xjust = 1, yjust = 1)
975     }
976  invisible(beta)
977}
978
979oc.curves.p <- function(object, nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
980{
981  if (!(object$type=="p" | object$type=="np"))
982     stop("not a `qcc' object of type \"p\" or \"np\".")
983
984  size <- unique(object$sizes)
985  if (length(size) > 1)
986     stop("Operating characteristic curves available only for equal sample sizes!")
987
988  if (is.null(object$limits))
989     stop("the `qcc' object does not have control limits!")
990  limits <- object$limits
991  p <- seq(0, 1, length=101)
992
993  if (object$type=="p")
994     { UCL <- min(floor(size*limits[,2]), size)
995       LCL <- max(floor(size*limits[,1]), 0) }
996  else
997     { UCL <- min(floor(limits[,2]), size)
998       LCL <- max(floor(limits[,1]), 0) }
999  beta <- pbinom(UCL, size, p) - pbinom(LCL-1, size, p)
1000  names(beta) <- p
1001
1002  oldpar <- par(no.readonly = TRUE)
1003  if(restore.par) on.exit(par(oldpar))
1004  par(bg  = qcc.options("bg.margin"),
1005      cex = oldpar$cex * qcc.options("cex"),
1006      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))
1007
1008  plot(p, beta, type = "n",
1009       ylim = c(0,1), xlim = c(0,1),
1010       xlab = expression(p),
1011       ylab = "Prob. type II error ")
1012  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
1013       col = qcc.options("bg.figure"))
1014  box()
1015  mtext(paste("OC curves for", object$type, "Chart"),
1016        side = 3, line = par("mar")[3]/3,
1017        font = par("font.main"),
1018        cex  = qcc.options("cex"),
1019        col  = par("col.main"))
1020
1021  lines(p, beta)
1022  lines(rep(p[which.max(beta)], 2), c(0, max(beta)), lty = 2)
1023
1024  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.")
1025
1026  if (identify)
1027     { labels <- paste("p=", formatC(p, 2, flag="-"),
1028                       ": beta=", formatC(beta, 4, flag="-"),
1029                       ", ARL=", formatC(1/(1-beta), 2, flag="-"), sep="")
1030       i <- identify(p, beta, labels, pos=4, offset=0.2)
1031       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
1032     }
1033  invisible(beta)
1034}
1035
1036oc.curves.c <- function(object, nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
1037{
1038  type <- object$type
1039  if (!(object$type=="c" | object$type=="u"))
1040     stop("not a `qcc' object of type \"c\" or \"u\".")
1041
1042  size <- unique(object$sizes)
1043  if (length(size) > 1)
1044     stop("Operating characteristic curves available only for equal sample size!")
1045
1046  if (is.null(object$limits))
1047     stop("the `qcc' object does not have control limits!")
1048  limits <- object$limits
1049  CL  <- object$center
1050  std.dev <- object$std.dev
1051  if (object$type=="c")
1052     { max.lambda <- ceiling(CL+10*std.dev)
1053       UCL <- floor(limits[1,2])
1054       LCL <- floor(limits[1,1])
1055     }
1056  else
1057     { max.lambda <- ceiling(CL*size+10*std.dev)[1]
1058       UCL <- floor(size*limits[1,2])
1059       LCL <- floor(size*limits[1,1])
1060     }
1061  lambda <- seq(0, max.lambda)
1062  beta <- ppois(UCL, lambda) - ppois(LCL-1, lambda)
1063  names(beta) <- lambda
1064
1065  oldpar <- par(no.readonly = TRUE)
1066  if(restore.par) on.exit(par(oldpar))
1067  par(bg  = qcc.options("bg.margin"),
1068      cex = oldpar$cex * qcc.options("cex"),
1069      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))
1070
1071  plot(lambda, beta, type = "n",
1072       ylim = c(0,1), xlim = range(lambda),
1073       xlab = "Mean",
1074       ylab = "Prob. type II error ")
1075  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
1076       col = qcc.options("bg.figure"))
1077  box()
1078  mtext(paste("OC curves for", object$type, "Chart"),
1079        side = 3, line = par("mar")[3]/3,
1080        font = par("font.main"),
1081        cex  = qcc.options("cex"),
1082        col  = par("col.main"))
1083
1084  lines(lambda, beta)
1085  lines(rep(lambda[which.max(beta)], 2), c(0, max(beta)), lty = 2)
1086
1087  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.")
1088
1089  if (identify)
1090     { labels <- paste("lambda=", formatC(lambda, 0, flag="-"),
1091                       ": beta=", formatC(beta, 4, flag="-"),
1092                       ", ARL=", formatC(1/(1-beta), 2, flag="-"), sep="")
1093       i <- identify(lambda, beta, labels, pos=4, offset=0.2)
1094       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
1095     }
1096  invisible(beta)
1097}
1098
1099oc.curves.R <-
1100function(object, n, c = seq(1, 6, length=101), nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
1101{
1102# Draw the operating-characteristic curves for the R-chart with nsigmas
1103# limits. The values on the vertical axis give the probability of not detecting
1104# a change from sigma to c*sigma on the first sample following the change.
1105
1106  if (!(object$type=="R"))
1107     stop("not a `qcc' object of type \"R\".")
1108
1109  size <- unique(object$sizes)
1110  if (length(size) > 1)
1111     stop("Operating characteristic curves available only for equal sample sizes!")
1112  if (missing(n))
1113     n <- unique(c(size, c(2,5,10,15,20)))
1114  if (is.null(nsigmas))
1115    { tail.prob <- (1 - object$confidence.level) / 2
1116      beta.fun1 <- function(c, n, p)
1117      {
1118        lcl <- qtukey(p, n, Inf)
1119        ucl <- qtukey(p, n, Inf, lower.tail = FALSE)
1120        ptukey(ucl / c, n, Inf) - ptukey(lcl / c, n, Inf)
1121      }
1122      beta <- outer(c, n, beta.fun1, tail.prob)
1123    }
1124  else
1125    { exp.R.unscaled <- qcc.options("exp.R.unscaled")
1126      se.R.unscaled <- qcc.options("se.R.unscaled")
1127      Rtab <- min(length(exp.R.unscaled), length(se.R.unscaled))
1128      if (any(n > Rtab))
1129          stop(paste("group size must be less than",
1130                      Rtab + 1, "when giving nsigmas"))
1131      beta.fun2 <- function(c, n, conf)
1132      {
1133        d2 <- exp.R.unscaled[n]
1134        d3 <- se.R.unscaled[n]
1135        lcl <- pmax(0, d2 - conf * d3)
1136        ucl <- d2 + conf * d3
1137        ptukey(ucl / c, n, Inf) - ptukey(lcl / c, n, Inf)
1138      }
1139      beta <- outer(c, n, beta.fun2, nsigmas)
1140    }
1141
1142  colnames(beta) <- paste("n=",n,sep="")
1143  rownames(beta) <- c
1144
1145  oldpar <- par(no.readonly = TRUE)
1146  if(restore.par) on.exit(par(oldpar))
1147  par(bg  = qcc.options("bg.margin"),
1148      cex = oldpar$cex * qcc.options("cex"),
1149      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))
1150
1151  plot(c, beta[,1], type="n",
1152       ylim = c(0,1), xlim = c(1,max(c)),
1153       xlab = "Process scale multiplier",
1154       ylab = "Prob. type II error ")
1155  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
1156       col = qcc.options("bg.figure"))
1157  box()
1158  mtext(paste("OC curves for", object$type, "Chart"),
1159        side = 3, line = par("mar")[3]/3,
1160        font = par("font.main"),
1161        cex  = qcc.options("cex"),
1162        col  = par("col.main"))
1163
1164  matlines(c, beta, lty = 1:length(n), col = 1)
1165
1166  names(dimnames(beta)) <- c("scale multiplier", "sample size")
1167
1168  if (identify)
1169     { cs <- rep(c,length(n))
1170       betas <- as.vector(beta)
1171       labels <- paste("c=", formatC(cs, 2, flag="-"),
1172                       ": beta=", formatC(betas, 4, flag="-"),
1173                       ", ARL=", formatC(1/(1-betas), 2, flag="-"), sep="")
1174       i <- identify(cs, betas, labels, pos=4, offset=0.2)
1175       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
1176     }
1177  else
1178     { legend(max(c), 1, legend = paste("n =", n),
1179              bg = qcc.options("bg.figure"),
1180              lty = 1:length(n), xjust = 1, yjust = 1)
1181     }
1182  invisible(beta)
1183}
1184
1185oc.curves.S <- function(object, n, c = seq(1, 6, length=101), nsigmas = object$nsigmas, identify=FALSE, restore.par=TRUE)
1186{
1187# Draw the operating-characteristic curves for the S-chart with nsigmas
1188# limits. The values on the vertical axis give the probability of not detecting
1189# a change from sigma to c*sigma on the first sample following the change.
1190
1191  if (!(object$type=="S"))
1192     stop("not a `qcc' object of type \"S\".")
1193
1194  size <- unique(object$sizes)
1195  if (length(size) > 1)
1196     stop("Operating characteristic curves available only for equal sample sizes!")
1197  if (missing(n))
1198     n <- unique(c(size, c(2,5,10,15,20)))
1199  if (is.null(nsigmas))
1200    { tail.prob <- (1 - object$confidence.level) / 2
1201      beta.fun1 <- function(c, n, p)
1202      {
1203        ucl <- sqrt(qchisq(1 - p, n - 1) / (n - 1))
1204        lcl <- sqrt(qchisq(p, n - 1) / (n - 1))
1205        pchisq((n - 1) * (ucl / c)^2, n - 1) - pchisq((n - 1)* (lcl / c)^2, n - 1)
1206      }
1207      beta <- outer(c, n, beta.fun1, tail.prob)
1208    }
1209  else
1210    { c4 <- .qcc.c4
1211      beta.fun2 <- function(c, n)
1212      {
1213        center <- c4(n)
1214        tol <- sqrt(1 - c4(n)^2)
1215        lcl <- pmax(0, center - nsigmas * tol)
1216        ucl <- center + nsigmas * tol
1217        pchisq((n - 1) * (ucl / c)^2, n - 1) - pchisq((n - 1) * (lcl / c)^2, n - 1)
1218      }
1219      beta <- outer(c, n, beta.fun2)
1220    }
1221
1222  colnames(beta) <- paste("n=",n,sep="")
1223  rownames(beta) <- c
1224
1225  oldpar <- par(no.readonly = TRUE)
1226  if(restore.par) on.exit(par(oldpar))
1227  par(bg  = qcc.options("bg.margin"),
1228      cex = oldpar$cex * qcc.options("cex"),
1229      mar = pmax(oldpar$mar, c(4.1,4.1,2.1,2.1)))
1230
1231  plot(c, beta[,1], type="n",
1232       ylim = c(0,1), xlim = c(1,max(c)),
1233       xlab = "Process scale multiplier",
1234       ylab = "Prob. type II error ")
1235  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
1236       col = qcc.options("bg.figure"))
1237  box()
1238  mtext(paste("OC curves for", object$type, "Chart"),
1239        side = 3, line = par("mar")[3]/3,
1240        font = par("font.main"),
1241        cex  = qcc.options("cex"),
1242        col  = par("col.main"))
1243  matlines(c, beta, lty = 1:length(n), col = 1)
1244
1245  names(dimnames(beta)) <- c("scale multiplier", "sample size")
1246
1247  if (identify)
1248     { cs <- rep(c,length(n))
1249       betas <- as.vector(beta)
1250       labels <- paste("c=", formatC(cs, 2, flag="-"),
1251                       ": beta=", formatC(betas, 4, flag="-"),
1252                       ", ARL=", formatC(1/(1-betas), 2, flag="-"), sep="")
1253       i <- identify(cs, betas, labels, pos=4, offset=0.2)
1254       apply(as.matrix(labels[i$ind]), 1, cat, "\n")
1255     }
1256  else
1257     { legend(max(c), 1, legend = paste("n =", n),
1258              bg = qcc.options("bg.figure"),
1259              lty = 1:length(n), xjust = 1, yjust = 1)
1260     }
1261  invisible(beta)
1262}
1263
1264#-------------------------------------------------------------------#
1265#                                                                   #
1266# Miscellaneous functions                                           #
1267#                                                                   #
1268#-------------------------------------------------------------------#
1269
1270qcc.groups <- function(data, sample)
1271{
1272  if(length(data)!=length(sample))
1273    stop("data and sample must be vectors of equal length")
1274  x <- lapply(split(data, sample), as.vector)
1275  lx <- sapply(x, length)
1276  for(i in which(lx != max(lx)))
1277      x[[i]] <- c(x[[i]], rep(NA, max(lx)-lx[i]))
1278  x <- t(sapply(x, as.vector))
1279  return(x)
1280}
1281
1282qcc.overdispersion.test <- function(x, size,
1283                            type=ifelse(missing(size), "poisson", "binomial"))
1284{
1285  type <- match.arg(type, c("poisson", "binomial"))
1286  if (type=="binomial" & missing(size))
1287     stop("binomial data require argument \"size\"")
1288  if (!missing(size))
1289     if (length(x) != length(size))
1290        stop("arguments \"x\" and \"size\" must be vector of same length")
1291
1292  n <- length(x)
1293  obs.var <- var(x)
1294  if (type=="binomial")
1295     { p <- sum(x)/sum(size)
1296       theor.var <- mean(size)*p*(1-p) }
1297  else if (type=="poisson")
1298          { theor.var <- mean(x) }
1299       else
1300          stop("invalid \"type\" argument. See help.")
1301
1302  D <- (obs.var * (n-1)) / theor.var
1303  p.value <- 1-pchisq(D, n-1)
1304
1305  out <- matrix(c(obs.var/theor.var, D, signif(p.value,5)), 1, 3)
1306  rownames(out) <- paste(type, "data")
1307  colnames(out) <- c("Obs.Var/Theor.Var", "Statistic", "p-value")
1308  names(dimnames(out)) <- c(paste("Overdispersion test"), "")
1309  return(out)
1310}
1311
1312
1313