--- title: "Monitoring univariate dispersion with R and S charts" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Monitoring univariate dispersion with R and S charts} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) library(IQCC) ``` ## Choosing a dispersion statistic An R chart monitors each subgroup through its sample range \[ R = X_{(n)} - X_{(1)}, \] whereas an S chart monitors its sample standard deviation. The range is easy to obtain and remains useful for the small, fixed subgroup sizes common in routine process control. The standard deviation uses all observations and is generally preferable as subgroup size increases. Both charts assume independent subgroups from a stable process; the exact probability calculations below also assume normal observations within each subgroup. The range and standard deviation are nonnegative and have skewed sampling distributions at small sample sizes. A symmetric normal approximation can therefore misplace one or both control limits. IQCC exposes conventional and distribution-based versions through `cchart.R()` and `cchart.S()`. ## The relative range For a normal process with standard deviation \(\sigma\), Barbosa, Gneri, and Meneguetti (2013) define the relative range \[ W = \frac{R}{\sigma}. \] Its distribution function can be evaluated with the Tukey studentized-range distribution at infinite denominator degrees of freedom: \[ F_W(w;n) = F_{\mathrm{Tukey}}(w;n,\infty). \] Consequently, `stats::ptukey(w, nmeans = n, df = Inf)` evaluates \(F_W\), and `stats::qtukey(p, nmeans = n, df = Inf)` obtains its \(p\) quantile. The IQCC functions `d2(n)` and `d3(n)` calculate the mean and standard deviation of \(W\), respectively. ```{r relative-range-moments} sizes <- c(2, 5, 10, 15) data.frame( n = sizes, d2 = d2(sizes), d3 = d3(sizes) ) ``` ## Classical and exact R limits The conventional three-sigma limits are \[ LCL_R = \max\{0,d_2(n)-3d_3(n)\}\hat\sigma, \qquad UCL_R = \{d_2(n)+3d_3(n)\}\hat\sigma. \] For nominal two-sided false-alarm probability \(\alpha\), equal-tail exact probability limits are \[ LCL_R = F_W^{-1}(\alpha/2;n)\hat\sigma, \qquad UCL_R = F_W^{-1}(1-\alpha/2;n)\hat\sigma. \] `cchart.R(type = "tukey")` uses \(\alpha=0.0027\), so its probabilities are 0.00135 and 0.99865. The scale estimate is obtained from the Phase I reference subgroups supplied as `y`; the subgroups in `x` are then plotted as Phase II. The chart signals below the lower limit or above the upper limit. ## Reproducing Table 2 Table 2 of Barbosa et al. (2013) prints \(W\) quantiles to five decimal places. The following calculation reproduces small, intermediate, and large rows directly from the stated Tukey relationship. ```{r published-table-2} table2 <- data.frame( n = c(2, 5, 10), published_q_001 = c(0.00177, 0.36739, 1.08458), published_q_99865 = c(4.53274, 5.37740, 5.87416) ) table2$calculated_q_001 <- stats::qtukey( 0.001, nmeans = table2$n, df = Inf ) table2$calculated_q_99865 <- stats::qtukey( 0.99865, nmeans = table2$n, df = Inf ) table2$error_q_001 <- table2$calculated_q_001 - table2$published_q_001 table2$error_q_99865 <- table2$calculated_q_99865 - table2$published_q_99865 table2 ``` All discrepancies are below \(5\times10^{-6}\), half a unit at the fifth decimal place. The package's tabular helper produces the same quantiles for a chosen nominal alpha. ```{r tukey-table-helper} table.qtukey(alpha = 0.0027, n = 5) ``` ## False-alarm risk and subgroup size For the classical limits, the actual in-control risk is \[ 1 - \left[F_W\{d_2(n)+3d_3(n)\} - F_W\{\max(0,d_2(n)-3d_3(n))\}\right]. \] `alpha.risk()` performs this exact evaluation. It does not report a nominal design value; it reports the probability induced by the conventional limits. ```{r false-alarm-risk} risk <- alpha.risk(sizes) data.frame( n = sizes, nominal_alpha = 0.0027, actual_alpha = risk, inflation = risk / 0.0027, in_control_arl = 1 / risk ) ``` The inflation is most severe for the smallest subgroups and remains material over this range of \(n\). Exact probability limits instead allocate the chosen tail probabilities by construction, conditional on the process scale used to form them. ## Phase I and Phase II with an R chart Phase I reference data establish the process scale and the fixed limits. Phase II data are prospective subgroups evaluated against that design. The exact IQCC wrapper makes this split explicit through `y` and `x` for the limits. ```{r exact-r-chart} data(pistonrings) phase1 <- pistonrings[1:25, ] phase2 <- pistonrings[26:40, ] r_chart <- cchart.R( x = phase2, n = 5, type = "tukey", y = phase1 ) r_chart$limits r_chart$statistics ``` The returned `qcc` object's center line is still the mean range calculated from `x`, not the Phase I expected range based on `y`. Thus the current wrapper freezes the Phase I scale for its limits but does not freeze every displayed design quantity. This distinction does not change the limit-crossing signals, but it matters when interpreting the center line. The conventional wrapper delegates scale estimation and limits to `qcc` using the supplied `x`. It is appropriate for a retrospective chart of a stable reference data set, but its interface does not accept a separate Phase I scale for prospective Phase II monitoring. ```{r conventional-r-chart} r_chart_classical <- cchart.R( x = phase1, n = 5, type = "norm" ) r_chart_classical$limits ``` ## S charts For normal samples, \[ \frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1}. \] This gives equal-tail probability limits \[ \hat\sigma\sqrt{\frac{\chi^2_{\alpha/2,n-1}}{n-1}} \quad\text{and}\quad \hat\sigma\sqrt{\frac{\chi^2_{1-\alpha/2,n-1}}{n-1}}. \] `cchart.S(type = "n")` requests the normalized `qcc` chart, while `cchart.S(type = "e", m = n)` uses the chi-square formula with \(\alpha=0.0027\). In the current interface, both variants estimate the scale from `x`; there is no separate Phase I/Phase II argument. For comparison, the three-sigma construction underlying a classical S chart uses \(E(S)=c_4(n)\sigma\) and \(SD(S)=\sigma\sqrt{1-c_4(n)^2}\). With a process-scale estimate \(\hat\sigma\), its limits have the form \[ \max\{0,c_4(n)-3\sqrt{1-c_4(n)^2}\}\hat\sigma \quad\text{and}\quad \{c_4(n)+3\sqrt{1-c_4(n)^2}\}\hat\sigma. \] As with R, these moment-based limits are not equal-tail probability limits for small \(n\). ```{r s-charts} data(softdrink) s_chart_normalized <- cchart.S(softdrink, type = "n") s_chart_exact <- cchart.S(softdrink, type = "e", m = 10) rbind( normalized = s_chart_normalized$limits, exact = s_chart_exact$limits ) ``` ## Interpretation and limitations A point above the upper limit indicates increased within-subgroup dispersion; a point below a positive lower limit indicates decreased dispersion. A signal is evidence against the in-control model, not a diagnosis of its cause. Check data quality, rational subgrouping, independence, and process context before acting on it. The exact calibration relies on normal within-subgroup observations and a fixed known scale. Operational charts replace that scale by a Phase I estimate. Their limits are therefore plug-in limits and do not account for additional uncertainty from estimating \(\sigma\), especially when the number of reference subgroups is small. The current R and S graphical wrappers also do not expose all displayed limits through dedicated pure numerical limit functions; use the returned `qcc` object's `limits` component when auditing a fitted chart. ## References Barbosa, E. P., Gneri, M. A., and Meneguetti, A. (2013). Range control charts revisited: Simpler Tippett-like formulae, its practical implementation, and the study of false alarm. *Communications in Statistics - Simulation and Computation*, 42(2), 247--262. doi: [10.1080/03610918.2011.639967](https://doi.org/10.1080/03610918.2011.639967). Montgomery, D. C. (2009). *Introduction to Statistical Quality Control*, 6th ed. Wiley. ```{r execution-marker, results='asis'} cat("\n") ```