How can I plot multiple control limits in an EWMA control chart with qcc?

Question

I can plot an EWMA control chart with 3 sigma limits (nsigmas=3). Can anyone help me with plotting additional control lines on the same chart such as 1 and 2 sigma limits?

The only way I can think is to create separate qcc objects with each of these limits, somehow extract their values, and then plot them onto the EWMA chart. Surely there's a simpler way?

library(qcc)
LRR <- c(-0.1, -0.1, -0.09, -0.07, -0.27, -0.18,
-0.8, -0.86, -0.82,  0.01,  0.02)

q1 <- ewma(LRR, center =  -0.3, std.dev = 0.1 , lambda = 0.2,plot=F)

plot(q1, add.stats = F, label.limits = c("LCL", "UCL"), 
     xlab="Group", ylab= "LRR", ylim=c(-1.0,0.5),nsigmas = 3)

Now I'd just like to add 1 & 2 sigma control limits.

Thanks.

Answer 1

Thank you so much Paul. Your answer is certainly more sophisticated than mine was!

I ended up solving it by;

'escalc' to find the LRR, then

    rma(yi=LRR,vi=LRR_var)

est = 'estimate' from this rma

    q1 <- ewma(LRR, center =  est,  lambda = 0.2,labels=year,plot=F)
    LCL1 <- est - q1$sigma
    LCL2 <- est - 2*q1$sigma
    UCL1 <- est + q1$sigma
    UCL2 <- est + 2*q1$sigma
    LCL3 <- est - 3*q1$sigma
    UCL3 <- est + 3*q1$sigma

Then using plot

    abline(h=center,lty=1)
    lines(spline(1:4,LCL1,n=10), lty = "dotted", col = "grey")
    lines(spline(1:4,UCL1,n=10), lty = "dotted", col = "grey")
    lines(spline(1:4,LCL2,n=10), lty = "dotted", col = "black")
    lines(spline(1:4,UCL2,n=10), lty = "dotted", col = "black")
    lines(spline(1:4,LCL3,n=10), lty = "longdash", col = "black")
    lines(spline(1:4,UCL3,n=10), lty = "longdash", col = "black")

Answer 2

One way to go is to adjust the original plot.ewma.qcc function which you can find here. With a bit of effort you can work out how the limits are handled in this function. You could create your own plot.ewma.qcc2 function as shown below. The function is lengthy but halfway through, below the comment add more control limits you will see I added four lines() function calls to calculate and display the 1 and 2 sigma control lines.

Here's how this looks like:

plot.ewma.qcc2 <- 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
  if ((missing(object)) | (!inherits(object, "ewma.qcc"))) 
    stop("an object of class `ewma.qcc' is required")
  type <- object$type
  data.name <- object$data.name
  center <- object$center
  std.dev <- object$std.dev
  stats <- object$statistics
  limits <- object$limits
  newstats <- object$newstats
  newdata.name <- object$newdata.name
  violations <- object$violations
  nviolations <- length(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("EWMA Chart\nfor", data.name)
    else if (chart.all) 
      main.title <- paste("EWMA Chart\nfor", data.name, 
                          "and", newdata.name)
    else main.title <- paste("EWMA 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(indices, statistics, type = "n", ylim = if (!missing(ylim)) 
    ylim
    else range(statistics, limits), ylab = ifelse(missing(ylab), 
                                                  "Group Summary Statistics", ylab), xlab = ifelse(missing(xlab), 
                                                                                                   "Group", 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"))
  abline(h = center, lty = 1)
  lines(indices, limits[indices, 1], lty = 2)
  lines(indices, limits[indices, 2], lty = 2)
  
# add more control limits
  
  lines(indices, center + object$sigma, lty = 2, col = "red")
  lines(indices, center - object$sigma, lty = 2, col = "red")
  
  lines(indices, center + 2*object$sigma, lty = 2, col = "green")
  lines(indices, center - 2*object$sigma, lty = 2, col = "green")
  
  points(indices, statistics, pch = 3, cex = 0.8)
  lines(indices, object$y[indices], type = "o", pch = 20)
  mtext(label.limits, side = 4, at = limits[nrow(limits), ], 
        las = 1, line = 0.1, col = gray(0.3), cex = par("cex"))
  if (nviolations > 0) {
    if (is.null(qcc.options("beyond.limits"))) 
      stop(".qcc.options$beyond.limits undefined. See help(qcc.options).")
    points(violations, object$y[violations], col = qcc.options("beyond.limits")$col, 
           pch = qcc.options("beyond.limits")$pch)
  }
  if (chart.all & !is.null(newstats)) {
    len.obj.stats <- length(stats)
    len.new.stats <- length(newstats)
    abline(v = len.obj.stats + 0.5, lty = 3)
    mtext("Calibration data", cex = par("cex") * 
            0.8, at = len.obj.stats/2, line = 0, adj = 0.5)
    mtext("New data", cex = par("cex") * 0.8, 
          at = len.obj.stats + len.new.stats/2, line = 0, adj = 0.5)
  }
  if (add.stats) {
    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.15, 0.6)
    top.line <- 4.5
    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"))
    if (length(center) == 1) {
      mtext(paste("Center = ", signif(center[1], 
                                      digits = 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"))
    }
    mtext(paste("StdDev = ", signif(std.dev, digits = 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"))
    mtext(paste("Smoothing parameter = ", signif(object$lambda, 
                                                 digits = digits)), side = 1, line = top.line, adj = 0, 
          at = at.col[2], font = qcc.options("font.stats"), 
          cex = par("cex") * qcc.options("cex.stats"))
    mtext(paste("Control limits at ", object$nsigmas, 
                "*sigma", 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"))
    mtext(paste("No. of points beyond limits =", nviolations), 
          side = 1, line = top.line + 2, adj = 0, at = at.col[2], 
          font = qcc.options("font.stats"), cex = par("cex") * 
            qcc.options("cex.stats"))
  }
  invisible()
}