<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/estimating-control-chart-constants-with-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://r-bar.net/estimating-control-chart-constants-with-r/"> R – R-BAR</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In this post, I will show you how a very basic R code can be used to estimate quality control constants needed to construct X-Individuals, X-Bar, and R-Bar charts. The value of this approach is that it gives you a mechanical sense of where these constants come from and some reinforcement on their application.</p>
<p>If you work in a production or quality control environment, chances are you've made or seen a control chart. If you're new to control charting or need a refresher check out <a rel="nofollow" href="https://amzn.to/2PxDb5m" target="_blank"><em>Understanding Statistical Process Control</em></a>, Wheeler et. al. If you want to dive in and start making control charts with R, check out R packages</p>
<ul>
<li><a rel="nofollow" href="http://ggqc.r-bar.net/" target="_blank">ggQC</a>: ggplot based QC charting</li>
<li><a rel="nofollow" href="https://cran.r-project.org/web/packages/qcc/vignettes/qcc_a_quick_tour.html" target="_blank">qcc</a>: base graphics QC charting </li>
</ul>
<p>If your familiar with control charts, you've likely encountered cryptic alpha-numeric constants like d2, A2, E2, d3, D3, and asked, </p>
<blockquote>
<p>“What are they and where do they come from?”</p>
</blockquote>
<p><b>Short (not so satisfying) answer</b>: They are constants you plug into formulas to determine your control limits. Their value depends on how many samplings you do at a time and the type of chart you are making. For example, if you measure the size of <b>5</b> widgets per lot of 50, then your subgroup size, n, is 5 and you should be using a set of control chart constants for n = 5. </p>
<blockquote>
<p>So where do they come from and how are they calculated?</p>
</blockquote>
<p>Read on. </p>
<h3>X-Bar and X-Individuals Constants</h3>
<p>Often, control charts represent variability in terms of the mean range, <span>R</span>, observed over several subgroup rather than the mean standard deviation. The table below should make the idea of subgroup range and mean range more clear. </p>
<p><img src="https://i2.wp.com/cdn.r-bar.net/PID278/rangetable.png?w=456" data-recalc-dims="1"></p>
<p>Why range? My guess is that, historically, employees at all levels would have understood the concept of range. Range requires no special computation, just (max-min). Speculation aside, we begin our quest to understand where control constants come from with the relationship shown in <strong>Eq. 1</strong> that the mean subgroup range is proportional to standard deviation of the individual values. </p>
<p><img src="https://i1.wp.com/cdn.r-bar.net/PID278/EQ1.png?w=456" data-recalc-dims="1"></p>
<p>The proportionality constant between <span>R</span>(X<sub>Sub_Grp_Indv</sub>) and S(X<sub>indv</sub>) is d<sub>2</sub>, the first constant we'll be estimating. The relationship is expressed in <strong>Eq.2</strong></p>
<p><img src="https://i1.wp.com/cdn.r-bar.net/PID278/EQ2.png?w=456" data-recalc-dims="1"></p>
<p>To estimate d<sub>2</sub> for n = 2 (i.e, the subgroup size is 2), we start by drawing two samples from a normal distribution with mean = 0 and sd = 1. Why you ask? Because it makes the math really simple. Consider <strong>Eq. 2</strong>, if S(X<sub>indv</sub>) = 1 then <span>R</span>(X<sub>Sub_Grp_Indv</sub>) = d<sub>2</sub>. </p>
<p>So all we need to do to determine d<sub>2</sub> when n = 2 is: </p>
<ol>
<li> Draw 2 individuals from the normal distribution,
<li> Determine the range of the 2 samples.
<li> Repeat many, many times
<li> Take the average of all the ranges you’ve calculated, <span>R</span>
<li> <span>R</span> = d<sub>2</sub> (when the mean = 0 and sd = 1)
</ol>
<p>The R code for the process is shown below. </p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
</pre></td><td><pre> 
```r
<span>require</span><span>(</span>magrittr<span>)</span>  <span>#Bring in the Pipe Function</span>
reps <span>&lt;-</span> 1E6
<span>set.<span>seed</span></span><span>(</span><span>5555</span><span>)</span>  
<span>replicate</span><span>(</span>reps, <span>rnorm</span><span>(</span>n<span>=</span><span>2</span>, <span>mean</span> <span>=</span> <span>0</span>, <span>sd</span> <span>=</span> <span>1</span><span>)</span> <span>%>%</span> <span>#Draw Two From the Normal Distibution </span>
    <span>range</span><span>(</span><span>)</span> <span>%>%</span> <span>#Determine the Range Vector = (Max, Min)</span>
    <span>diff</span><span>(</span><span>)</span> <span>%>%</span> <span>#Determine the Difference of the Range Vector</span>
    <span>abs</span><span>(</span><span>)</span> <span>#Take the Absolution Value to make sure the Result is positive</span>
    <span>)</span> <span>%>%</span> <span># Replicate the above proceedure 1,000,000 Times </span>
    <span>mean</span><span>(</span><span>)</span> <span>-></span> R_BAR <span>-></span> d2 <span>#Take the mean of the 1,000,000 ranges</span>
d2
```
 
```
<span>## [1] 1.12804</span>
```</pre></td></tr></table></div>
<p>The pipes make the above code easy to read but slow things down quite a bit. The following code does the same thing about 12.5 times faster. </p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
</pre></td><td><pre> 
```r
reps <span>&lt;-</span> 1E6
<span>set.<span>seed</span></span><span>(</span><span>5555</span><span>)</span>
d2 <span>&lt;-</span> R_BAR <span>&lt;-</span> <span>mean</span><span>(</span><span>replicate</span><span>(</span>reps, <span>abs</span><span>(</span><span>diff</span><span>(</span><span>range</span><span>(</span><span>rnorm</span><span>(</span><span>2</span><span>)</span><span>)</span><span>)</span><span>)</span><span>)</span><span>)</span>
d2  
```
 
```
<span>## [1] 1.12804</span>
```</pre></td></tr></table></div>
<p>Once you have d<sub>2</sub>, calculating E2 (3σ for the individuals) and A2 (3σ for the sub-group means) is straight forward as shown in <strong>Eq.3 – Eq.6</strong>. A2 and E3 are the coefficients to the left of <span>R</span>. </p>
<p><img src="https://i2.wp.com/cdn.r-bar.net/PID278/EQ-3-4-5-6.png?w=456" data-recalc-dims="1"></p>
<p>The code below gives the expected results for all the control constants need to construct X-Bar and X-Individual charts. </p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
</pre></td><td><pre> 
```r
<span>c</span><span>(</span>N<span>=</span><span>2</span>, d2 <span>=</span> d2, E2 <span>=</span> <span>3</span><span>/</span>d2, A2 <span>=</span> <span>3</span><span>/</span><span>(</span>d2<span>*</span><span>sqrt</span><span>(</span><span>2</span><span>)</span><span>)</span><span>)</span>
```
 
```
<span>## N d2 E2 A2 </span>
<span>## 2.000000 1.128040 2.659480 1.880536</span>
```</pre></td></tr></table></div>
<h3>R-Bar Constants</h3>
<p>The constants for <span>R</span> charts are d<sub>3</sub> (1σ around <span>R</span>,), D3 (Lower 3σ limit of <span>R</span>) and D4 (Upper 3σ limit of <span>R</span>). To get these constants, we start with the assumption that the standard deviation of <span>R</span> is proportional to the standard deviation of the individual X's. The proportionality constant is d<sub>3</sub> shown in <strong>Eq.7</strong>. Notice that Eq. 7 has the same form as Eq. 2.</p>
<p><img src="https://i0.wp.com/cdn.r-bar.net/PID278/EQ7.png?w=456" data-recalc-dims="1"></p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
</pre></td><td><pre> 
```r
reps <span>&lt;-</span> 1E6
<span>set.<span>seed</span></span><span>(</span>seed<span>)</span>
d3 <span>&lt;-</span> <span>sd</span><span>(</span><span>replicate</span><span>(</span>reps, <span>abs</span><span>(</span><span>diff</span><span>(</span><span>range</span><span>(</span><span>rnorm</span><span>(</span><span>2</span><span>)</span><span>)</span><span>)</span><span>)</span><span>)</span><span>)</span>
d3  
```
 
```
<span>## [1] 0.8529419</span>
```</pre></td></tr></table></div>
<p>Notice in the R code above, the only difference between the calculation of d<sub>3</sub> and d<sub>2</sub> is that we use <strong>standard deviation</strong> rather than the <strong>mean</strong> of the R<sub>Sub_Grp_Indv</sub>. Now we have d<sub>3</sub>, but we need to do a little simple algebra to express the S(R<sub>Sub_Grp_Indv</sub>) in terms of <span>R</span>. Remember, historically the employee doesn't need to worry about standard deviation – just ranges. We can define the above expression in term of <span>R</span> by combining <strong>Eq.2</strong> and <strong>Eq.7</strong>, yielding <strong>Eq.8</strong>.</p>
<p><img src="https://i2.wp.com/cdn.r-bar.net/PID278/EQ8.png?w=456" data-recalc-dims="1"></p>
<p>OK almost to D3 and D4. The lower 3σ limit of <span>R</span> can be expressed as Eq.8:</p>
<p><img src="https://i1.wp.com/cdn.r-bar.net/PID278/EQ9.png?w=456" data-recalc-dims="1"></p>
<p>Factoring out the <span>R</span> terms on the right-hand side of the expression yields </p>
<p><img src="https://i0.wp.com/cdn.r-bar.net/PID278/EQ10.png?w=456" data-recalc-dims="1"></p>
<p>The expression inside the parentheses is D3. D4, the upper limit of <span>R</span> is evaluated analogously. The only difference is a “+” sign in the final expression. Final expressions for D3 and D4 are:</p>
<p><img src="https://i1.wp.com/cdn.r-bar.net/PID278/EQ11-12.png?w=456" data-recalc-dims="1"></p>
<p>All done! Here is the R code summarizing the constants for <span>R</span> using n=2.</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
</pre></td><td><pre> 
```r
<span>c</span><span>(</span>N  <span>=</span> <span>2</span>, 
  d3 <span>=</span> d3,
  D3 <span>=</span> <span>ifelse</span><span>(</span><span>1</span> <span>-</span> <span>3</span><span>*</span>d3<span>/</span>d2 <span>&lt;</span> <span>0</span>, <span>0</span>, <span>1</span> <span>-</span> <span>3</span><span>*</span>d3<span>/</span>d2<span>)</span>, 
  D4 <span>=</span> <span>1</span> <span>+</span> <span>3</span><span>*</span>d3<span>/</span>d2
  <span>)</span>
```
 
```
<span>## N d3 D3 D4 </span>
<span>## 2.0000000 0.8529419 0.0000000 3.2683817</span>
```</pre></td></tr></table></div>
<p>Notice for D3 the value is 0, this is because the value calculated was negative. Such values are rounded to zero per the R code above.</p>
<h3>Summary</h3>
<p>In this post, we used R to estimate the control chart constants needed to produce X-Individuals, X-Bar, and R-Bar charts. All the constants together are shown below. In addition, the constants for n = 7 have also been presented. </p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
</pre></td><td><pre> 
```r
reps <span>&lt;-</span> 1E6
<span>set.<span>seed</span></span><span>(</span><span>5555</span><span>)</span>
FUN_d2 <span>&lt;-</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span><span>mean</span><span>(</span><span>replicate</span><span>(</span>reps, <span>abs</span><span>(</span><span>diff</span><span>(</span><span>range</span><span>(</span><span>rnorm</span><span>(</span>x<span>)</span><span>)</span><span>)</span><span>)</span><span>)</span><span>)</span><span>}</span>
FUN_d3 <span>&lt;-</span> <span>function</span><span>(</span>x<span>)</span> <span>{</span><span>sd</span><span>(</span><span>replicate</span><span>(</span>reps, <span>abs</span><span>(</span><span>diff</span><span>(</span><span>range</span><span>(</span><span>rnorm</span><span>(</span>x<span>)</span><span>)</span><span>)</span><span>)</span><span>)</span><span>)</span><span>}</span>
 
Ns <span>&lt;-</span> <span>c</span><span>(</span><span>2</span>,<span>7</span><span>)</span>
d2 <span>&lt;-</span> <span>sapply</span><span>(</span>Ns, FUN_d2<span>)</span>
d3 <span>&lt;-</span> <span>sapply</span><span>(</span>Ns, FUN_d3<span>)</span>
 
<span>round</span><span>(</span><span>data.<span>frame</span></span><span>(</span>
  N <span>=</span> Ns,
  d2 <span>=</span> d2,  
  E2 <span>=</span> <span>3</span><span>/</span>d2, 
  A2 <span>=</span> <span>3</span><span>/</span><span>(</span>d2<span>*</span><span>sqrt</span><span>(</span>Ns<span>)</span><span>)</span>,
  d3 <span>=</span> d3,
  D3 <span>=</span> <span>ifelse</span><span>(</span><span>1</span> <span>-</span> <span>3</span><span>*</span>d3<span>/</span>d2 <span>&lt;</span> <span>0</span>, <span>0</span>, <span>1</span> <span>-</span> <span>3</span><span>*</span>d3<span>/</span>d2<span>)</span>, 
  D4 <span>=</span> <span>1</span> <span>+</span> <span>3</span><span>*</span>d3<span>/</span>d2  
  <span>)</span>, digits <span>=</span> <span>3</span><span>)</span>
```
 
```
<span>## N d2 E2 A2 d3 D3 D4</span>
<span>## 1 2 1.128 2.659 1.881 0.854 0.000 3.272</span>
<span>## 2 7 2.704 1.109 0.419 0.833 0.076 1.924</span>
```</pre></td></tr></table></div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://r-bar.net/estimating-control-chart-constants-with-r/"> R – R-BAR</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/YtSB_jSmV1o" height="1"></a></div>