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

1 0) Software setup

pkgs <- c("qcc","ggplot2","dplyr","tidyr","knitr")
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)
library(knitr)

set.seed(2026)

1.1 0.1 A simple table helper (safe during knitting)

k_tbl <- function(x, caption=NULL, digits=4){
  x <- as.data.frame(x)
  knitr::kable(x, caption = caption, digits = digits)
}

2 1) Practical 1 — Variables control charts (new datasets)

2.1 1.1 X̄–R chart using ChickWeight (real data)

2.1.1 What is the dataset?

ChickWeight contains body weights of chicks measured over time under different diets.

data(ChickWeight)
str(ChickWeight)
## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   578 obs. of  4 variables:
##  $ weight: num  42 51 59 64 76 93 106 125 149 171 ...
##  $ Time  : num  0 2 4 6 8 10 12 14 16 18 ...
##  $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
##  $ Diet  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "formula")=Class 'formula'  language weight ~ Time | Chick
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "outer")=Class 'formula'  language ~Diet
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "labels")=List of 2
##   ..$ x: chr "Time"
##   ..$ y: chr "Body weight"
##  - attr(*, "units")=List of 2
##   ..$ x: chr "(days)"
##   ..$ y: chr "(gm)"
summary(ChickWeight$weight)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    35.0    63.0   103.0   121.8   163.8   373.0

2.1.2 Step A: Create rational subgroups (simple approach)

For teaching, we take Diet 1 only, order observations by time, then form consecutive subgroups of size \(n=5\).

cw <- ChickWeight |> filter(Diet == 1) |> arrange(Time, Chick)
x <- as.numeric(cw$weight)
x <- x[is.finite(x)]

n <- 5
m <- floor(length(x)/n)
m_use <- min(m, 25)               # keep it small and classroom-friendly
Xmat <- matrix(x[1:(m_use*n)], nrow=m_use, ncol=n, byrow=TRUE)

dim(Xmat)
## [1] 25  5
k_tbl(head(Xmat, 5), caption="First 5 subgroups (each row is a subgroup of size 5)")
First 5 subgroups (each row is a subgroup of size 5)
V1 V2 V3 V4 V5
39 41 41 41 42
41 41 42 42 43
42 41 43 43 42
41 40 41 41 41
35 45 49 48 51

2.1.3 Step B: Visualize raw measurements (to understand variation)

ggplot(cw |> slice(1:(m_use*n)), aes(x = 1:(m_use*n), y = weight)) +
  geom_line(color="#2c7fb8", linewidth=0.5) +
  geom_point(color="#2c7fb8", size=1.6) +
  labs(title="Raw measurement stream used for subgrouping (Diet 1)",
       x="Observation index (time-ordered)", y="Weight") +
  theme_minimal(base_size = 12)

2.1.4 Step C: R chart (variation first)

qR <- qcc(Xmat, type="R", title="R Chart (ChickWeight, Diet 1)", plot=TRUE)

summary(qR)
## 
## Call:
## qcc(data = Xmat, type = "R", plot = TRUE, title = "R Chart (ChickWeight, Diet 1)")
## 
## R chart for Xmat 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    8.00   12.00   20.32   23.00   68.00 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  20.32
## Standard deviation:  8.736028 
## 
## Control limits:
##  LCL      UCL
##    0 42.96602

2.1.5 Step D: X̄ chart (mean)

qX <- qcc(Xmat, type="xbar", title="X-bar Chart (ChickWeight, Diet 1)", plot=TRUE)

summary(qX)
## 
## Call:
## qcc(data = Xmat, type = "xbar", plot = TRUE, title = "X-bar Chart (ChickWeight, Diet 1)")
## 
## xbar chart for Xmat 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   40.80   47.80   62.60   65.08   78.80  104.00 
## 
## Group sample size:  5
## Number of groups:  25
## Center of group statistics:  65.08
## Standard deviation:  8.736028 
## 
## Control limits:
##       LCL      UCL
##  53.35939 76.80061

2.1.6 Step E: Table of subgroup statistics (mean + range)

sub_mean <- rowMeans(Xmat)
sub_range <- apply(Xmat, 1, function(v) max(v) - min(v))
stats_tbl <- data.frame(
  subgroup = 1:m_use,
  xbar = sub_mean,
  R = sub_range
)
k_tbl(head(stats_tbl, 10), caption="Subgroup means and ranges (first 10 subgroups)", digits=3)
Subgroup means and ranges (first 10 subgroups)
subgroup xbar R
1 40.8 3
2 41.8 2
3 42.2 2
4 40.8 1
5 45.6 16
6 48.0 7
7 47.8 12
8 47.6 7
9 54.2 10
10 57.0 9

2.2 1.2 Individuals chart (I chart) using faithful (real data)

2.2.1 Dataset description

faithful contains waiting time between eruptions and eruption duration for Old Faithful geyser.

We monitor waiting time as a single observation stream (no subgrouping).

data(faithful)
str(faithful)
## 'data.frame':    272 obs. of  2 variables:
##  $ eruptions: num  3.6 1.8 3.33 2.28 4.53 ...
##  $ waiting  : num  79 54 74 62 85 55 88 85 51 85 ...
summary(faithful$waiting)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    43.0    58.0    76.0    70.9    82.0    96.0

2.2.2 Step A: Individuals chart

xI <- as.numeric(faithful$waiting)
qI <- qcc(xI, type="xbar.one", title="Individuals Chart (faithful: waiting time)", plot=TRUE)

summary(qI)
## 
## Call:
## qcc(data = xI, type = "xbar.one", plot = TRUE, title = "Individuals Chart (faithful: waiting time)")
## 
## xbar.one chart for xI 
## 
## Summary of group statistics:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 43.00000 58.00000 76.00000 70.89706 82.00000 96.00000 
## 
## Group sample size:  1
## Number of groups:  272
## Center of group statistics:  70.89706
## Standard deviation:  18.19175 
## 
## Control limits:
##       LCL      UCL
##  16.32181 125.4723

2.2.3 Step B: Moving Range (MR) chart (correct, size-2 subgroup method)

Moving range \(MR_t = |x_t - x_{t-1}|\) equals the range of a 2-point subgroup \((x_t, x_{t-1})\).

pair_mat <- cbind(xI[-1], xI[-length(xI)])   # (x_t, x_{t-1})
qMR <- qcc(pair_mat, type="R", title="Moving Range (MR) Chart (faithful)", plot=TRUE)

summary(qMR)
## 
## Call:
## qcc(data = pair_mat, type = "R", plot = TRUE, title = "Moving Range (MR) Chart (faithful)")
## 
## R chart for pair_mat 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  9.0000 22.0000 20.5203 31.0000 47.0000 
## 
## Group sample size:  2
## Number of groups:  271
## Center of group statistics:  20.5203
## Standard deviation:  18.19175 
## 
## Control limits:
##  LCL      UCL
##    0 67.04588

3 2) Practical 1 (continued) — Attribute control charts (new datasets)

3.1 2.1 p chart (defectives) using Titanic (real data)

3.1.1 Dataset description

Titanic is a 4-way table (Class × Sex × Age × Survived) with passenger counts. We create a binary outcome by expanding to individual records (still manageable: 2201 passengers).

data(Titanic)
Titanic_df <- as.data.frame(Titanic)

# Expand counts into individual-level rows (2201 rows)
titanic_ind <- Titanic_df[rep(1:nrow(Titanic_df), Titanic_df$Freq), c("Class","Sex","Age","Survived")]
titanic_ind$defective <- as.integer(titanic_ind$Survived == "No")  # define "defect" = death

nrow(titanic_ind)
## [1] 2201
table(titanic_ind$defective)
## 
##    0    1 
##  711 1490

3.1.2 Make simple batches (20 batches of size 50)

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

k_tbl(batches, caption="Batch defectives (Titanic: death=defective)", digits=0)
Batch defectives (Titanic: death=defective)
t samp d
1 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0 34
2 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1 36
3 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1 38
4 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0 29
5 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1 28
6 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1 38
7 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1 39
8 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0 32
9 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1 36
10 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1 36
11 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1 38
12 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0 29
13 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1 39
14 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 34
15 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1 37
16 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1 35
17 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1 32
18 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0 30
19 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1 33
20 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1 31

3.1.3 p chart

qP <- qcc(batches$d, type="p", sizes=rep(n, T), title="p Chart (Titanic: death as defective)", plot=TRUE)

summary(qP)
## 
## Call:
## qcc(data = batches$d, type = "p", sizes = rep(n, T), plot = TRUE,     title = "p Chart (Titanic: death as defective)")
## 
## p chart for batches$d 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.560   0.635   0.690   0.684   0.745   0.780 
## 
## Group sample size:  50
## Number of groups:  20
## Center of group statistics:  0.684
## Standard deviation:  0.4649129 
## 
## Control limits:
##           LCL       UCL
##     0.4867542 0.8812458
##     0.4867542 0.8812458
## ...                    
##     0.4867542 0.8812458

3.2 2.2 c chart (defects per unit) using discoveries (real data)

discoveries is a time series of counts (important discoveries per year). We treat the count as “defects per year” for chart demonstration.

data(discoveries)
xC <- as.numeric(discoveries)
summary(xC)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     2.0     3.0     3.1     4.0    12.0
qC <- qcc(xC, type="c", title="c Chart (discoveries time series counts)", plot=TRUE)

summary(qC)
## 
## Call:
## qcc(data = xC, type = "c", plot = TRUE, title = "c Chart (discoveries time series counts)")
## 
## c chart for xC 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     2.0     3.0     3.1     4.0    12.0 
## 
## Group sample size:  1
## Number of groups:  100
## Center of group statistics:  3.1
## Standard deviation:  1.760682 
## 
## Control limits:
##  LCL      UCL
##    0 8.382045

3.3 2.3 u chart (defects per opportunity) using Titanic aggregated by Class

Here: defects = deaths; opportunity = total passengers in the class.

ucb <- as.data.frame(Titanic) |> 
  group_by(Class) |> 
  summarise(
    deaths = sum(Freq[Survived=="No"]),
    total  = sum(Freq),
    u = deaths/total,
    .groups="drop"
  ) |> mutate(t = row_number())

k_tbl(ucb, caption="u chart data: deaths per passenger (by Class)", digits=4)
u chart data: deaths per passenger (by Class)
Class deaths total u t
1st 122 325 0.3754 1
2nd 167 285 0.5860 2
3rd 528 706 0.7479 3
Crew 673 885 0.7605 4 
qU <- qcc(ucb$deaths, type="u", sizes=ucb$total, title="u Chart (Titanic: deaths per passenger)", plot=TRUE)

summary(qU)
## 
## Call:
## qcc(data = ucb$deaths, type = "u", sizes = ucb$total, plot = TRUE,     title = "u Chart (Titanic: deaths per passenger)")
## 
## u chart for ucb$deaths 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.3753846 0.5333198 0.6669201 0.6174192 0.7510195 0.7604520 
## 
## Summary of group sample sizes:                         
##   sizes  285 325 706 885
##   counts   1   1   1   1
## 
## Number of groups:  4
## Center of group statistics:  0.676965
## Standard deviation:  0.8227788 
## 
## Control limits:
##        LCL       UCL
##  0.5400463 0.8138837
##  0.5307534 0.8231767
##  0.5840679 0.7698621
##  0.5939928 0.7599372

4 3) Practical 2 — Acceptance sampling (simple + more explanation)

4.1 3.1 Quick theory

A single sampling plan \((n,c)\): - Inspect \(n\) items. - Let \(X\) = number of defectives in the sample. - Accept if \(X \le c\), else reject.

Under the Binomial model \(X \sim \mathrm{Bin}(n,p)\), the OC curve is: \[ P_a(p) = \mathbb{P}(X \le c). \]

4.2 3.2 OC curve with table + plot

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

n <- 60
c <- 1
pgrid <- seq(0, 0.20, by=0.01)

oc_tbl <- data.frame(p = pgrid, Pa = Pa(pgrid, n, c))
k_tbl(oc_tbl, caption="OC table for plan (n=60, c=1)", digits=4)
OC table for plan (n=60, c=1)
p Pa
0.00 1.0000
0.01 0.8788
0.02 0.6619
0.03 0.4592
0.04 0.3022
0.05 0.1916
0.06 0.1179
0.07 0.0709
0.08 0.0418
0.09 0.0242
0.10 0.0138
0.11 0.0077
0.12 0.0043
0.13 0.0023
0.14 0.0013
0.15 0.0007
0.16 0.0004
0.17 0.0002
0.18 0.0001
0.19 0.0000
0.20 0.0000 
ggplot(oc_tbl, aes(p, Pa)) +
  geom_line(color="#2c7fb8", linewidth=0.9) +
  geom_point(color="#2c7fb8", size=2) +
  labs(title="OC Curve (n=60, c=1)", x="p (lot fraction defective)", y="Pa(p)=P(accept)") +
  theme_minimal(base_size=12)

4.3 3.3 AQL and LTPD (quick risk check)

AQL <- 0.01
LTPD <- 0.08

Pa_AQL <- Pa(AQL, n, c)
Pa_LTPD <- Pa(LTPD, n, c)

risk_tbl <- data.frame(
  quantity = c("Pa(AQL)", "Producer risk alpha = 1-Pa(AQL)", "Pa(LTPD)", "Consumer risk beta = Pa(LTPD)"),
  value = c(Pa_AQL, 1-Pa_AQL, Pa_LTPD, Pa_LTPD)
)
k_tbl(risk_tbl, caption="Risks at AQL and LTPD", digits=6)
Risks at AQL and LTPD
quantity value
Pa(AQL) 0.878767
Producer risk alpha = 1-Pa(AQL) 0.121233
Pa(LTPD) 0.041771
Consumer risk beta = Pa(LTPD) 0.041771

4.4 3.4 Real-lot demonstration using ToothGrowth (real data)

We treat the whole dataset as a “lot” (real values), define “defect” = outside illustrative specs, then repeatedly sample to estimate acceptance probability.

data(ToothGrowth)
lot <- as.numeric(ToothGrowth$len)
lot <- lot[is.finite(lot)]

# illustrative specs = 10th and 90th percentiles (teaching)
LSL <- as.numeric(quantile(lot, 0.10, names=FALSE))
USL <- as.numeric(quantile(lot, 0.90, names=FALSE))

p_true <- mean(lot < LSL | lot > USL)

n_plan <- 60
c_plan <- 1
B <- 1500

acc <- replicate(B, {
  samp <- sample(lot, n_plan, replace=TRUE)
  d <- sum(samp < LSL | samp > USL)
  as.integer(d <= c_plan)
})

k_tbl(data.frame(LSL=LSL, USL=USL, p_true=p_true, Pa_estimated=mean(acc)),
      caption="Real-lot acceptance simulation (ToothGrowth as a lot)", digits=6)
Real-lot acceptance simulation (ToothGrowth as a lot)
LSL USL p_true Pa_estimated
8.11 27.3 0.183333 0

5 4) Practical 3 — Capability (different dataset, more explanation)

5.1 4.1 Dataset: airquality Temperature (real)

We will use airquality$Temp (daily temperature in NYC, May–Sep 1973).

data(airquality)
temp <- as.numeric(airquality$Temp)
temp <- temp[is.finite(temp)]
summary(temp)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   56.00   72.00   79.00   77.88   85.00   97.00

5.1.1 Basic descriptive table

desc <- data.frame(
  N = length(temp),
  mean = mean(temp),
  sd = sd(temp),
  min = min(temp),
  Q1 = as.numeric(quantile(temp, 0.25)),
  median = median(temp),
  Q3 = as.numeric(quantile(temp, 0.75)),
  max = max(temp)
)
k_tbl(desc, caption="Descriptive summary: airquality$Temp", digits=2)
Descriptive summary: airquality$Temp
N mean sd min Q1 median Q3 max
153 77.88 9.47 56 72 79 85 97

5.2 4.2 Choose simple illustrative specs and compute Cp, Cpk

For teaching: set specs as 10th and 90th percentiles.

LSL <- as.numeric(quantile(temp, 0.10, names=FALSE))
USL <- as.numeric(quantile(temp, 0.90, names=FALSE))

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

k_tbl(data.frame(LSL=LSL, USL=USL, mean=mu, sd=s, Cp=Cp, Cpk=Cpk),
      caption="Capability indices (illustrative specs)", digits=4)
Capability indices (illustrative specs)
LSL USL mean sd Cp Cpk
64.2 90 77.8824 9.4653 0.4543 0.4267

5.2.1 Plots: histogram + QQ plot + boxplot by month

aq2 <- airquality |> filter(!is.na(Temp))

p1 <- ggplot(aq2, aes(Temp)) +
  geom_histogram(bins=20, fill="#a6cee3", color="white") +
  geom_vline(xintercept=c(LSL,USL), linetype="dashed", color="#e31a1c", linewidth=1) +
  labs(title="Temp histogram with specs (dashed)", x="Temp", y="Count") +
  theme_minimal(base_size=12)

p2 <- ggplot(aq2, aes(sample=Temp)) +
  stat_qq(color="#2c7fb8") + stat_qq_line(color="#d95f0e") +
  labs(title="Normal QQ plot (Temp)", x="Theoretical", y="Sample") +
  theme_minimal(base_size=12)

p3 <- ggplot(aq2, aes(factor(Month), Temp, fill=factor(Month))) +
  geom_boxplot(alpha=0.7) +
  guides(fill="none") +
  labs(title="Temp by Month (boxplots)", x="Month", y="Temp") +
  theme_minimal(base_size=12)

p1; p2; p3

6 5) Practical 4 — DMAIC mini-case (easy, with more explanation)

We demonstrate the workflow, not claiming real engineering intervention here.

6.1 5.1 Define

CTQ: Temperature (Temp).
Defect: Temp outside \([LSL, USL]\).
Goal: reduce defect rate and improve \(C_{pk}\).

6.2 5.2 Measure

Baseline defect fraction:

p_def_before <- mean(temp < LSL | temp > USL)
k_tbl(data.frame(defect_fraction=p_def_before, defect_ppm=1e6*p_def_before),
      caption="Baseline defect fraction (illustrative)", digits=6)
Baseline defect fraction (illustrative)
defect_fraction defect_ppm
0.196078 196078.4

6.3 5.3 Analyze (simple)

We check if average temperature differs across months (simple ANOVA-style view via boxplot already shown). Also we compute the month-wise defect rate table.

month_tbl <- aq2 |> 
  group_by(Month) |> 
  summarise(
    n = n(),
    defects = sum(Temp < LSL | Temp > USL),
    defect_rate = defects/n,
    mean_temp = mean(Temp),
    sd_temp = sd(Temp),
    .groups="drop"
  )

k_tbl(month_tbl, caption="Month-wise defect rates and summaries", digits=4)
Month-wise defect rates and summaries
Month n defects defect_rate mean_temp sd_temp
5 31 14 0.4516 65.5484 6.8549
6 30 2 0.0667 79.1000 6.5986
7 31 3 0.0968 83.9032 4.3155
8 31 5 0.1613 83.9677 6.5853
9 30 6 0.2000 76.9000 8.3557

6.4 5.4 Improve (demo)

We create an “after” dataset by reducing SD by 15% and centering to midpoint (teaching demonstration).

temp_after <- (temp - mean(temp))*0.85 + mean(temp)
temp_after <- temp_after - (mean(temp_after) - (LSL+USL)/2)

mu2 <- mean(temp_after); s2 <- sd(temp_after)
Cp2  <- (USL-LSL)/(6*s2)
Cpk2 <- min((USL-mu2)/(3*s2), (mu2-LSL)/(3*s2))
p_def_after <- mean(temp_after < LSL | temp_after > USL)

compare <- data.frame(
  stage=c("Before","After (demo)"),
  mean=c(mu, mu2),
  sd=c(s, s2),
  Cp=c(Cp, Cp2),
  Cpk=c(Cpk, Cpk2),
  defect_fraction=c(p_def_before, p_def_after)
)

k_tbl(compare, caption="Before vs After (demo) summary", digits=4)
Before vs After (demo) summary
stage mean sd Cp Cpk defect_fraction
Before 77.8824 9.4653 0.4543 0.4267 0.1961
After (demo) 77.1000 8.0455 0.5345 0.5345 0.1111

6.5 5.5 Control

Individuals chart + Moving Range chart for the “after” stream.

qI2 <- qcc(temp_after, type="xbar.one", title="Individuals Chart (After demo)", plot=TRUE)

pair_mat2 <- cbind(temp_after[-1], temp_after[-length(temp_after)])
qMR2 <- qcc(pair_mat2, type="R", title="Moving Range (MR) Chart (After demo)", plot=TRUE)

summary(qI2)
## 
## Call:
## qcc(data = temp_after, type = "xbar.one", plot = TRUE, title = "Individuals Chart (After demo)")
## 
## xbar.one chart for temp_after 
## 
## Summary of group statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   58.50   72.10   78.05   77.10   83.15   93.35 
## 
## Group sample size:  1
## Number of groups:  153
## Center of group statistics:  77.1
## Standard deviation:  3.267019 
## 
## Control limits:
##       LCL      UCL
##  67.29894 86.90106
summary(qMR2)
## 
## Call:
## qcc(data = pair_mat2, type = "R", plot = TRUE, title = "Moving Range (MR) Chart (After demo)")
## 
## R chart for pair_mat2 
## 
## Summary of group statistics:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.000000  1.700000  2.550000  3.685197  5.100000 15.300000 
## 
## Group sample size:  2
## Number of groups:  152
## Center of group statistics:  3.685197
## Standard deviation:  3.267019 
## 
## Control limits:
##  LCL      UCL
##    0 12.04063

7 6) What students should submit (short checklist)

For each practical: 1. dataset description (2–3 lines)
2. code + plot(s)
3. key numbers (limits, Cp/Cpk, Pa(AQL), etc.)
4. 3–6 lines interpretation (very important)

End of Set AA.