wer_all.R

library(plyr)
library(qcc)

data("pistonrings")
attach(pistonrings)
x <- qcc.groups(diameter, sample)
detach(pistonrings)

find_zones <- function(x) {

  tmp <- qcc(x, type="xbar", plot = FALSE)
  xbars <- tmp$statistics
  x.mean <- mean(xbars)
  # x.mean <- mean(x)
  # x.sd <- sd(x)
  x.sd <- diff(as.numeric(tmp$limits))/6
  
  boundaries <- seq(-4, 4)
  # creates a set of zones for each point in x
  zones <- sapply(boundaries, function(i) {
    i * rep(x.sd, length(xbars))
  })
  zones + x.mean
}

head(find_zones(x))
# [,1]     [,2]    [,3]    [,4]    [,5]     [,6]     [,7]
# [1,] 73.98559 73.99009 73.9946 73.9991 74.0036 74.00811 74.01261
# [2,] 73.98559 73.99009 73.9946 73.9991 74.0036 74.00811 74.01261
# [3,] 73.98559 73.99009 73.9946 73.9991 74.0036 74.00811 74.01261
# [4,] 73.98559 73.99009 73.9946 73.9991 74.0036 74.00811 74.01261
# [5,] 73.98559 73.99009 73.9946 73.9991 74.0036 74.00811 74.01261
# [6,] 73.98559 73.99009 73.9946 73.9991 74.0036 74.00811 74.01261
# [,8]     [,9]
# [1,] 74.01712 74.02162
# [2,] 74.01712 74.02162
# [3,] 74.01712 74.02162
# [4,] 74.01712 74.02162
# [5,] 74.01712 74.02162
# [6,] 74.01712 74.02162


evaluate_zones <- function(x) {
  zones <- find_zones(x)
  colnames(zones) <- paste("zone", -4:4, sep="_")

  tmp <- qcc(x, type="xbar", plot = FALSE)
  xbars <- tmp$statistics
  x.zones <- rowSums(xbars > zones) - 4
  x.zones
}

evaluate_zones(x)
# [1]  2  0  1  0  0 -1  0 -1  1 -1 -2  0 -1 -2  1 -1  0  1 -1  2  0
# [22]  0  0  1 -1  2  0 -2  0 -1  1  1 -1  2  2  1  3  4  5  3

find_violations <- function(x.zones, i) {
  rule1=0; rule2=0; rule3=0; rule4=0
  if (i > 7){
  values <- x.zones[max(i-7, 1):i]
  rule4 <- ifelse(all(values > 0), 1,
                  ifelse(all(values < 1), -1,
                         0))
  }
  if (i > 4){
  values <- x.zones[max(i-4, 1):i]
  rule3 <- ifelse(sum(values > 1) >= 4, 1,
                  ifelse(sum(values < 0) >= 4, -1,
                         0))
  }
  if (i > 2){
  values <- x.zones[max(i-2, 1):i]
  rule2 <- ifelse(sum(values > 2) >= 2, 1,
                  ifelse(sum(values < -1) >= 2, -1,
                         0))
  }
  values <- x.zones[i]
  rule1 <- ifelse(any(values > 3), 1,
                  ifelse(any(values < -2), -1,
                         0))
  
  c("rule1"=rule1, "rule2"=rule2, "rule3"=rule3, "rule4"=rule4)
}

find_violations(evaluate_zones(x), 38)
# rule1 rule2 rule3 rule4 
# 1     1     1     0 

compute_violations <- function(x, start=1) {
  x.zones <- evaluate_zones(x)
  tmp <- qcc(x, type="xbar", plot = FALSE)
  xbars <- tmp$statistics
  results <- ldply(start:length(xbars), function(i) {
    find_violations(x.zones, i)
  })
  results$color <- ifelse(results$rule1!=0, "pink",
                          ifelse(results$rule2!=0, "red",
                                 ifelse(results$rule3!=0, "orange",
                                        ifelse(results$rule4!=0, "yellow",
                                               "black"))))
  results
}

tail(compute_violations(x))
# rule1 rule2 rule3 rule4 color
# 35     0     0     0     0 black
# 36     0     0     0     0 black
# 37     0     0     0     0 black
# 38     1     1     1     0  pink
# 39     1     1     1     0  pink
# 40     0     1     1     0   red

plot.wer <- function(x, holdout = dim(x)[1] - 1) {
  
  tmp <- qcc(x, type="xbar", plot = FALSE)
  xbars <- tmp$statistics
  
  wer <- compute_violations(x, start = length(xbars) - holdout)
  bands <- find_zones(x)
  plot.data <- xbars[(length(xbars) - holdout):length(xbars)]
  plot(plot.data, col=wer$color, type='b', pch=19,
       ylim=c(min(bands), max(bands)),
       main="Western Electric Rule Ouput",
       xlab="", ylab="")
  cols = c("green","red", "orange", "yellow", "black", "yellow", "orange", "red","green")
  
  for (i in 2:8) {
    lines(bands[,i], col=cols[i], lwd=0.75, lty=2)
  }
}

plot.wer(x)