© 2026 Dr. Debashis Chatterjee. All rights reserved.
For educational use only


1 0) Install & load packages (run once)

pkgs <- c("qcc","ggplot2","dplyr","tidyr")
to_install <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
if(length(to_install) > 0) install.packages(to_install, dependencies = TRUE)

library(qcc)
library(ggplot2)
library(dplyr)
library(tidyr)

set.seed(123)

2 1) Practical 1 — Control charts (very simple)

2.1 1.1 Short theory (what students must remember)

A control chart checks whether a process is stable over time.

  • Center line (CL) = average level of the process
  • UCL/LCL = upper/lower control limits (usually “3-sigma limits”)
  • If a point goes beyond limits → possible special cause → investigate.

First check variability (R or S chart), then check mean ( chart).


2.2 1.2 X̄–R chart (real dataset from qcc)

We use a standard real dataset already included in the qcc package.

data(pistonrings, package="qcc")
dim(pistonrings)
## [1] 200   3
head(pistonrings)
##   diameter sample trial
## 1   74.030      1  TRUE
## 2   74.002      1  TRUE
## 3   74.019      1  TRUE
## 4   73.992      1  TRUE
## 5   74.008      1  TRUE
## 6   73.995      2  TRUE

2.2.1 Step A: R chart (variation)

qR <- qcc(pistonrings, type="R", title="R Chart (pistonrings)", plot=TRUE)

summary(qR)
## 
## Call:
## qcc(data = pistonrings, type = "R", plot = TRUE, title = "R Chart (pistonrings)")
## 
## R chart for pistonrings 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 72.96700 72.99875 73.00900 73.37860 74.00100 74.03600 
## 
## Group sample size:  3
## Number of groups:  200
## Center of group statistics:  73.3786
## Standard deviation:  43.34235 
## 
## Control limits:
##  LCL      UCL
##    0 188.8907

2.2.2 Step B: X̄ chart (mean)

qX <- qcc(pistonrings, type="xbar", title="X-bar Chart (pistonrings)", plot=TRUE)

summary(qX)
## 
## Call:
## qcc(data = pistonrings, type = "xbar", plot = TRUE, title = "X-bar Chart (pistonrings)")
## 
## xbar chart for pistonrings 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 25.33067 28.58142 31.83467 31.70953 34.75025 38.00967 
## 
## Group sample size:  3
## Number of groups:  200
## Center of group statistics:  31.70954
## Standard deviation:  43.34235 
## 
## Control limits:
##        LCL      UCL
##  -43.36162 106.7807

2.2.3 Nice simple colored plot of subgroup means (for understanding)

df <- data.frame(
  subgroup = 1:nrow(pistonrings),
  xbar = qX$statistics
)

ggplot(df, aes(subgroup, xbar)) +
  geom_line(color="#2c7fb8", linewidth=0.6) +
  geom_point(color="#2c7fb8", size=2) +
  geom_hline(yintercept=qX$limits[2], linetype="dashed", color="#d95f0e") +
  geom_hline(yintercept=qX$center,     linetype="solid",  color="#1b9e77") +
  geom_hline(yintercept=qX$limits[1], linetype="dashed", color="#d95f0e") +
  labs(title="X-bar chart (simple view)", x="Subgroup", y="Subgroup mean") +
  theme_minimal(base_size=12)


2.3 1.3 p chart (defectives) — simple real dataset (airquality)

We use airquality (real).
Define defective day = high ozone day (Ozone > 80).

data(airquality)
aq <- airquality |> filter(!is.na(Ozone)) |> mutate(defective = as.integer(Ozone > 80))
table(aq$defective)
## 
##   0   1 
## 100  16

Create 20 batches of size 10 (simple).

T <- 20
n <- 10
batches <- tibble(t=1:T) |>
  rowwise() |>
  mutate(
    samp = list(sample(aq$defective, n, replace=TRUE)),
    d = sum(unlist(samp))
  ) |> ungroup()

batches
## # A tibble: 20 × 3
##        t samp           d
##    <int> <list>     <int>
##  1     1 <int [10]>     2
##  2     2 <int [10]>     4
##  3     3 <int [10]>     2
##  4     4 <int [10]>     1
##  5     5 <int [10]>     0
##  6     6 <int [10]>     2
##  7     7 <int [10]>     2
##  8     8 <int [10]>     2
##  9     9 <int [10]>     1
## 10    10 <int [10]>     3
## 11    11 <int [10]>     0
## 12    12 <int [10]>     1
## 13    13 <int [10]>     1
## 14    14 <int [10]>     0
## 15    15 <int [10]>     1
## 16    16 <int [10]>     2
## 17    17 <int [10]>     2
## 18    18 <int [10]>     0
## 19    19 <int [10]>     2
## 20    20 <int [10]>     1

Now build a p chart (we give d and n).

qP <- qcc(batches$d, type="p", sizes=rep(n, T), title="p Chart (airquality: Ozone>80)", plot=TRUE)

summary(qP)
## 
## Call:
## qcc(data = batches$d, type = "p", sizes = rep(n, T), plot = TRUE,     title = "p Chart (airquality: Ozone>80)")
## 
## p chart for batches$d 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.100   0.150   0.145   0.200   0.400 
## 
## Group sample size:  10
## Number of groups:  20
## Center of group statistics:  0.145
## Standard deviation:  0.3521008 
## 
## Control limits:
##     LCL       UCL
##       0 0.4790322
##       0 0.4790322
## ...              
##       0 0.4790322

2.4 1.4 c chart (defects per unit) — real dataset warpbreaks

Here breaks is the count of breaks (defects).

data(warpbreaks)
qC <- qcc(warpbreaks$breaks, type="c", title="c Chart (warpbreaks)", plot=TRUE)

summary(qC)
## 
## Call:
## qcc(data = warpbreaks$breaks, type = "c", plot = TRUE, title = "c Chart (warpbreaks)")
## 
## c chart for warpbreaks$breaks 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 10.00000 18.25000 26.00000 28.14815 34.00000 70.00000 
## 
## Group sample size:  1
## Number of groups:  54
## Center of group statistics:  28.14815
## Standard deviation:  5.305483 
## 
## Control limits:
##      LCL     UCL
##  12.2317 44.0646

3 2) Practical 2 — Single sampling plan (very simple)

3.1 2.1 Short theory

A single sampling plan is \((n,c)\):

  • inspect \(n\) items,
  • accept if number of defectives \(X \le c\).

The OC curve shows \(P_a(p)=\mathbb{P}(X\le c)\) as a function of lot defective fraction \(p\).


3.2 2.2 OC curve in R (easy code)

Pa <- function(p, n, c) pbinom(c, size=n, prob=p)

n <- 50
c <- 2
p <- seq(0, 0.20, by=0.005)

oc <- data.frame(p=p, Pa=Pa(p,n,c))

ggplot(oc, aes(p, Pa)) +
  geom_line(color="#2c7fb8", linewidth=0.9) +
  geom_point(color="#2c7fb8", size=2) +
  labs(title="OC Curve (n=50, c=2)", x="p (lot fraction defective)", y="Pa(p)=P(accept)") +
  theme_minimal(base_size=12)

Check AQL and LTPD quickly:

AQL <- 0.01
LTPD <- 0.07
Pa(AQL,n,c)
## [1] 0.9861827
Pa(LTPD,n,c)
## [1] 0.3107886

4 3) Practical 3 — Process capability (simple)

4.1 3.1 Short theory

Given specs LSL and USL and data \(x\):

\[ C_p=\frac{USL-LSL}{6s},\qquad C_{pk}=\min\left(\frac{USL-\bar{x}}{3s},\frac{\bar{x}-LSL}{3s}\right) \]


4.2 3.2 Capability using pistonrings (simple demo)

x <- as.numeric(unlist(pistonrings, use.names=FALSE))
x <- x[is.finite(x)]

# For teaching: choose "spec limits" as 1% and 99% quantiles (illustration)
LSL <- as.numeric(quantile(x, 0.01, names=FALSE))
USL <- as.numeric(quantile(x, 0.99, names=FALSE))

mu <- mean(x); s <- sd(x)
Cp <- (USL-LSL)/(6*s)
Cpk <- min((USL-mu)/(3*s), (mu-LSL)/(3*s))

data.frame(LSL=LSL, USL=USL, mean=mu, sd=s, Cp=Cp, Cpk=Cpk)
##   LSL      USL     mean       sd        Cp      Cpk
## 1   0 74.02603 31.70953 31.72384 0.3889085 0.333183

Simple plot: histogram + spec lines.

df <- data.frame(x=x)
ggplot(df, aes(x)) +
  geom_histogram(bins=25, fill="#a6cee3", color="white") +
  geom_vline(xintercept=c(LSL,USL), linetype="dashed", color="#e31a1c", linewidth=1) +
  labs(title="Capability demo: histogram with (illustrative) specs", x="Measurement", y="Count") +
  theme_minimal(base_size=12)


5 4) Practical 4 — DMAIC mini-demo (very easy)

We do a simple Before vs After comparison using the same real data idea.

  • Before: original pistonrings values
  • After (demo): reduced variation + centered (only to show the method)
x_before <- x

# Simple improvement demo: reduce SD by 20% and center at midpoint of specs
x_after <- (x_before - mean(x_before))*0.80 + mean(x_before)
x_after <- x_after - (mean(x_after) - (LSL+USL)/2)

summ <- data.frame(
  stage = c("Before","After (demo)"),
  mean  = c(mean(x_before), mean(x_after)),
  sd    = c(sd(x_before), sd(x_after)),
  defect_fraction = c(mean(x_before<LSL | x_before>USL), mean(x_after<LSL | x_after>USL))
)
summ
##          stage     mean       sd defect_fraction
## 1       Before 31.70953 31.72384            0.01
## 2 After (demo) 37.01301 25.37908            0.00

Simple plot to compare distributions.

df2 <- tibble(
  value = c(x_before, x_after),
  stage = rep(c("Before","After (demo)"), each=length(x_before))
)

ggplot(df2, aes(value, fill=stage)) +
  geom_histogram(bins=30, alpha=0.6, position="identity") +
  geom_vline(xintercept=c(LSL,USL), linetype="dashed", color="#e31a1c") +
  labs(title="DMAIC demo: Before vs After distributions", x="Value", y="Count") +
  theme_minimal(base_size=12)


6 What students should write in their lab copy

For each practical, students should write:

  1. Dataset used
  2. Code and plots
  3. Control limits / key numbers
  4. Short interpretation (2–5 lines)

End of Set A (Simple).