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