Statistical Quality Control with the qcr Package

Abstract The R package qcr for Statistical Quality Control (SQC) is introduced and described. It includes a comprehensive set of univariate and multivariate SQC tools that completes and increases the SQC techniques available in R. Apart from integrating different R packages devoted to SQC (qcc, MSQC), qcr provides nonparametric tools that are highly useful when Gaussian assumption is not met. This package computes standard univariate control charts for individual measurements, x̄, S, R, p, np, c, u, EWMA, and CUSUM. In addition, it includes functions to perform multivariate control charts such as Hotelling T2, MEWMA and MCUSUM. As representative features, multivariate nonparametric alternatives based on data depth are implemented in this package: r, Q and S control charts. The qcr library also estimates the most complete set of capability indices from first to the fourth generation, covering the nonparametric alternatives, and performing the corresponding capability analysis graphical outputs, including the process capability plots. Moreover, Phase I and II control charts for functional data are included.


Introduction
Throughout the last decades, there has been an increasing interest to measure, improve and control the quality of products, services and procedures. This is connected to the strong relationship between quality, productivity, prestige, trust and brand image. In fact, implementing procedures of statistical quality control (SQC) are currently related to increasing companies' competitiveness.
The concept of quality control has been extended from the first definitions based on the idea of adjusting production, to a standard model to satisfy customer requirements and include all participants. Nowadays, SQC is not only applied to manufactured products but to all industrial and service processes. features (CTQs), thus performing process changes when actually required. Furthermore, control charts provide estimates of the natural range of process variation (natural control limits), allowing us to compare this range with those limits specified by standards, company managers or customers (specification limits). Hence, the process monitoring can be carried out by comparing each new observation with these natural limits, preventing defects in the final product. Briefly, a control chart is a two dimensional graph whose axis represents the variable or attribute that is being monitored (CTQ variables). The estimation of natural control limits of the CTQ variables is developed by a process composed of two phases: In Phase I, the natural control limits are estimated using a preliminary sample (calibration sample) where we assume that the causes of variation are only random. In Phase II, each new observation is plotted on the control chart along with the natural limits obtained in the previous step. The set of new observations (twhich are not used to calculate the natural control limits) make up the so-called monitoring sample. Patterns, observations of out of control limits, runs of more than six observations on one side of the central line, among others, are some of the different criteria to identify out of control states in a specific process, providing also valuable information about the detection of any assignable causes of variation in the monitoring.
Anomaly detection-quality control approach to improve and ensure energy efficiency and hygrothermal comfort in buildings The most used control charts are based on the assumptions of normality and independence of the studied CTQ variables. These charts are used to control position and dispersion of CTQ attributes and variables. Figure 2 shows some of the most important types of control charts. These can be classified according to the type of feature that is being controlled (attribute or variable), the variable dimension (univariate or multivariate), and assuming or not a parametric distribution of the variable (parametric or nonparametric). The qcr package provides charts for the mean (x), standard deviation (s), range (R), individual measurements (I), moving ranges (M R), proportion of nonconforming units (p), number of nonconforming units (np), number of defects per unit (c), mean number of defects per control unit (u), exponentially weighted moving average (EWMA), and cumulative sum control chart (CUSUM). The last two techniques are also called memory control charts and they are the specially designed to detect shifts of less than two standard deviations, both when using rational samples or individual measurements. On the other hand, new control charts based on the concept of data depth and developed by ? are implemented in qcr. Those are the r, Q and S control charts, the nonparametric alternatives for individual measurements, mean control chart and CUSUM control chart, respectively. When more than one variable defines the process quality, multivariate control charts are applied. If the Gaussian assumption is met, Hotelling T 2 control chart can be applied. If we want to detect small deviations, multivariate EWMA (MEWMA) and multivariate CUSUM (MCUSUM) can be implemented. When no parametric distribution is assumed, r, Q and S charts can be used.
Another interesting SQC tool, which is very useful in industry, is the Process Capability Analysis (PCA). It estimates how well a process meets the tolerances defined by the company, customers, standards, etc., by comparing the specification tolerances with respect to the natural range of variation of CTQ features. The process capability is measured using capability indicators, Process Capability Ratio (PCR) is a numerical score that helps the manufacturers know whether the output of a process meets the engineering specifications. Large PCR values show that the industrial or service process is capable of meeting the customer requirements. There have been many different PCRs developed in the last four decades that require the Gaussian assumption for the CTQ variable (?). But many processes in industry and real applications do not meet this hypothesis, thus we could innacuratelly estimate the capability using PCR, hence many authors have studied different nonparametric alternatives to traditional PCR (?).
The qcr package has been developed in R (?) under the GNU license. Nowadays there are other R packages that currently provide quality control tools for users. The use of each one is shown in Figure 3.
The qcc package (?) was developed by Professor Luca Scrucca of the Department of Economics, Finance and Statistics at the University of Perugia. It enables us to perform Shewhart quality control charts for variables and attributes, as well as the CUSUM and EWMA charts for detecting small changes in the CTQ variable. Multivariate analysis is performed applying the Hotelling T 2 control chart. Additionally, it has functions implemented to obtain the operating characteristic curves (OC) and to estimate process capability analysis indices. Pareto and Ishikawa diagrams are also implemented. Otherwise, the IQCC package (?) is maintained by Professor Emanuel P. Barbosa of the Institute of Mathematics in the State University of Campinas. It has a smaller number of control charts implemented but it incorporates multivariate graphics. The qualityTools package (?) was developed to aid learning in quality sciences. The Figure 3 shows some of its utilities, e.g. capability analysis (providing a comprehensive set of parametric distributions) and design of experiments. In addition, the SixSigma library (??) provides alternative functions to qualityTools and qcc packages, and the possibility of implementing process maps.  Furthermore, there are other libraries specifically focused on control chart applications. Namely, the spcadjust (?) that allows us to estimate the calibration control limits of Shewhart, CUSUM and EWMA control charts, and the spc (?) which provides tools for the evaluation of EWMA, CUSUM and Shiryaev-Roberts control charts by using Average Run Length and RL quantiles criteria. Moreover, the MSQC package (?) is a set of tools for multivariate process control, mainly control charts. It contains the main alternatives for multivariate control charts such as Hotelling (T 2 ), Chi squared, MEWMA, MCUSUM and Generalized Variance control charts. It also includes some tools to evaluate the multivariate normal assumption. The corresponding multivariate capability analysis can be performed using the MPCI library (?) that provides different multivariate capability indices. It is also interesting to mention the edcc package (?) for its ecomomic design of control charts by minimizing the expected cost per hour of the studied process.
It is important to emphasize that the qcr package also includes new applications such as nonparametric approaches of control charts and capability indices (also covering the capability plots), which are currently unavailable in the R software.

Datasets in the qcr package
The qcr package contains new databases (see table 1) based on study cases tackled by the authors during their professional activity as well as well-known datasets implemented on other packages focused on statistical quality control such as: • archery1: It consists of a stage in which the archer shoots 72 arrows. The information is given in x and y coordinates. It is implemented in the MSQC package (?).
• circuit: Number of nonconformities observed in 26 successive samples of 100 printed circuit boards. It is implemented in the qcc package (?).
• dowel1: Diameter and length of a dowel pin. It is implemented in the MSQC package (?).
• orangejuice: Frozen concentrated orange juice is packed in 6-oz cartons. These cartons are formed on a machine by spinning them from a cardboard stock and attaching a metal bottom panel. A can is then inspected to determine whether, when filled, the liquid could possibly leak either on the side seam or around the bottom joint. If this occurs a can is considered nonconforming. The data were collected as 30 samples of 50 cans each at half-hour intervals over a three-shift period in which the machine was in continuous operation. It is implemented in the qcc package (?).
• pcmanufact: A personal computer manufacturer counts the number of nonconformities per unit on the final assembly line. He collects data on 20 samples of 5 computers each. It is implemented in the qcc package (?).
• pistonrings: Piston rings for an automotive engine are produced by a forging process. The inside diameter of the rings manufactured by the process is measured on 25 samples, each of size 5, drawn from a process being considered 'in control'. It is implemented in the qcc package (?).

Univariate and multivariate parametric control charts in qcr
The construction of a control chart is equivalent to the plotting of the acceptance regions of a sequence of hypothesis tests over time. Namely, thex chart is a control chart used to monitor the process mean µ. It plots the sample means,X's, corresponding to subgroups of the {X 1 , X 2 , ...} observations and is equivalent to test the hypotheses H 0 : µ = µ 0 versus Hα : µ = µ 0 (for some target value µ 0 ) conducted over time, usingx as the test statistic. Here we assume that {X 1 , X 2 , ...} are the sample measurements of a particular CTQ feature that follows the F distribution with mean µ and standard deviation σ. When there is insufficient evidence to reject H 0 , we can state that the process is under control; otherwise, the process is out of control. In other words, processes are under control when their sources of variation are only the sources common to the process (?). The decision to reject or not H 0 is based on the value of the sample meanx observed at each time interval (?). The control charts are easy to construct, visualize, and interpret, and most important, have proven their effectiveness in practice since the 1920's.
Control charts are defined, on the one hand, by a center line that represents the average value of the the CTQ feature corresponding to the in-control state and, on the other hand, two horizontal lines, called the upper control limit (UCL) and the lower control limit (LCL). The region between the control limits corresponds to the region where H 0 is not rejected (defined in the previous section). As a consequence, the process will be out of control when an observed rational sample or an individual measurement falls outside the limits. Let w be a sample statistic that measures a quality characteristic of interest, and suppose that the mean of w is µw and the standard deviation of w is σw. Then the center line, the upper control limit, and the lower control limit become: Where L is the "distance" of the control limits from the center line, expressed in standard deviation units.
When several random variables characterize the quality of a process/service, applying statistical multivariate quality control technique becomes necessary. In fact, if we analyze each variable separately, the probability that an observation of a variable will fall within the calculated limits when it is known that the process is actually under control, will no longer be 0.9973 for 6σ amplitude.

Name Description counters
A water supply company wants to control the performance of the water counters installed throughout a city. For this purpose, 60 rational samples have been taken, each one composed by 3 measurements, from the same age (10 years) and caliber water counters corresponding to two different brands, and during a period of 5 years. This dataset is based on a study case of A Coruña's water supply company, Empresa Municipal de Aguas de La Coruña (Emalcsa).

employment
A Spaniard-Argentinian hotel company wants to control the level of occupancy (measured in %) in their establishments through the application of a continuous control. For this purpose, 48 subsamples have been taken from six hotels corresponding to two different countries.

oxidation
This database contains information about the resistance against the oxidation of olive oil of the Picual variety. Five measurements of the Onset Oxidation Temperature (OOT, index that measures the resistance against the oxidation) are obtained from 50 batches of Picual olive oil produced in chronological order. It is important to note that OOT decreases as the oil is progressively mixed with other olive oil varieties defined by a lower OOT.

plates
A chemical company is developing a patent for a new variant of artificial stone mostly made of quartz (93wt% and polyester resin). This company is launching a pilot plant where it begins to produce plates of this material to industry scale. The CTQ variable of this product is the Vickers hardness. In order to measure the hardness level and hardness homogeneity of the product, 50 plates have been measured 5 times in different sections. The characteristic learning curves, through gradual level change, can be observed.

presion
A shipyard of recreational boats is intended to optimize and control the mechanical properties of the yacht hulls made of a composite based on epoxy resin. At this regard, the modulus of elasticity due to tensile efforts is measured after applying two different curing pressures: 0.1 and 10 MPa. Overall, 60 subsamples, composed of three measurements, obtained from 60 vessels, have been taken. Table 1: Some of the specific datasets included in the qcr package Assuming independence, it will be 0.9973 p , where p is the number of CTQ features, while the probability of type I will actually lead to α = 1 − (1 − α) p . Therefore, the control limits are different from those drawn assuming the control of each CTQ variable independently from the others. Moreover, if the variables are dependent, the calculation of α becomes more complex. This subject is particularly important today, as automatic inspection procedures make it customary to measure many parameters of each product over time. The more common multivariate parametric control charts are the Hotelling T 2 (to identify big shifts), and the multivariate CUSUM (MCUSUM) and EWMA (MEWMA) for identifying small shifts.
The functions that compute the quality control statistics for the different univariate control charts (involving continuous, attribute or count data) are shown in Table 2. For the sake of simplicity, and taking into account that these types of control charts are implemented in other packages, the use of these functions is not shown in this work. More details are given in the help of qcr package (?).

Nonparametric control charts based on data depth
The control charts presented in this section were proposed by ? as an alternative to those described in previous section. The main idea of its control graphs is to reduce each multivariate measure to a univariate index, that is, its relative center-exterior classification induced by a depth of data. This approach is completely nonparametric, therefore, these control charts are not defined by any parametric assumption regarding the process model. Thus, they are applicable in a wider number of case studies than those counterparts such as T 2 , MSUSUM and MEWMA control charts. In addition, these graphs allow the simultaneous detection of the change of location (shift of the mean)  Liu proposed and justified three types of control charts, the r, Q, and S charts which can be considered as data-depth-based multivariate generalizations of the univariate X,x, and CUSUM charts, respectively.

Data depth
In multivariate analysis, the term depth refers to the degree of centrality of a point regarding to a data cloud or a probability distribution. Therefore, it is possible to define a rank in the multidimensional Euclidean space through the calculation of observation depth. According to ? and ?, the depth function can be defined as a bounded function Dp : R d −→ R, with P the distribution set in R d , that assigns at each point of R d its degree of centrality with respect to P . Depth functions with which control charts can be performed are the • Simplicial depth (?), • Mahalanobis depth (?), • Halfspace or Tukey depth (?), • Likelihood depth (?), and • Random projection depth (?).

Statistics derived from data depth
Let G a k-dimensional distribution, and let Y 1 , . . . , Ym be m random observations from G. The sample Y 1 , . . . , Ym is generally the reference sample of a CTQ variable in the context of quality control, composed of measurements from products obtained by an under control process. If X 1 , X 2 , . . . are the new observations from the manufacturing process, assuming that the different X i values follow an F distribution if the quality of the studied product has been deteriorated or, in other words, if the process is out of control. Otherwise they follow a G distribution. Let D G (·) denote a notion of depth, and assume that G and F are two continuous distributions. Thus, if all the D G (Y i ) values are sorted in increasing order, and Y [j] denotes the sample value associated with the jth smallest depth value, then being the most central point. Therefore, the smaller the order (or the rank) of a point, the farther that point will be from the underlying distribution G(·).
? defines the rank statistic as When G is unknown, the empirical distribution Gm of the sample {Y 1 , . . . , Ym} can be used instead, and the statistic is defined by In the same way that r G and r Gm , the Q statistics can be also defined as follows whereby Fn(·) denotes the empirical distribution of the sample {X 1 , . . . , Xn}. The control charts corresponding to these statistics can be developed as described in the following sections.

The r chart
. . , Ym are available. As a result, the r chart consists of plotting the rank statistic in regard to time. The control chart central line is CL = 0.5, whereas the lower limit is LCL =α, with α accounting for the false alarm rate. The process will be out of control if r G (·) falls under LCL. A small value of the rank statistic r Gm (X) means that only a very small proportion of Y i values are more outlying than X. Therefore, assuming that X ∼ F , then a small value of r Gm (X) suggests a possible deviation from G to F . This may be due to a shifting in the location and/or an increase in the scale of the studied CTQ variable. Taking into account that the UCL is not defined for the r chart, the CL line serves as a reference to identify emerging patterns, runs or trends. If r Gm (X) is greater than 0.5, there is an evidence of scale decreasing, and also could take place a negligible location shift. This case should be tackled as an improvement in quality given a gain in accuracy, and thus the process should not considered as out of control.

The Q chart
The idea behind the Q chart is similar to the one behind thex chart. If X 1 , X 2 , . . . are univariate and G is a normal distribution, thex chart plots the averages of consecutive subsets of the de different X i . A goal of this type of chart is that it can prevent the identification of false alarm when the process is actually in control (even when some individual sample points fall out of control limits due to random fluctuations). The Q chart is the nonparametric alternative tox chart. It is performed by plotting the averages of consecutive subsets of size n corresponding to the rank statistic (r G (X i ) or r Gm (X i )), given by Q G, F j n or Q Gm, F j n , whereas F j n is the empirical distribution of the X i 's in the jth subset, j = 1, 2, . . . . Accordingly, if only {Y 1 , Y 2 , . . . , Ym} are available, the Q chart plots the sequence

The S control chart
The S control chart is based on the CUSUM univariate control chart, which is basically the plot of n i=1 (X − µ) which reflects the pattern of the total deviation from the expected value. As mentioned above, it is more effective than the X chart or thex chart in detecting small process changes. The nonparametric CUSUM chart based on data depth suggests plotting Sn(G) and Sn(Gm), defined by with control limits CL = 0 and LCL = −Zα n 12 1 2 , and And, if only Y 1 , . . . Ym are available, the control limits are CL = 0 and LCL = −Zα n 2 ( 1 m + 1 n )

12
. The LCL control limits in both cases constitute a curve instead of a straight line; if n is large, the control chart S should be standardized as follows Therefore, this S * chart is defined by CL = 0 and LCL = −Zα.

Examples of r, Q and S control charts applied using synthetic data
A bivariate data set is used to illustrate how the previously discussed control charts arise. In fact, a synthetic dataset composed of 540 observations of a bidimensional standard Gaussian variable has been simulated, in addition to 40 individuals corresponding to another bidimensional Gaussian variable with mean and standard deviation equal to 2.

r chart
The r control chart can be obtained by applying the npqcs.r() function to the npqcd object and plotting the result R> res.npqcs <-npqcs.r(data.npqcd, method = "Tukey", alpha = 0.025) R> plot(res.npqcs, title = " r Control Chart") The resulting chart is shown in Figure 4, where it can be observed that the process is out of control from the 42nd observation, as expected, taking into account that most of the r Gm (X i ) values are falling below the LCL.

Q chart
In this case, the dataset is assumed to be composed of rational samples of size 4, thus, the Q nonparametric alternative ofx chart is proposed and applied to control the bidimensional process R> n <-4 # samples R> m <-20 # measurements R> k <-2 # number of variables R> x.a <-array( , dim = c(n, k, m)) R> for (i in 1:m) { + x.a[, , i] <-x[(1 + (i -1) * n):(i * n), ] + } R> data.npqcd <-npqcd(x.a, G) R> res.npqcs <-npqcs.Q(data.npqcd, method = "Tukey", alpha = 0.025) R> plot(res.npqcs, title = "Q Control Chart") Figure 5 clearly shows that the process is out of control in the second half, from the 20th rational sample. We can also see that the high random fluctuations of the r chart are attenuated in the Q chart due to the averaging effect.

S chart
Finally, the nonparametric counterpart of CUSUM control chart is performed from the multivariate individual observations R> data.npqcd <-npqcd(x, G) R> res.npqcs <-npqcs.S(data.npqcd, method = "Tukey", alpha = 0.05) R> plot(res.npqcs, title = "S Control Chart") Figure 6 shows that the process is out of control from the 48th observation. Note that the S graph performs better in identifying small changes in a process. In this case, he performance of Q chart is better than the corresponding to the S chart.

Control charts for functional data based on data depth
In the paradigm of Industry 4.0, processes and services are many times described by continuously monitored data of hourly, daily, monthly curves or smooth functions. When the processes are defined by functional data, the authors encourage to apply control charts based on Functional Data Analysis (FDA) in order to implement control and improvement tasks. In this section, the use of the control charts presented in Flores et al. (?). Summarizing, this methodology consist on the proposal of new Phase I and Phase II control charts to be applied in those case study in which the datum unit is a curve. Phase I control chart is based on the computation of functional data depth (specifically Fraiman and Muniz (?), Mode (?), and random projections (?) data depth) from which a data depth control chart is developed. Once the in control calibration sample is obtained, the Phase II control chart based on functional data depth and rank nonparametric control chart can be applied. In addition to the Phase I functional data depth and Phase II rank control charts, plots of functional envelopes from the original curves are provided in order to help to identify the possible assignable causes of out of control states.

Estimating a Phase I control chart for functional data (calibration)
A dataset is simulated in order to illustrate the use of FDA control charts for Phase I and II. A functional mean, mu 0 , and a functional standard deviation, sigma, are defined as shown in (?). A n 0 = 100 hundred curves composed of m = 30 points are simulated. They accounts for the calibration or retrospective sample.
R> fdchart <-fdqcd(mdata) R> plot(fdchart,type="l",col="gray",main="Functional data") The following step is to identify those curves that account for the in control process. This task is done by the application of a Phase I control chart for functional data. This method is implemented in the qcr package by the fdqcs.depth function. Specifically, the arguments and default values for this functions are R> fdqcs.depth.default <-function(x, data.name=NULL,func.depth = depth.mode,nb=200, + type = c("trim","pond"),ns = 0.01, + plot = TRUE, trim = 0.025, smo =0.05, + draw.control = NULL,...) where func.depth is the type of depth measure, by default depth.mode, nb the number of bootstrap resamples, type accounts for the method used to trim the data, trim or pond (?), ns is the quantile to determine the cutoff from the bootstrap procedure (?), plot a logical value indicating that it should be plotted, trim the percentage of the trimming, smo the smoothing parameter for the bootstrap resampling (?), whereas draw.control specifies the col, lty and lwd for the fdataobj, statistic, IN and OUT objects. When the fdqcs.depth function is applied to the curves of calibration sample, the fddep object of fdqcs.depth class is obtained. It is composed by the original data, the depth corresponding to each curve, the lower control limit of the depth chart, the index of those curves out of control, the curves that accounts for the limits of the envelope composed by the deepest cures, and the deepest curve or functional median.  Figure 8 shows the control chart for the depth of the curves (right panel). The LCL is estimated by smoothed bootstrap procedure (?). In order to provide a tool to identify the assignable cause of each out of control curve, the original curves with the envelope with the 99% of the deepest curves are also shown (left panel). The analysis of the shape and magnitude of the curves in and out of bounds can help to associate each curve out of control to an assignable cause, allowing for processes control, maintenance and improvement. The Phase I ends when a calibration sample without curves out of control is obtained. The iterative procedure to obtain an in control calibration sample in shown in the following lines.

Estimating a Phase II control chart for functional data (monitoring)
The next step is to perform the Phase II of process control. The monitoring phase is performed by the application of Phase II control charts for functional data based on multivariate nonparametric control charts. Firstly, a monitoring sample composed by 50 curves is simulated by the following code.
R> fdchart_a <-fdqcd(mdata_a,"Monitoring curves") R> phase2.chart <-fdqcs.rank(fdchart,fdchart_a) R> plot(phase2.chart) R> summary(phase2.chart) Figure 10 accounts for the FDA chart with the calibration sample and its envelope composed by the deepest curves. Moreover, the monitoring sample is also included and compared with the calibration sample by using the FDA chart. In addition, the Phase II rank control chart for functional data is shown including both calibration and monitoring samples or only the ranks corresponding to the monitoring sample. The second population that corresponds with the monitoring sample is identified by the control chart from the first monitored curve (panels below in Figure 10).

Process capability analysis
The analysis of the capability of a process in the case of statistical quality control is done through the calculation of the so called capability. These indices measure whether a process is capable or not of meeting the corresponding technical specifications set by the customer, or the manufacturer, by comparing those with the natural variability of the CTQ variable that characterizes the process. The interpretation of these indices is associated with the result of this relation. Capability indices are generally calculated as the ratio between the length of the specification interval and the natural variability of the process in terms of σ. Large values of these indices mean that the corresponding process is capable of producing articles that meet the requirements of the client, and manufacturers. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 50 100 150 0.0 0.2 0.4 0.6 0.8 1.0

Phase II: Rank Chart
Index q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q In other words, the larger the value of the capability index, the smaller the number of products outside the specification limits.
In this section, we describe the capability indices for processes whose underlying distribution is normal and not normal (exponential, Weibull, etc.). However, it is important to note that the development of programming tools for nonparametric capability analysis is one of the main goals and contributions of the qcr package. In addition to the estimation of capability indices, a graphical output is provided. Based on the proposal of qualityTools package, the qcr graphical output for capability analysis includes a normality test for the CTQ variable, a Q-Q plot, a histogram with the theoretical Gaussian distribution density, parametric and nonparametric estimates of capability indices and a contour capability control chart. In the following lines, parametric and nonparametric capability analysis utilities are described using different examples of application.

Assuming a normal distribution
The most widely used capability indexes in the industry analyze the process capability under the assumptions of stabilized process (in control) and a Gaussian distributed CTQ variable. Table 3 shows the main parametric (assuming Gaussian distribution) indices, namely Cp, C pk , Cpm, and C pmk .
? proposed a general formulation of these indices by an expression that depends on the nonnegative parameters u and v: Whereby d = (U SL − LSL)/2, m = (LSL + U SL)/2, U SL is the upper specification limit, the LSL is the lower specification limit, σ is the theoretical standard deviation, µ accounts for the theoretical mean of the CTQ variable, and T is the specification target (by default the mean between the LSL and U SL). The indices shown in Table 3 are obtained from this expression just considering values of 0 or 1 for u and v: Cp (0, 0) = Cp, Cp (1, 0) = C pk , Cp (0, 1) = Cpm, Cp (1, 1) = C pmk .
The piston rings data set is used to illustrate the calculation of the capability indices using the qcs.cp() function based on the expressions previously described in the Table 3. From the statistics obtained from thex control chart of pistonrings dataset, the γ and β values are estimated and the corresponding capability index is computed R> data("pistonrings") R> xbar <-qcs.xbar (pistonrings[1:125, ], plot = FALSE) R> limits <-c(lsl = 73.99, usl = 74.01) Actual capability with respect to the specification limitsĈ Shifting of the mean with respect to the targetĈ pm =Ĉ p 1+ μ−T σ 2 C pk correction for detecting deviations with respect to the targetĈ pkm =Ĉ pk 1+ μ−T σ 2 Table 3: PCR from first to fourth generation, U SL is the upper specification limit, LSL is the lower specification limit, µ is the real mean,μ is the estimated mean andσ is the estimated standard deviation. Consequently, the obtained results are Cp = 0.3407, C pk = 0.3006, Cpm = 0.3382, and C pmk = 0.2984, respectively. The argument parameters accounts for the u and v values, while object is the type of control chart from which the σ is estimated, limits are the specification control limits, and contour the parameter that indicates when the process capability contour chart is plotted.

Process capability plot
In ? and ?, a graphical method (based on common capbility indices) to analyze the capability of a process is proposed. The goal of using this type of plot (if compared with respect to only capability indices calculation) is to provide immediate information of the location and spread of the CTQ feature, and about the capability to meet the specifications of the corresponding process. When using this chart, a process will be capable if the process capability index is higher than a certain value k, with k > 1. The most used values for k are k = 1, k = 4/3, or k = 5/3, even 2 at a Six Sigma level, taking into account the usual index limits for which a process could be assumed capable. It will also be assumed that the target value matches the center of the specification interval, that is, T = (U SL+LSL) 2 = m. Then, one of the indices defined by the Cp (u, v) family is used, e.g., C pk or Cpm, and the process will be defined as capable if Cp (u, v) > k, given the values of u, v, and k. Also note that if µ = T , all the Cp (u, v) indices are defined by the same expression as the Cp. Moreover, different setting for u, v, and k impose different constraints on process parameters (µ, σ). This can be easily seen through a process capability plot. This graph is a contour plot of Cp (u, v) = k as a function of µ and σ, but it can also be defined as a function of δ and γ, with δ = µ−T d and γ = σ d . The contour line is obtained by rewriting the index Cp (u, v) as a function of δ and γ as follows Cp (u, v) = 1−u|δ| 3 √ γ 2 +v(δ) 2 . Therefore, the Cp (u, v) = k equation is solved, plotting γ depending on the values of δ. The resulting expressions are: When u = v = 0, that is, when we consider the index Cp = k, we have γ = 1 3k , and |δ| ≤ 1. It is important to highlight that the γ axis accounts for the process spread, whereas the δ axis accounts for the process location. The values of the parameters µ and σ which provide values (δ, γ) within the region bounded by the contour line Cp (u, v) = k and the δ axis will provide a larger Cp (u, v) value than k, leading a capable process. Furthermore, values of µ and σ which provide values (δ, γ) outside this region will provide a value Cp (u, v) smaller than k, i.e., a non-capable process. In the case of the process not being capable, this type of plot is useful to understand if the corrective actions have to be performed to decrease the process spread, or the process location (deviation with respect to target) or even when both changes are needed to improve the process capability. This can be observed by observing the distance with respect to the x and y axis. Below are some examples of capability plot application which can be generated through the application of the qcs.cp function with contour=TRUE and k=1 (default values) R> oldpar <-par(mfrow = c(2, 2)) R> qcs.cp(object = xbar, parameters = c(0, 0), limits = limits, + ylim = c(0, 1)) R> qcs.cp(object = xbar, parameters = c(1, 0), limits = limits, + ylim = c(0, 1)) R> qcs.cp(object = xbar, parameters = c(0, 1), limits = limits, + ylim = c(0, 1)) R> qcs.cp(object = xbar, parameters = c(1, 1), limits = limits, + ylim = c(0, 1)) R> par(oldpar) The result is shown in Figure 11. In all the cases the points in red are out of the area defined by the line in blue and the δ axis. Thus, the corresponding process is not capable no matter the capability index that is used. In any case, note that the Cp index is useless in identifying non capable processes due to location shifts with respect to the target. In the same way, the C pk index assumes as capable processes that are far from the target as long as it were close to the specification limits (as shown in Figure 11). Thus, the use of the Cpm and C pmk are recommended due to they take into account both shifts from the target and the spread. In the present case the process is not capable due to the spread rather than the target shift. Therefore, the process changes could be due to decreases in the variability process.

Estimated process capability plot
In practice, the process parameters are unknown and we need to estimate them. We can perform a decision rule based on the sample statistics that provide a sample estimate of the capability index and, finally the so called estimated process capability plot, also called γ * − δ * plot (?). It allows us to decide whether a process is capable or not assuming that µ and σ parameters are unknown and estimated byμ =X andσ 2 = 1 n n i=1 X 2 i −X 2 . They are the maximum likelihood estimators when the CTQ variable of the process is normally distributed, and X 1 , X 2 , . . . , Xn a random sample of a normal distribution with µ mean and σ 2 variance. The qcr package only provides the γ * − δ * plot corresponding to the Cpm index taking into account that the other capability indices do not consider shifts from the target value in their calculations. For the general case, see the work of ?. In order to obtain an appropriate decision rule for the case of Cpm index, we test the hypotheses H 0 : Cpm ≤ k 0 versus H 1 : Cpm > k 0 , usinĝ as test statistic. The null hypothesis will be rejected ifĈpm > cα, where the constant cα is determined by previously defining a signifition test level α . ? showed that the null hypothesis H 0 : Cpm ≤ k 0 can be reduced to H 0 : Cpm = k 0 . Thus, for given values of α and n, the process will be considered capable ifĈpm > cα, with cα > k 0 . ? proved that, when the Cpm index is used, the critical value for a given α is obtained as where χ 2 α,n is the quantile α of a χ 2 distribution with n degrees of freedom. The qcr package includes the qcs.hat.cpm() function to obtain both the theoretical capability plot and the estimated capability plot from sample statistics. Among other options, the user can indicate the control chart from which the estimatesμ andσ are obtained (alternatively,μ andσ can be introduced through mu and std.dev), and the specification limits using limits. Furthermore, the signification level and the capability limit can be modified, as they are set to α = 0.05 and k 0 = 1 by default. The following code illustrates its application to pistonrings data : Process capability plots using the Cp, C pk , Cpm, and C pmk indexes.

Nonparametric capability analysis
Traditional assumptions about data such as normality or independence are frequently violated in many real situations. Thus, in scenarios in which assumptions of normality are not verified, the indices defined in the previous sections are not valid. ? and ? proposed generalizations of Cp (u, v) for the case of arbitrary distributions of data where Fα is the percentile α% of the corresponding distribution and M the median of the process. However, the distribution of the underlying process is always unknown. ? calculated estimates for F 99.865 , F 0.135 and M based on the sample percentiles.
? proposed the following estimator where Up is an estimator for F 99.865 , Lp is an estimator for F 0.135 andM is an estimator for M , obtained from the tables developed by ?.
The qcs.cpn() function of qcr calculates C N p , C N pk , C N pm and C N pmk using the formulation described by ?. The code that illustrates its use is shown below. To obtain the nonparametric capability indices it is necessary to indicate the u and v parameters.

Tools for a comprehensive processs capability analysis
Function qcs.ca() provides a comprehensive information of the capability of a process, summarized trough a graphica output. This function calculates the process capability indices Cp, C pk , C pL , C pU , Cpm, C pmk from a qcs object, assuming a Gaussian distribution. Moereover, it computes confidence limits for Cp using the method described by ?. Approximate confidence limits for C pl , Cpu and C pk are also estimated using the method described in ?, while the confidence limits for Cpm are based on the aproximated method of ?, that assumes the target is the mean of the specification limits. Moreover, the C N p , C N pk , C N pm , and C N pmk nonparametric capability indices are also obtained. There is also an specific box within the summary plot that shows the proportion of observations, and expected observations under the Gaussian assumption, out of the specification limits (nonconforming observations). Further, a histogram of data sample is provided, in addition to the corresponding Gaussian density curves obtained from the sample estimates (one per standard deviation estimate procedure). They are displayed along with the specification limits, a quantile-quantile plot for the specified distribution and a process capability plot obtained from the Cpm index (both using theoretical and sample alternatives). In order to describe the qcs.ca() performance, the following code corresponds to the analysis of the first 125 observations of the pistonrings dataset (the corresponding output is shown in Figure 13).

Conclusions
The qcr package has been developed to provide users with a comprehensive set of functions that manage statistical process control, ranging from univariate parametric analysis to multivariate and FDA nonparametric statistics. This package includes the main types of control charts and capability indices. It combines the main features of reputed SQC packages in R such as qcc and qualityTools with the proposal of a new graphical appearance and the implementation of new SQC tools with increasing importance in Industry 4.0 such as multivariate and nonparametric analysis. In addition to some utilities provided by reference R packages such as qcc, SixSigma and qualityTools, qcr implements very important statistical techniques of Control and Analysis tasks of the Six Sigma procedure that are not included in other libraries. In the case of multivariate control charts, these tools are the MEWMA and MCUSUM multivariate control charts, on the one hand, and the r, Q and S nonparametric control charts based on data depth, on the other hand. In addition, Phase I and Phase II control charts for functional data (monthly, dialy, hourly curves) based on functional data depth, bootstrap procedures and nonparametric rank charts have been also implemented in the qcr package. These control charts for functional data provide tools to control and improve processes when their CTQ variables are obtained as hourly, monthly, daily, yearly smooth curves. It is also very important to note that qcr provides functions to perform nonparametric capability analysis. In addition, the new implementation of the process capability plots for the main parametric capability indices allows us to analyze if improvements in process spread or/and process location are needed to obtain a capable process. The comparison between suppliers, machines, etc., is enabled through capability plots. All these utilities intend to make qcr a useful tool for users of a wide variety of industries, providing a competitive alternative to commercial software.