Efficient power computation for r out of m runs rules schemes

In this article we develop a power computation code in the R language which provides an easy to use tool to researchers in designing Shewhart control charts. It enables researchers to use different existing and newly introduced sensitizing rules and runs rules schemes designed for Shewhart-type control charts for location and spread. The code provides researchers to compute the power for different options of r out of m rules/schemes. The code is flexible to apply for any sample size, false alarm rate, type of control limits (one- or two-sided), amount of shift in the process parameters and a variety of popular distributions for commonly used Shewhart-type control charts (i.e.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{{X}} ,R ,S}$$\end{document} and S2 charts). These mentioned benefits of the developed functional code are only partially found in features of the existing software packages and these programs may be enhanced by adding the features of the developed code as a function in their libraries dealing with quality control charting.


Introduction
Processes show variations in its output.These variations are mainly of two types, namely natural and unnatural.The natural variations in a process (in-control process) are inherent.For the unnatural variations (out-of-control process) there exist special reasons which need to be identified.A timely identification of these special causes of variations helps boosting the quality of a process.In order to address and differentiate these two types of variations we have a very useful method, called statistical process control (SPC), which contains powerful tools like Pareto charts, check sheets, cause and effect diagrams, and control charts.In this SPC tool kit the control chart is the most important and commonly used tool due to its statistical framework.Shewharttype control charts are very popular because of their simplicity.The charts are good in detecting large shifts in the process parameters.We may enhance their detection ability by supporting their design structures with some extra sensitizing rules and runs rules schemes.
Many authors have proposed runs rules schemes for this purpose; see for example Alwan (2000), Klein (2000), Khoo (2004), Koutras et al. (2007), Antzoulakos and Rakitzis (2008) and the references therein.Implementation of such schemes increases the sensitivity of the control chart for smaller shifts, but at the same time complicates their design structure.As a result performance evaluations of these types of charts become a difficult task.A popular performance evaluation criterion for comparing control charts is the power, which is the probability of declaring a process as out-ofcontrol when it is actually out-of-control.
The most popular Shewhart-type control charts (i.e.X , R, S and S 2 charts) generally work with only one sensitizing rule: a process is said to be out of control if the charting statistic exceeds one of the control limits (this is usually denoted by 1 out of 1 or (1/1)).This 1/1 is mainly sensitive for large shifts.By applying extra runs rules in the basic design structure of a control chart gives a push to their detection ability for smaller shifts.However, this will result in some complications.We list the following issues in this respect: (i) difficulties in the power computation due to complicated design structures; (ii) inflation in the false alarm rate in case of simultaneous application of more rules; (iii) biasedness in case of separate use of each rule; (iv) limited availability in software packages for power evaluation of different r/m rules; (v) non-flexibility for other choices of false alarm rates.
In order to overcome these issues with the runs rules Riaz et al. (2011) provided an extended and compact form of the runs rules and recommended their application in industrial processes for efficient process monitoring.
There are software packages available which accommodate a few of the aforementioned runs rules or schemes and even these suffer some of the problems highlighted in (i)-(v) above.Therefore, some efficient code is needed which addresses and overcomes these complications and is flexible for a variety of distributions useful for practical purposes.In this study we develop an efficient code for different normal and non-normal distributions in the R language which addresses issues (i) through (v), particularly with reference to Riaz et al. (2011).The choice of the R language is made due to its friendly environment, easy access, interactive software environment for statistical computations and graphics, its de-facto standard behavior and implementation link with the S programming language, and its wide application among statisticians and researchers of other fields.Some more similar motivations may be found in Croux and Rousseeuw (1992), Nadarajah and Kotz (2007), Höhle (2007), Berens (2009), Harner et al. (2009), Hornik et al. (2009), Tierney (2009) and Spiring (2010).
Before discussing the said code in the R language we review first the methodological details of all the rules and procedures of Riaz et al. (2011).
We define a runs rules scheme as: "At least r out of m consecutive (r/m) points fall outside its respective signaling (control) limits h r/m which we may set on one side (lower or upper) or on both the sides of the sampling distribution of a control charting statistic.In this study we consider m = 1, 2, 3, . . ., 9 and for each m we consider the choices r = m − 2, m − 1, m; where 1 ≤ r ≤ m.In this way we have 24 possible rules of decision.In general terminology we refer these rules as r/m rules.In this set of rules the phrase "at least" is in fact the suggested modification for these runs rules schemes which helps to overcome the above mentioned problems (ii), (iii) and (v) and gives an attractive independence to each rule.The signaling limit h r/m is set according to a pre-specified value of false alarm rate for a given value of n and it is obtained using the mathematical expression α = r ) , where α is the pre-specified false alarm rate and p is the probability of a single point falling outside the respective signaling limits depending upon r and m.For more details the interested reader is referred to Riaz et al. (2011).
Different software packages, like MINITAB, SPSS, MATLAB, MAPLE, R, and S, do not have built-in functions for power computations of the different runs rules schemes on different types of control charts under a variety of distributions and also do not allow flexibility for fixing α at the desired level.An R package for statistical quality control namely "qcc" is given by Scrucca (2004) to monitor process characteristics.It is an efficient package for control charting of the charts like X , R, S and S 2 charts, but it works with only one sensitizing rule (i.e.1/1) and misses the application of the other rules which are really needed to address the smaller shifts.The power computations accommodated by Scrucca (2004) (and similarly for many other packages e.g.MINITAB) are limited to very few rules/schemes (like 1/1, 2/3, or 4/5) and with a limited scope (such as: not providing flexibility of picking any choice of 'r ' and 'm'; not accommodative for choosing probability limits at any desired false alarm rate; not allowing to compute the detection power for any amount of shift in the parameter of interest), while the other rules/schemes, which may be more powerful at detecting different shifts, are missing in their environments.Therefore, we need an efficient code for power computations for all the rules in general, which may work for the commonly used distributional environments.Practitioners generally prefer statistical techniques (e.g. a control chart) which have a high power and they use it for their research proposals (cf.Mahoney and Magel 1996).So the power computation code of this study would be of great value for their future studies.Literature supporting the power evaluation criterion may be also found in Albers and Kallenberg (2006), Riaz (2008) and Montgomery (2009).
The functional code is accommodative for any sample size (n), false alarm rate (α), amount of shift in the parameters of interest, any type of control limits (one-or two-sided) for the design structures of X , R, S and S 2 charts under a variety of distributional setups.The R code is written by using some built-in and some user-defined functions to program the aforesaid r/m runs rules schemes for all the mentioned choices of r and m to achieve the desired objectives.The code is flexible to cover different popular probability distributions of practical importance including Normal, t, Gamma, Chisquared, Weibull, Logistic, and Lognormal (cf.Schoonhoven and Does 2010; Abbasi and Miller 2011 and the references therein).It may be easily modified for other distributions as well.The code is available on the website: http://www.ibisuva.nl/en/Research/computer-codes.html

Description of the code
In this section we provide a detailed description of the code in terms of its development and functionality for the desired purposes.This description will help researchers to implement it easily for their desired use.The code mainly consists of three userdefined functions, namely "onepoint", "controllimit" and "power" to meet the desired objectives.
The function "onepoint" (onepoint = function(alpha, r, m)) provides the probability of a single point falling outside the respective signaling limits.It enables the researcher to fix α at the desired level by applying the index search approach using the arguments "alpha", "r "and "m".The argument "alpha" is the pre-specified false alarm rate and it may take any value between 0 and 1, the argument "r " refers to the favorable points for an out-of-control signal and the argument "m" refers to the total consecutive points considered in a given rule.
The function "controllimit"(controllimit = function(chart,r,m,side,n,simu, mean,sd,alpha, dist)) results in the signaling (control) limit(s) for a given r/m rule.It depends upon the arguments "chart", "r ", "m" "side", "n", "simu", "mean", "sd", "alpha" and "dist".The argument "chart" can take one of the four possible choices, namely "xbar", "S", "R", "S2".The argument "side" relates to the nature of control limits and it may take one of the three possible options namely "L" for lower sided, "U" for upper sided and "T" for two sided limits.The argument "n" is the sample size and it may take any value.The argument "simu" shows the number of simulations to be used to compute the desired power of a control chart.The argument "dist" stands for the particular probability distribution under investigation.The arguments "mean" and "sd" represent the specified values of standard deviation and mean respectively of the normal distribution (which is the default distribution in our code and the default values for "mean" and "sd" are taken as 0 and 1 respectively).The other arguments are as defined earlier.
The function "power" power = function(chart, r, m, side, n, simu, mean, sd, alpha, delta, dist, df, shape, scale, meanlog, sdlog, location) is the front end function for the researcher.The arguments used in this function are "chart", "r", "m", "side", "n", "simu", "mean", "sd", "alpha", "delta", "dist", "df", "shape", "scale", "meanlog", "sdlog" and "location".The argument "delta" represents the amount of shifts in terms of sigma units (i.e.delta*sd) which we want to detect for a given chart.In the case of the X chart, delta = 0 refers to an in-control situation and otherwise out-of-control.For the other three charts delta = 1 refers to an in-control situation and otherwise out-of-control.The argument "dist" represents the particular probability distribution we are concerned and once it is entered as an argument then its corresponding parameter values need to be specified as follow up arguments.The following table may be used to enter the arguments about distributions and their corresponding parameters.By entering the corresponding parameters of a specific distribution the other parameter related arguments automatically become null.The arguments "mean", "sd" and "dist" by default are set as zero, one and normal respectively in the power function because of its popularity in literature.All the other arguments used in "power" function are as defined earlier and it calls the function "controllimit" to use the control limits.The "power" function obtains the corresponding statistic values according to the choice of "chart" and then these statistics values are consecutively compared with the values of the control limit obtained from the function "controllimit".The number of times a point goes out of control is counted.The function "power" ends with an output command (statement) in its coded version which provides the summary of the final output of the desired control chart in the form of probability of detecting out-of-control signals (i.e.power).A summary of the output is shown in the R editor or console in an attractive format for the user.The order in which the code should be carried out is "onepoint" -> "controllimit" -> "power".

Execution procedure of the code
In this section we list the necessary steps which help the users and researchers to adopt our developed code for computing the control limits of a control chart for a given sample size and finally evaluating the detection (power) ability for different probability distributions.
(i) Run the code (as provided on the web link http://www.ibisuva.nl/en/Research/computer-codes.html) in the R editor or console and save the workspace so that it may be loaded before using the "power" function next time (or alternatively "power" may be added in the R library as a function of quality control package).(ii) To get the final output in the form of probability distribution, control chart, control limit, rule type, false alarm rate, delta and power, execute the following function "power" after giving the desired options/values to the arguments used in it e.g.power(chart = "xbar",r = 1,m = 1,side = "T",n = 5,simu = 10000, alpha = 0.0027, delta = 0, dist = "cauchy",location = 10,scale = 10) is for the X chart using the 1/1 rule at α = 0.0027 using two sided control limits for sample size n = 5 and an in-control situation (delta = 0).The other options for these arguments may easily be entered following the description of Sect. 2 as per requirements.(iii) After the execution of the "power" function we will get the final output of the data which researchers and users may want in the form of summary display in the R-console providing the information of distribution (normal, t, chisquared, logistic, weibull, lognormal, gamma), control chart type, control limits type (one-or two-sided), runs rules type (r/m), false alarm rate, delta and power.This summary in the display is easy to understand and is very attractive for researchers to quickly have an idea of the chart's ability.

Illustrations and demonstrations of the code's results
To exemplify the application of our developed code for the said purposes we apply our code for a few of the runs rules schemes with different distributions at different false alarm rates using some choices of delta (shifts) for the control charts considered in this study.Also we show what kind of output is obtained by our code.These results may be compared with those cases for which theoretical results are also available for the validation of our power computation code.
For the 1/1 rule on the X chart we upload the code (as provided on the web link http:// www.ibisuva.nl/en/Research/computer-codes.html) in the R console and execute the following statement: power chart = "xbar , r = 1, m = 1, side = "T , n = 5, simu = 10000, alpha = 0.0027, delta = 0) which results in the following output in the R console: Summary for out-of-control probability We can see in the above summary that the code has given the same false alarm rate as desired and the value of mean, sd and dist are not provided in power function as it the default case of our code as mentioned in Sect. 2.  , 2, 3); (iii) in general a higher choice of "m" results into more efficiency and for a specific choice of "m" a superior performance is experienced when "r " is chosen as "m − 2" followed by "m − 1" and finally "m" for X , R, S charts under different distributional setups (cf.Figs. 1, 2, 3).

Summary, conclusions, recommendations and future proposals
It has never been an easy task to fix the false alarm rate at any level with different runs rules schemes to compute the power of a control chart for a specific probability distribution.The existing software packages generally work with a few sensitizing rules/runs rules schemes for very few distributional environments of practical importance and at very limited choices of false alarm rates.We have developed a functional code which can handle this for many rules or schemes and for any choice of the false alarm rate under a variety of commonly used distributions (this we have done for the X , R, S and S 2 charts at the moment and it may be done for other charts as well).
The code provides an attractive summary display of the power along with the control limits of the corresponding control charts at any false alarm rate and for any amount of shift in the R console directly for different distributions including Normal, t, Gamma, Chisquared, Weibull, Logistic, and Lognormal This code may be easily modified for many other distributions other than those considered in this article.Moreover, it may be extended for the EWMA and CUSUM charts.Also the code may be modified to obtain the average run length (ARL) of the control charts.Further, the code may be enhanced to accommodate a list of values for delta simultaneously and then display the power/ARL results in the form of tables or graphical displays, i.e. their corresponding curves.
Keeping in mind the above mentioned features of the proposed functional code we suggest its inclusion in the library of the R language as a quality control tool/package with the name "power" just like the "qcc" package of Scrucca (2004).In general, it may be considered as a contribution for any package in the form of functional code for the said purposes.The application of this code may be particularly of more benefit for the researchers who rely on the power criterion for the selection of a statistical procedure for their research proposals and hence may use this functional code in their future studies (cf.Mahoney and Magel 1996).