Control Chart and Process Capability

The code that produced every figure in this report, exactly as it ran. Fingerprint 3477eec8ea75b3a8. The same code on the same data gives the same report.

`standard_control_charts_v2` <- function(pf) {
  `%||%` <- function(a, b) if (!is.null(a)) a else b
  tidy <- function(x) {
    x <- as.numeric(x)
    ifelse(is.na(x), NA_real_,
      ifelse(abs(x) >= 1000, round(x, 0),
        ifelse(abs(x) >= 100, round(x, 1),
          ifelse(abs(x) >= 1, round(x, 2), signif(x, 3)))))
  }
  p_cell <- function(p) if (is.na(p)) NA_real_ else max(p, 1e-12)
  p_text <- function(p) if (is.na(p)) "" else if (p < 1e-4) "<0.0001" else format(signif(p, 3))

  inputs <- pf$taskList$inputs
  params <- inputs$module_parameters %||% list()
  question <- (inputs$userContext %||% list())$objective %||%
    "Is this process stable over time, and can it hold its specification?"

  col_map <- inputs$column_mapping %||% list()
  df <- renderObject.taskFunction.init(inputs, col_map)
  human <- function(sem) {
    v <- col_map[[sem]]
    if (is.null(v) || !nzchar(as.character(v))) sem else as.character(v)
  }
  n_in <- nrow(df)
  if (!("measurement" %in% names(df))) stop("column_mapping must map the measured value to measurement (one row per measurement, in the order it was taken).")

  raw_names <- local({
    ds <- inputs$dataset %||% inputs$df
    if (is.data.frame(ds)) return(names(ds))
    if (is.list(ds) && length(ds) > 0) {
      rows <- ds[seq_len(min(length(ds), 50))]
      nm <- unique(unlist(lapply(rows, function(r) if (is.list(r)) names(r) else NULL)))
      if (length(nm)) return(nm)
    }
    character(0)
  })
  mapped_actual <- unique(as.character(unlist(col_map)))
  ignored_cols  <- setdiff(raw_names, unique(c(mapped_actual, make.names(mapped_actual))))

  #' ## The measurement (95% rule)
  v <- df$measurement
  if (!is.numeric(v)) {
    ch <- trimws(as.character(v)); nb <- !is.na(ch) & ch != ""
    conv <- suppressWarnings(as.numeric(ch))
    if (sum(nb) == 0 || sum(!is.na(conv[nb])) < 0.95 * sum(nb))
      stop(sprintf("The measurement column '%s' is not numeric, so it cannot be charted.", human("measurement")))
    v <- conv
  }
  v <- as.numeric(v); v[!is.finite(v)] <- NA

  #' ## Order
  #' Every rule on this chart is about SEQUENCE, so the order is stated rather than assumed.
  has_order <- "order_by" %in% names(df)
  ord_note <- "the order the rows arrived in"
  idx <- seq_len(n_in)
  if (has_order) {
    o <- df$order_by
    op <- suppressWarnings(as.POSIXct(as.character(o), tz = "UTC"))
    key <- if (sum(!is.na(op)) >= 0.95 * sum(!is.na(o))) as.numeric(op) else suppressWarnings(as.numeric(as.character(o)))
    if (all(is.na(key))) { key <- seq_len(n_in); ord_note <- sprintf("the order the rows arrived in ('%s' could not be read as a date or a number)", human("order_by")) }
    else ord_note <- sprintf("'%s'", human("order_by"))
    idx <- order(key, na.last = TRUE)
  }
  v <- v[idx]
  sub_all <- if ("subgroup" %in% names(df)) as.character(df$subgroup)[idx] else NULL
  keep <- !is.na(v)
  n_missing <- sum(!keep)
  x <- v[keep]
  sub <- if (!is.null(sub_all)) sub_all[keep] else NULL
  n <- length(x)
  if (n < 20) stop(sprintf("Only %d usable measurements; a control chart needs at least 20 to estimate its own limits.", n))
  if (isTRUE(stats::sd(x) == 0)) stop("Every measurement is the same value, so the process has no variation to chart.")

  #' ## Chart geometry: sigma from the average moving range (d2 = 1.128 for consecutive pairs)
  cl <- mean(x)
  mr <- abs(diff(x))
  mrbar <- mean(mr)
  sigma <- mrbar / 1.128
  if (!is.finite(sigma) || sigma <= 0) stop("The measurements never change between consecutive readings, so no control limits can be estimated.")
  ucl <- cl + 3 * sigma; lcl <- cl - 3 * sigma
  two_up <- cl + 2 * sigma; two_dn <- cl - 2 * sigma
  mr_ucl <- 3.267 * mrbar   # D4 for consecutive pairs

  same_side_run <- function(vv, centre, k) {
    side <- sign(vv - centre); flag <- rep(FALSE, length(vv)); run <- 0; last <- 0
    for (i in seq_along(vv)) {
      if (side[i] == 0) { run <- 0; last <- 0; next }
      if (side[i] == last) run <- run + 1 else { run <- 1; last <- side[i] }
      if (run >= k) flag[(i - k + 1):i] <- TRUE
    }
    flag
  }
  trend_run <- function(vv, k) {
    dv <- sign(diff(vv)); flag <- rep(FALSE, length(vv)); run <- 0; last <- 0
    for (i in seq_along(dv)) {
      if (dv[i] == 0) { run <- 0; last <- 0; next }
      if (dv[i] == last) run <- run + 1 else { run <- 1; last <- dv[i] }
      if (run >= k - 1) flag[(i - k + 2):(i + 1)] <- TRUE
    }
    flag
  }
  two_of_three <- function(vv, up, dn) {
    flag <- rep(FALSE, length(vv))
    for (i in seq_len(max(0, length(vv) - 2))) {
      w <- i:(i + 2); hi <- vv[w] > up; lo <- vv[w] < dn
      if (sum(hi) >= 2) flag[w[hi]] <- TRUE
      if (sum(lo) >= 2) flag[w[lo]] <- TRUE
    }
    flag
  }
  f1 <- (x > ucl) | (x < lcl)
  f2 <- same_side_run(x, cl, 9)
  f3 <- trend_run(x, 6)
  f5 <- two_of_three(x, two_up, two_dn)
  any_flag <- f1 | f2 | f3 | f5
  n_flagged <- sum(any_flag)
  in_control <- n_flagged == 0

  fmt_points <- function(w) if (!length(w)) "none" else {
    if (length(w) <= 8) paste(w, collapse = ", ") else paste0(paste(utils::head(w, 8), collapse = ", "), " and ", length(w) - 8, " more")
  }
  rule_defs <- list(
    list(rule = "A point beyond the limits", description = "A single measurement outside the 3-sigma control limits", flag = f1),
    list(rule = "Nine in a row on one side", description = "Nine consecutive measurements on the same side of the process mean, which is a shift rather than nine separate problems", flag = f2),
    list(rule = "Six steadily rising or falling", description = "Six consecutive measurements moving in one direction", flag = f3),
    list(rule = "Two of three beyond 2 sigma", description = "Two measurements out of three consecutive ones past the same 2-sigma line", flag = f5))
  rules_df <- do.call(rbind, lapply(rule_defs, function(r) data.frame(
    rule = r$rule, description = r$description,
    points_flagged = fmt_points(which(r$flag)), count = sum(r$flag), stringsAsFactors = FALSE)))

  #' ## Specification limits: a parameter wins; otherwise a mapped column holding ONE value.
  resolve_spec <- function(pname, sem) {
    pv <- params[[pname]]
    if (!is.null(pv) && nzchar(as.character(pv))) {
      num <- suppressWarnings(as.numeric(pv))
      if (is.na(num)) stop(sprintf("module_parameters$%s must be a number.", pname))
      return(list(value = num, source = sprintf("the %s parameter", pname)))
    }
    if (sem %in% names(df)) {
      cv <- suppressWarnings(as.numeric(as.character(df[[sem]])))
      u <- unique(cv[!is.na(cv)])
      if (length(u) == 1) return(list(value = u, source = sprintf("the column '%s'", human(sem))))
      if (length(u) > 1) stop(sprintf("The column '%s' holds %d different specification limits; a capability study needs one limit for the whole series. Filter to one part or pass the limit as a parameter.", human(sem), length(u)))
    }
    list(value = NA_real_, source = "")
  }
  lo <- resolve_spec("lower_spec", "lower_spec")
  hi <- resolve_spec("upper_spec", "upper_spec")
  lsl <- lo$value; usl <- hi$value
  has_spec <- is.finite(lsl) || is.finite(usl)

  #' ## Spread: short-term (within) and overall
  sd_overall <- stats::sd(x)
  sigma_within <- sigma
  sigma_source <- "the average moving range between consecutive measurements divided by 1.128"
  if (!is.null(sub)) {
    tb <- table(sub)
    good <- names(tb)[tb >= 2]
    if (length(good) >= 2) {
      ss <- 0; dfree <- 0
      for (g in good) { gv <- x[sub == g]; ss <- ss + sum((gv - mean(gv))^2); dfree <- dfree + length(gv) - 1 }
      if (dfree > 0 && ss > 0) {
        sigma_within <- sqrt(ss / dfree)
        sigma_source <- sprintf("the pooled spread within the %d subgroups of '%s'", length(good), human("subgroup"))
      }
    }
  }
  sw <- if (n >= 3 && n <= 5000) stats::shapiro.test(x) else NULL
  normal_ok <- is.null(sw) || sw$p.value >= 0.05

  cap_df <- NULL; dist_df <- NULL; cpk <- NA_real_
  if (has_spec) {
    cpu <- if (is.finite(usl)) (usl - cl) / (3 * sigma_within) else NA_real_
    cpl <- if (is.finite(lsl)) (cl - lsl) / (3 * sigma_within) else NA_real_
    ppu <- if (is.finite(usl)) (usl - cl) / (3 * sd_overall) else NA_real_
    ppl <- if (is.finite(lsl)) (cl - lsl) / (3 * sd_overall) else NA_real_
    cpk <- suppressWarnings(min(c(cpu, cpl), na.rm = TRUE))
    ppk <- suppressWarnings(min(c(ppu, ppl), na.rm = TRUE))
    cp  <- if (is.finite(lsl) && is.finite(usl)) (usl - lsl) / (6 * sigma_within) else NA_real_
    ppv <- if (is.finite(lsl) && is.finite(usl)) (usl - lsl) / (6 * sd_overall) else NA_real_
    idx_rows <- list(
      c("Cp (short-term width)", cp), c("Cpk (short-term, centred)", cpk),
      c("Pp (overall width)", ppv),   c("Ppk (overall, centred)", ppk))
    keep_idx <- vapply(idx_rows, function(r) is.finite(as.numeric(r[2])), logical(1))
    cap_df <- data.frame(
      index = vapply(idx_rows[keep_idx], function(r) r[1], character(1)),
      value = tidy(vapply(idx_rows[keep_idx], function(r) as.numeric(r[2]), numeric(1))),
      stringsAsFactors = FALSE)
    dist_df <- data.frame(value = tidy(x), series = rep("Measurements", n), stringsAsFactors = FALSE)
  }

  #' ## Places
  chart_df <- rbind(
    data.frame(point = seq_len(n), value = tidy(x),  series = "Measurement", stringsAsFactors = FALSE),
    data.frame(point = seq_len(n), value = tidy(cl), series = "Process mean", stringsAsFactors = FALSE),
    data.frame(point = seq_len(n), value = tidy(ucl), series = "Upper control limit", stringsAsFactors = FALSE),
    data.frame(point = seq_len(n), value = tidy(lcl), series = "Lower control limit", stringsAsFactors = FALSE))
  mr_df <- rbind(
    data.frame(point = 2:n, moving_range = tidy(mr), series = "Moving range", stringsAsFactors = FALSE),
    data.frame(point = 2:n, moving_range = tidy(mr_ucl), series = "Upper limit", stringsAsFactors = FALSE))

  stats_rows <- list(
    c("Measurements charted", format(n)),
    c("Process mean (centre line)", format(tidy(cl))),
    c("Sigma (short-term)", format(tidy(sigma_within))),
    c("Standard deviation (overall)", format(tidy(sd_overall))),
    c("Upper control limit", format(tidy(ucl))),
    c("Lower control limit", format(tidy(lcl))),
    c("Average moving range", format(tidy(mrbar))),
    c("Measurements a rule flagged", sprintf("%d of %d", n_flagged, n)))
  if (has_spec) {
    if (is.finite(lsl)) stats_rows[[length(stats_rows) + 1]] <- c("Lower specification limit", format(tidy(lsl)))
    if (is.finite(usl)) stats_rows[[length(stats_rows) + 1]] <- c("Upper specification limit", format(tidy(usl)))
    n_out <- sum((is.finite(lsl) & x < lsl) | (is.finite(usl) & x > usl))
    stats_rows[[length(stats_rows) + 1]] <- c("Measurements outside specification", sprintf("%d of %d", n_out, n))
  }
  process_df <- data.frame(
    statistic = vapply(stats_rows, function(r) r[1], character(1)),
    value     = vapply(stats_rows, function(r) r[2], character(1)),
    stringsAsFactors = FALSE)

  checks_list <- list(
    list(check = "Enough measurements for limits",
         statistic = sprintf("%d measurements", n), p = NA_real_,
         verdict = if (n >= 25) "holds" else "strained",
         note = "Limits estimated from fewer than about 25 measurements move a lot as more arrive."),
    list(check = "The measurements are in a real order",
         statistic = sprintf("ordered by %s", ord_note), p = NA_real_,
         verdict = if (has_order) "holds" else "strained",
         note = "Every rule here is about sequence. Without a column saying what the order is, the chart trusts the order the rows arrived in."),
    list(check = "The process was stable over the series",
         statistic = sprintf("%d of %d measurements flagged by a rule", n_flagged, n), p = NA_real_,
         verdict = if (in_control) "holds" else "strained",
         note = "A short-term index describes a process that holds still. When rules have fired, the short-term figures describe a process that was not one process over this period."))
  if (!is.null(sw)) checks_list[[length(checks_list) + 1]] <- list(
    check = "Measurements are normally distributed",
    statistic = sprintf("Shapiro-Wilk W %s", format(tidy(sw$statistic))), p = p_cell(sw$p.value),
    verdict = if (normal_ok) "holds" else "strained",
    note = if (has_spec)
      "The capability indices assume a normal distribution and overstate or understate the defect rate badly when it fails. Read them with this row." else
      "Noted for the capability indices; no specification limit was supplied, so none were computed.")
  checks_df <- data.frame(
    check     = vapply(checks_list, function(c) c$check, character(1)),
    statistic = vapply(checks_list, function(c) c$statistic, character(1)),
    p_value   = vapply(checks_list, function(c) c$p, numeric(1)),
    verdict   = vapply(checks_list, function(c) c$verdict, character(1)),
    note      = vapply(checks_list, function(c) c$note, character(1)),
    stringsAsFactors = FALSE)

  spec_text <- if (has_spec)
    sprintf(" Specification limits came from %s.", paste(Filter(nzchar, c(lo$source, hi$source)), collapse = " and ")) else
    " No specification limit was supplied, so no capability index was computed."
  method <- sprintf(
    "Individuals and moving-range control chart of %d measurements of '%s' (%d supplied%s), ordered by %s. The centre line is the process mean and the short-term sigma is %s, giving 3-sigma natural process limits. Four special-cause rules were run: a point beyond the limits, nine in a row on one side of the centre line, six steadily rising or falling, and two of three beyond the same 2-sigma line.%s%s The control limits describe what this process does and are not a specification; a stable process is not thereby a good one.",
    n, human("measurement"), n_in,
    if (n_missing > 0) sprintf(", %d without a measurement", n_missing) else "",
    ord_note, sigma_source, spec_text,
    if (has_spec) " Cp and Cpk use the short-term sigma, Pp and Ppk the overall standard deviation across the whole series; the gap between the two pairs is the drift or shift over the period." else "")

  assumptions <- list(
    "Control limits come from the process's own variation, not from anyone's specification. A process can be perfectly stable and still fail its specification, and an unstable one can sit well inside it.",
    "Every rule here is about the order of the measurements. If the rows are not in the order they were taken, the rules are meaningless.",
    "A rule that fires identifies a pattern, not a set of individually faulty measurements.",
    "The capability indices assume the measurements are normally distributed, and mislead when they are not.",
    "A short-term index describes a process that is holding still; when the chart shows the process moved, it describes a state the process was not in.")

  answer <- if (in_control)
    sprintf("No special-cause rule fired across %d measurements: the process looks stable, with a mean of %s and limits of %s to %s.",
            n, format(tidy(cl)), format(tidy(lcl)), format(tidy(ucl)))
  else
    sprintf("%d of %d measurements were flagged by a special-cause rule, so the process was not stable over this period; the mean is %s and the limits are %s to %s.",
            n_flagged, n, format(tidy(cl)), format(tidy(lcl)), format(tidy(ucl)))

  results <- list()
  summary_vals <- list(points = n, centre = tidy(cl), sigma = tidy(sigma_within),
                       ucl = tidy(ucl), lcl = tidy(lcl), flagged = n_flagged)
  if (has_spec && is.finite(cpk)) summary_vals$cpk <- tidy(cpk)
  results$summary_metrics <- place_metric(summary_vals, lead = "points", place = "summary_metrics")
  results$control_chart <- place_trend(chart_df, x = "point", y = "value", series = "series", place = "control_chart")
  results$moving_range  <- place_trend(mr_df, x = "point", y = "moving_range", series = "series", place = "moving_range")
  if (has_spec) {
    results$capability_distribution <- place_distribution(dist_df, x = "value", series = "series", place = "capability_distribution")
    results$capability_indices <- place_comparison(cap_df, category = "index", value = "value", place = "capability_indices")
  } else {
    results$capability_distribution <- place_dropped("no specification limit was supplied, so there is nothing to place the measurements against; the control chart above still shows what the process does", place = "capability_distribution")
    results$capability_indices <- place_dropped("no specification limit was supplied, so no capability index was computed; supply a lower_spec, an upper_spec or both and these appear", place = "capability_indices")
  }
  results$process_table     <- place_table(process_df, place = "process_table")
  results$rule_violations   <- place_table(rules_df, place = "rule_violations")
  results$assumption_checks <- place_table(checks_df, place = "assumption_checks")
  results$spc_method <- list(kind = "metric", values = list(
      method = method, n_in = n_in, n_used = n,
      excluded = as.list(unname(ignored_cols)),
      assumptions = assumptions,
      x_column = sprintf("%s in order", human("measurement")),
      y_column = if (has_spec) sprintf("%d rules and the capability indices", length(rule_defs)) else sprintf("%d special-cause rules", length(rule_defs))),
    value_order = list("n_used", "n_in"))

  objects <- list()
  list(answer = answer, method = method, n = n, results = results, objects = objects,
       json_output = list(answer = answer, method = method, n = n))
}