Phase I control chart for individual autocorrelated data: application to prescription opioid monitoring

Abstract

Phase I or retrospective process monitoring plays a key part in an overall statistical process monitoring (SPM) regime and is increasingly emphasized in the recent literature. At present, a lot of the data in a variety of settings (public and private sector organizations) are collected individually and sequentially and thus are serially correlated (or autocorrelated). Though a reasonable amount of work is available in the control charting literature for prospective (Phase II) autocorrelated data monitoring, very little work exists for the retrospective phase (Phase I). In this article, we present a Shewhart-type control chart for Phase I monitoring of individual autocorrelated data, assuming normality, with estimated parameters. The methodology, while developed and presented for the first-order autoregressive (AR(1)) model for simplicity, may be adapted to more general time series models. The correct charting constants, adjusted for autocorrelation and parameter estimation, are derived, and tabulated for a nominal in-control (IC) false alarm probability (FAP). Simulation results show that the proposed chart is favorably IC FAP robust and effective for reasonably small sample sizes, moderate autocorrelation, and some model miss-specifications, compared to other approaches. An illustration using some public health data involving prescription fentanyl transactions is provided to show the potential for broader areas of applications of the proposed methodology. Along with a summary and recommendations, some future research areas are indicated. An R package is developed and made available for implementing the proposed methodology on demand.

Keywords:

• model miss-specification
• parameter estimation
• phase I monitoring
• retrospective control charts
• robustness
• serial correlation
• statistical process monitoring

1. Introduction

Statistical process monitoring (SPM) is commonly divided into two phases: Phase I, the retrospective phase, and Phase II, the prospective phase. In Phase I, often, historical (retrospective) data are collected, studied, and analyzed to glean information about the phenomenon of interest. The aim of this phase, sometimes referred to as the exploratory stage, is to better understand the process, establish process stability and detect outliers using a control chart while maintaining an acceptable overall rate of false alarms. Phase I control charts are often applied on a trial and error basis (Montgomery 2005) to flag and remove potential out of control (OOC) points, identify assignable causes, and take any corrective actions. This process continues with control limits being updated as necessary for the remaining data, until all plotted points fall within the control limits and display a random pattern. At this point, the process is declared to be in-control (IC) and the final set of data is called the reference or baseline data, which is then used to fit a parametric distribution and estimate any of its unknown parameters. In case a parametric distribution can’t be assumed, that is in the non-parametric case, the reference data can be used to estimate the underlying unknown distribution. With this initial phase completed, the control chart is deployed and prospective (new, incoming) data are monitored, sequentially, in Phase II. Chakraborti, Human, and Graham (2008), Capizzi and Masarotto (2013), (denoted CM), Jones-Farmer et al. (2014), Chakraborti and Graham (2019b) and others have discussed the importance and utility of a Phase I analysis, where control charts play an important role. In fact, Jones-Farmer et al. (2014) stated that “In some cases, the process knowledge and resulting process improvement from a Phase I analysis is sufficient and there may be no need for on-going Phase II monitoring except to ensure that the process improvement is sustained.”

There has been a large amount of work done on control charts for autocorrelated data over the last few years. We summarize the current state of the art from an extensive literature review, in Table S1, provided in the Supplementary Materials, for the univariate case, citing some of the key papers including both normal distribution and distribution-free methods. Note that the data in Table S1 reveal that only a very few Phase I control charts are currently available and they are mostly for the case of independent and identically distributed (i.i.d.) data. While the i.i.d. case is an important step forward, a significant amount of the data collected for monitoring are autocorrelated in many areas of applications. As far back as 1992, Alwan (1992) stated “it appears that the assumed underlying ideal of statistical control, that is, iid behavior, is rarely attained; substantial autocorrelated behavior is very common.” Almost thirty years later, Perry and Wang (2020) noted that industry processes today “are often computerized, and therefore, data acquisition via sensor technology is quite common. Data captured in modern industry are often collected very rapidly, thus inducing autocorrelation into the data sequence.” The point is that the autocorrelation, when it exists and is not removable, should be properly accounted for while designing the control charts, otherwise, many more false alarms may be observed which may disrupt and erode the efficacy of the overall monitoring process. Findings from Maragah and Woodall (1992), Köksal et al. (2008), Capizzi (2015) and Dasdemir et al. (2016) support the same premise: when the effects of autocorrelation are ignored, the odds of false alarms and degraded detection capability are increased.

While SPM with autocorrelated data has been studied in the literature, our literature review summarized in Table S1 of the Supplementary Materials indicates that the vast majority of these are for subgrouped data in Phase II and with i.i.d. observations. In fact, Jones-Farmer et al. (2014) stated “It is important to note that most Phase I statistical methods can only be applied with subgrouped observations; however, process-quality characteristics are often observed and recorded as individual observations (n = 1).” Even so, there have been only a few papers on monitoring individual autocorrelated data. For example, Dasdemir et al. (2016) and Weiß et al. (2018) considered the effects of autocorrelation, Shewhart-type charts and analyzed the impact of Phase I sample size on Phase II chart performance. However, in the current literature, there is no Phase I control chart focusing on (1) monitoring individual autocorrelated data with unknown parameters, designed using the false alarm probability (FAP) metric, which is defined as the probability of at least one false alarm and is recommended in Phase I (see Chakraborti and Graham 2019b; Jones-Farmer et al. 2014; Yao and Chakraborti 2021b), and (2) the area of public health, where a clear separation of Phases I and II is recommended (Woodall 2006). Motivated by this, we consider a Shewhart-type Phase I control chart for individual data following an autoregressive order 1 (AR(1)) process, which is one of the simpler and popular models for autocorrelated data, with unknown parameters. The charting constants of the proposed chart are determined for a nominal FAP value. We also provide an R package for an easy deployment of this methodology in practice, on demand. More details about our specific contributions are provided later.

At present, there are two main approaches for statistical monitoring of autocorrelated data. The first is to fit a class of time series models (such as ARMA: autoregressive moving average) to the data and monitor the resulting residuals (see, e.g., Alwan and Roberts 1988; Apley and Tsung 2002; Box, Jenkins, and Reinsel 2008; Lee and Apley 2011; Lu and Reynolds 1999; Montgomery and Mastrangelo 1991; Qiu, Li, and Li 2020; Wardell, Moskowitz, and Plante 1992, 1994; Zhang 1997, 1998, 2000) with a control chart. Alwan and Roberts (1988) called this type of charts the special-cause chart (SCC), also known as a residual chart. Although popular in the literature, questions have been raised about the residual chart, for example, about their interpretability, insufficiency in detecting some assignable causes, efficiency, effect of model misspecification, dependence among the residuals and loss of information (see, e.g., Faltin and Woodall 1991; Knoth and Schmid 2004; Macgregor 1991; Ryan 1991; Schmid 1995). Note however that, currently there is no available Phase I chart based on residuals. The second approach is to monitor the observations directly, assuming specific time series models, such as ARMA models, for the data, paying particular attention to the variance-covariance structure while constructing the control limits (see, e.g., Apley and Lee 2003, 2008 (denoted by AL); Alshraideh and Khatatbeh 2014; Garza-Venegas et al. 2018; Lu and Reynolds 1999; Vasilopoulos and Stamboulis 1978; Weiß et al. 2018). This type of chart has been referred to as the common-cause chart (CCC) in Alwan and Roberts (1988). We consider this direct approach and our charts may be viewed as a generalization of the charts in Champ and Jones (2004), Yao, Hilton, and Chakraborti (2017) and Yao and Chakraborti (2021a, 2021b), from the case of i.i.d. subgroup and individual normal data to the case of autocorrelated individual normal data from an AR(1) process with unknown parameters. Further generalizations can be made for observations from other time series processes.

The article is organized as follows. In Section 2, we first review the Phase I Shewhart-type control chart for i.i.d. individual normal data with unknown parameters and then consider the proposed Phase I chart for individual data from an AR(1) process, also with unknown parameters. In Section 3, we develop the methodology and determine the charting constants. In Section 4, we examine the effects of model misspecification on the proposed chart and adjusted versions of some existing charts in the literature. In Section 5, we illustrate the proposed methodology using a real-world dataset. Finally, in Section 6, we conclude with a summary and some recommendations.

2. Retrospective (phase I) monitoring of individual observations (measurements) with unknown parameters

We start with a quick review of the existing charts for i.i.d. individual data and then discuss the rationale of the proposed methodology by examining the implications for autocorrelated individual data from an AR(1) process. To this end, note that once a statistical model is chosen for the process output that is being monitored, one of the important considerations for setting up control charts is to determine whether the underlying parameters are known (specified) or unknown (to be estimated; most common situation). In addition to the choice and the adequacy of the model, it is now well-understood that parameter estimation plays a vital role in the success of the monitoring regime. Estimation generally degrades chart performance (more false alarms, for example) relative to the nominal, sometimes significantly, which can create major issues with chart implementation and interpretation. Two important assumptions are independence and the data distribution, typically assumed to be normal. Both may be of concern for individual data, which arises in many modern monitoring environments. It is well-known that when the data are dependent or not normally distributed, the performance of standard (traditional) control charts may be compromised in terms of the false alarm rate and shift detection ability.

2.1. Basic Shewhart model: independent and identically distributed (i.i.d.) normal observations

Let $Y_1,Y_2,\dots .,Y_m$ denote a sample of $m$ individual observations, where the i-th observation can be written as $$Y_i={\mu }_i+{\epsilon }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m$$ with ${\mu }_i$ being the mean and ${\epsilon }_i$ being the random error. The random error terms are assumed to be i.i.d. normally distributed with mean 0 and variance ${\sigma }^2$ (${\epsilon }_i\sim \mathit{iid} N\left(0,{\sigma }^2\right)$). In the IC case, ${\mu }_1={\mu }_2=\dots = {\mu }_m={\mu }_0$ (say) and the model reduces to

(1) $$Y_i={\mu }_0+{\epsilon }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m.$$(1)

Eq. [1(1) $$Y_i={\mu }_0+{\epsilon }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m.$$(1) ] is referred to as the basic (Shewhart) model in Montgomery (2005, p. 438). For this model, when the values ${\mu }_0$ and $\sigma ={\sigma }_0$ are known (denoted Case K), the standard Shewhart chart for the mean is given by (2) $${\mu }_0\mathrm{ }\mp 3{\sigma }_0$$(2)

When the parameters $\mu$ and $\sigma$ are unknown (denoted Case U), they are estimated from some data. Letting $\overline{Y}$ and $\widehat{\sigma }$ denote the estimates of $\mu$ and $\sigma ,$ respectively, the Shewhart chart for the mean in Case U is given by (3) $$\overline{Y} \mp a\widehat{\sigma }$$(3) where $a>0$ is the charting constant that needs to be determined.

It is now well-known that “3” is not always the correct charting constant, so that the 3-sigma chart is not the correct Shewhart control chart, in both phases I and II, in the sense that the constant 3 generally leads to a higher rate of false alarms (see, e.g., Quesenberry 1993. Jardim, Chakraborti, and Epprecht 2020) than a typical nominal value. Thus, in Case U, the correct charting constants need to be derived, adjusting for the effects of parameter estimation. This derivation depends on the type of data at hand (i.i.d. or not, subgroup or individual), the phase of monitoring, the monitoring framework (e.g., the type of control chart, change-point, or self-starting formulations), the assumed distribution, the estimators used and an appropriate IC performance metric.

2.2. Individual phase I data (observations) from a stationary AR(1) process

When individual observations are collected sequentially, in close physical (spatial) and/or temporal proximity, there is often a natural (serial) dependence (autocorrelation) among the observations. Given this fact, a time series model is often more appropriate for the analysis of such data. Also, the known (Case K) or the unknown (Case U) status of the parameters is an important consideration.

For simplicity, consider data from a stationary AR(1) process which has been a popular choice in the literature. In this case, in contrast to Eq. [1(1) $$Y_i={\mu }_0+{\epsilon }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m.$$(1) ], the $i$-th observation in Case K is given by (4) $$(Y_i-{\mu }_0)={\phi }_0(Y_{i-1}-{\mu }_0)+{\tau }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m,\mathrm{ }\left|{\phi }_0\right|<1\text{ and }{Y}_0={\mu }_0$$(4) where ${\mu }_0$ is the true IC process mean, ${\phi }_0$ is the true IC autoregressive parameter, and ${\tau }_i$ is the random error assumed to follow a normal distribution with mean 0 and variance ${\sigma }_0^2,$ the quantity ${\sigma }_0^2$ is the true IC variance of random errors. In the sequel, it is convenient to work with the standardized observation defined by (5) $$X_i=\frac{Y_i-{\mu }_{\mathrm{o}}}{\sqrt{{\gamma }_{\mathrm{o}}}},\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m,\mathrm{ }$$(5) where ${\gamma }_0=\frac{{\sigma }_0^2}{1-{\phi }_0^2}$ is the true IC process variance of $Y_i.$ The standardized IC AR(1) process in Case K may be rewritten as (6) $$X_i={\phi }_0{X}_{i-1}+{\epsilon }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m\text{ and }{X}_0=0$$(6) where ${\epsilon }_i=\frac{{\tau }_i}{\sqrt{{\gamma }_{\mathrm{o}}}}$ is the normal random error with mean 0 and variance $1-{\phi }_0^2.$ The Shewhart-type control chart is given by plotting $X_i$ with the control limits (7) $$\mathrm{LCL}=-c,\mathrm{ }\mathrm{CL}=0,\text{ and UCL}=c$$(7) where $c>0$ is the charting constant. The $X_i, i=1,2,\dots ,m,$ are referred to as the charting (plotting) statistics. Transforming back, the control limits for the original observations $Y_i,$ in Case K, may be calculated from (8) $$\mathrm{LCL}={\mu }_{\mathrm{o}}-c\sqrt{{\gamma }_{\mathrm{o}}},\mathrm{ }\mathrm{CL}={\mu }_{\mathrm{o}},\text{ and UCL}={\mu }_{\mathrm{o}}+c\sqrt{{\gamma }_{\mathrm{o}}}$$(8) and the control chart may be equivalently deployed by plotting the $Y_i$ against these limits.

In practice, the parameters are often unknown (Case U) and must be estimated to set up the control chart based monitoring regime. The related topic continues to be a very active research area in recent years. Some disadvantages of existing control charts for autocorrelated Phase I data are: (1) The control charts (see, e.g., Schmid 1995; Vasilopoulos and Stamboulis 1978) are mainly designed and adjusted using the IC average run length (ARL) as the performance metric instead of the false alarm probability (FAP), which is well-known to be more appropriate in Phase I (see, e.g., Chakraborti, Human, and Graham 2008; Jones-Farmer et al. 2014) and (2) The effect of the parameter estimation needs to be considered when building the control charts (Dasdemir et al. 2016; Garza-Venegas et al. 2018; Jones 2002). However, there are several choices for the estimators. For example, Schmid (1995) discussed a modified Shewhart chart by estimating ${\mu }_{\mathrm{o}}$ and ${\gamma }_{\mathrm{o}}$ using the sample mean and the sample standard deviation $${\widehat{\mu }}_0^{\mathrm{MOM}}=\overline{Y}=\frac{1}{m}{\sum }_{i=1}^m{Y}_i$$ and (9) $${\widehat{\gamma }}_0^{\mathrm{MOM}}=S=\frac{1}{m}{\sum }_{i=1}^m{(Y_i-\overline{Y})}^2,$$(9) respectively. Garza-Venegas et al. (2018) considered ${\mu }_0$ and ${\sigma }_0^2$ to be known and recommended ${\phi }_0$ to be estimated using the CSS (conditional sum of squares or least-square) estimator. Dasdemir et al. (2016) considered estimating the parameters using the MLE (maximum likelihood estimation) due to the robustness of the estimators.

Following these findings and considering the popularity of different estimators in the literature, we select and use a relatively robust method to estimate the unknown parameters. The method estimates ${\mu }_0$ and ${\gamma }_0$ using the method of moments (MOM) (Schmid 1995) and ${\phi }_0$ using the MLE due to its better robustness (Dasdemir et al. 2016) and being more accurate especially for small samples (Box, Jenkins, and Reinsel 2008, p. 245). Note that we also studied a method where the parameters are all estimated using the MLE and one where the ${\mu }_0$ and ${\gamma }_0$ are estimated using the MOM and ${\phi }_0$ using the CSS. Results for these two methods are not presented since they do not outperform the proposed method.

This article makes three major contributions. First, the Phase I Shewhart control chart for individual data is developed for the AR(1) process, “corrected” for parameter estimation and autocorrelation, under the normality assumption and using the FAP metric. This involves finding the charting constants currently unavailable in the literature. Second, the IC FAP-robustness of the proposed chart is studied with respect to sample size, magnitude of autocorrelation and model miss-specification. The proposed chart is compared with some existing charts originally proposed for Phase II applications, adapted for use in Phase I, using a Bonferroni adjustment for the FAP (see, e.g., Champ and Jones 2004). This includes the traditional I-MR control chart (Montgomery 2005; referred to as the classic chart hereafter), the modified EWMA chart with a smoothing constant $\lambda =1,$ which is a Shewhart-type chart (Apley and Lee 2008, referred to as the AL chart hereafter), and the residual control chart (see, e.g., Alwan and Roberts 1988). Additionally, a distribution-free Phase I chart for i.i.d. data, the recursive segmentation and permutation (RS/P) chart (Capizzi and Masarotto 2013, referred to as the CM chart hereafter), is included in the comparison. Third, the proposed methodology is illustrated using some opioid-related retail transaction data from the US Drug Enforcement Administration’s Automated Reports and Consolidated Ordering System (ARCOS) (Rich et al. 2020) database. This should help broaden the applicability of the proposed Phase I monitoring method in practice. Finally, some conclusions, recommendations and areas of future research are presented.

3. Methodology

A key to studying the properties of a control chart is to consider the so-called non-signaling event (NSE) and its probability. The NSE in Phase I, is the event where the control chart does not signal, that is, when the charting statistics all plot between the control limits. The $\mathrm{FAP}$ is equal to one minus the probability of the NSE.

3.1. Phase I control chart for the individual data from an AR(1) process in case U

In Case U, the three unknown parameters are estimated using the proposed method described above. Thus the NSE, in Case U, for the Phase I chart with the original observations is given by (10) $${\widehat{\mu }}_0^{\mathrm{MOM}}-c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}\le Y_i\le {\widehat{\mu }}_0^{\mathrm{MOM}}+c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}\mathrm{ }i=1,2,\dots ,m$$(10) where $c^{(U)}$ denotes the charting constant in Case U. The probability of the NSE is given by (11) $$P\left(({\widehat{\mu }}_0^{\mathrm{MOM}}-c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\le \underset{¯}{Y}\le ({\widehat{\mu }}_0^{\mathrm{MOM}}+c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\right)=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)=E_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m{\widehat{\phi }}_0^{\mathrm{MLE}}\right)\right)={\int }_{-1}^1\left[\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\dots \underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}{f}_{\mathit{mvn}\left(\underset{¯}{0},\widehat{\underset{¯}{\mathrm{\Sigma }}}\left({\widehat{\phi }}_0^{\mathrm{MLE}}\left(t\right)\right)\right)}\left(x_1,x_2,\dots ,x_m\right)d{x}_{1}d{x}_2\dots d{x}_m\right]f_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(t\right)dt$$(11) where ${\underset{¯}{X}}^{\left(U\right)}={\left[X_i^{\left(U\right)}\right]}^{\prime },\mathrm{ }i=1,2,\dots ,m$ denotes the vector of estimated standardized observations, $X_i^{\left(U\right)}=\frac{Y_i-{\widehat{\mu }}_0^{\mathrm{MOM}}}{\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}}$ is the $i$ th standardized observations. ${\widehat{\mu }}_0^{\mathrm{MOM}}$ and ${\widehat{\gamma }}_0^{\mathrm{MOM}}$ are the MOM estimators as in Eq. [7(7) $$\mathrm{LCL}=-c,\mathrm{ }\mathrm{CL}=0,\text{ and UCL}=c$$(7) ], ${\widehat{\phi }}_0^{\mathrm{MLE}}$ is the ML (maximum likelihood) estimator, $f_{{\widehat{\phi }}_0^{\mathrm{MLE}}}$ denotes the p.d.f. of ${\widehat{\phi }}_0^{\mathrm{MLE}},$ and ${\underset{¯}{1}}_m$ denotes a column vector with $m$ ones. Using the well-known Yule–Walker equation (see, e.g., Box, Jenkins, and Reinsel 2008, p. 60) or AR(1) correlation (see, e.g., Højsgaard, Halekoh, and Yan 2006), the estimated covariance matrix, which is a function of ${\widehat{\phi }}_0^{\mathrm{MLE}},$ is given by (12) $$\widehat{\underset{¯}{\mathrm{\Sigma }}}\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)=\left[\begin{array}{cc}\begin{array}{ccc}1& {\widehat{\phi }}_0^{\mathrm{MLE}}& {\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^2\\ {\widehat{\phi }}_0^{\mathrm{MLE}}& 1& {\widehat{\phi }}_0^{\mathrm{MLE}}\\ {\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^2& {\widehat{\phi }}_0^{\mathrm{MLE}}& 1\end{array}& \begin{array}{ccc}\dots & \dots & {\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^{m-1}\\ \dots & \dots & {\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^{m-2}\\ \dots & \dots & \dots \end{array}\\ \begin{array}{ccc}\dots & \dots & \dots \\ \dots & \dots & \dots \\ {\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^{m-1}& {\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^{m-2}& \dots \end{array}& \begin{array}{ccc}\dots & \dots & \dots \\ \dots & \dots & {\widehat{\phi }}_0^{\mathrm{MLE}}\\ \dots & {\widehat{\phi }}_0^{\mathrm{MLE}}& 1\end{array}\end{array}\right]$$(12)

Using Eq. [11(11) $$P\left(({\widehat{\mu }}_0^{\mathrm{MOM}}-c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\le \underset{¯}{Y}\le ({\widehat{\mu }}_0^{\mathrm{MOM}}+c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\right)=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)=E_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m{\widehat{\phi }}_0^{\mathrm{MLE}}\right)\right)={\int }_{-1}^1\left[\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\dots \underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}{f}_{\mathit{mvn}\left(\underset{¯}{0},\widehat{\underset{¯}{\mathrm{\Sigma }}}\left({\widehat{\phi }}_0^{\mathrm{MLE}}\left(t\right)\right)\right)}\left(x_1,x_2,\dots ,x_m\right)d{x}_{1}d{x}_2\dots d{x}_m\right]f_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(t\right)dt$$(11) ], we set (13) $$1-{\mathrm{FAP}}_0=1-\mathrm{FAP}=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)$$(13) where ${\mathrm{FAP}}_0$ is the nominal false alarm probability, typically taken as 0.05 or 0.1. The Phase I charting constant $c^{(U)}$ is obtained by solving Eq. [13(13) $$1-{\mathrm{FAP}}_0=1-\mathrm{FAP}=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)$$(13) ] numerically. As shown in Eq. [11(11) $$P\left(({\widehat{\mu }}_0^{\mathrm{MOM}}-c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\le \underset{¯}{Y}\le ({\widehat{\mu }}_0^{\mathrm{MOM}}+c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\right)=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)=E_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m{\widehat{\phi }}_0^{\mathrm{MLE}}\right)\right)={\int }_{-1}^1\left[\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\dots \underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}{f}_{\mathit{mvn}\left(\underset{¯}{0},\widehat{\underset{¯}{\mathrm{\Sigma }}}\left({\widehat{\phi }}_0^{\mathrm{MLE}}\left(t\right)\right)\right)}\left(x_1,x_2,\dots ,x_m\right)d{x}_{1}d{x}_2\dots d{x}_m\right]f_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(t\right)dt$$(11) ], the distribution of ${\widehat{\phi }}_0^{\mathrm{MLE}}$ is necessary to obtain the charting constant $c^{(U)}.$ It is well-known that ${\widehat{\phi }}_0^{\mathrm{MLE}},$ as the ML estimator, is asymptotically a normal random variable (see, e.g., Anderson 1994, p. 186). However, we do not consider the asymptotic normal distribution since the required sufficient length of many processes in practice, such as that in the application of the monthly-collected opioid process, may be unavailable, being excessively costly and time-consuming. Instead, the exact distribution of ${\widehat{\phi }}_0^{\mathrm{MLE}}$ is considered in our case. It is well-known that the ML estimators ${\widehat{\mu }}_0^{\mathrm{MLE}},$ ${\widehat{\gamma }}_0^{\mathrm{MLE}},\mathrm{ }$ and ${\widehat{\phi }}_0^{\mathrm{MLE}}$ have no closed analytical form and are solved using some nonlinear optimization methods, given by (see, e.g., Hamilton 1994, p. 122; Anderson 1994, p. 237, Hyndman et al. 2020) (14) $$\left({\widehat{\mu }}_0^{\mathrm{MLE}},{\widehat{\gamma }}_0^{\mathrm{MLE}}=\frac{{\left({\widehat{\sigma }}_0^{\mathrm{MLE}}\right)}^2}{1-{\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^2},{\widehat{\phi }}_0^{\mathrm{MLE}}\right)=\left({\widehat{\mu }}_0^{\mathrm{MLE}},{\left({\widehat{\sigma }}_0^{\mathrm{MLE}}\right)}^2,{\widehat{\phi }}_0^{\mathrm{MLE}}\right)=\mathrm{argmax}\left(\log \hspace{0.17em}\left(L\left({\mu }_0,{\sigma }_0^2\mathrm{ },\mathrm{ }{\phi }_0;\underset{¯}{Y}\right)\right)\right)$$(14) where $L\left({\mu }_0,{\sigma }_0^2, {\phi }_0;\underset{¯}{X}\right)$ is the likelihood function with respect to ${\mu }_0,{\sigma }_0^2, {\phi }_0$ using Eq. [4(4) $$(Y_i-{\mu }_0)={\phi }_0(Y_{i-1}-{\mu }_0)+{\tau }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m,\mathrm{ }\left|{\phi }_0\right|<1\text{ and }{Y}_0={\mu }_0$$(4) ]. A general formulation of the likelihood function can be found in Hyndman and Khandakar (2008), which can be written in the AR(1) case, as $$L\left({\mu }_0,{\sigma }_0^2,\mathrm{ }{\phi }_0;\underset{¯}{Y}\right)=\left(\frac{1}{\sqrt{2\pi {\sigma }_0^2}}{\mathrm{e}}^{-\frac{1}{2}\left(\frac{{(Y_1-{\mu }_0)}^2}{{\sigma }_0^2}\right)}\right)\dots \left(\frac{1}{\sqrt{2\pi {\sigma }_0^2}}{\mathrm{e}}^{-\frac{1}{2}\left(\frac{{\left((Y_m-{\mu }_0)-{\phi }_0(Y_{m-1}-{\mu }_0)\right)}^2}{{\sigma }_0^2}\right)}\right).$$

Recall that our method estimates the ${\mu }_0$ and ${\gamma }_0$ using the MOM and ${\phi }_0$ using the MLE. The exact distribution of ${\widehat{\phi }}_0^{\mathrm{MLE}}$ is approximated using a special bootstrapping method and involves the approximation of Eq. [11(11) $$P\left(({\widehat{\mu }}_0^{\mathrm{MOM}}-c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\le \underset{¯}{Y}\le ({\widehat{\mu }}_0^{\mathrm{MOM}}+c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\right)=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)=E_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m{\widehat{\phi }}_0^{\mathrm{MLE}}\right)\right)={\int }_{-1}^1\left[\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\dots \underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}{f}_{\mathit{mvn}\left(\underset{¯}{0},\widehat{\underset{¯}{\mathrm{\Sigma }}}\left({\widehat{\phi }}_0^{\mathrm{MLE}}\left(t\right)\right)\right)}\left(x_1,x_2,\dots ,x_m\right)d{x}_{1}d{x}_2\dots d{x}_m\right]f_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(t\right)dt$$(11) ]. The steps of the bootstrapping method are divided into two levels as follows

1. The first level of bootstrapping for generating the empirical distributions of the estimators: Generate $S_1$ samples of $\underset{¯}{W}={[W_1,W_2,\dots ,W_m]}^{\prime }$ where $X_i,i=1,2,\dots ,m$ is an i.i.d. normal random variable with mean $0,$ variance $1,$ and obtain $\underset{¯}{\widehat{Y}}={(\widehat{\underset{¯}{\mathrm{\Sigma }}}({\widehat{\phi }}_0^{\mathrm{MLE}}))}^{1/2}\underset{¯}{W}$ where $\widehat{\underset{¯}{\mathrm{\Sigma }}}\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)$ is an estimated covariance matrix using the MLE, as in Eq. [14(14) $$\left({\widehat{\mu }}_0^{\mathrm{MLE}},{\widehat{\gamma }}_0^{\mathrm{MLE}}=\frac{{\left({\widehat{\sigma }}_0^{\mathrm{MLE}}\right)}^2}{1-{\left({\widehat{\phi }}_0^{\mathrm{MLE}}\right)}^2},{\widehat{\phi }}_0^{\mathrm{MLE}}\right)=\left({\widehat{\mu }}_0^{\mathrm{MLE}},{\left({\widehat{\sigma }}_0^{\mathrm{MLE}}\right)}^2,{\widehat{\phi }}_0^{\mathrm{MLE}}\right)=\mathrm{argmax}\left(\log \hspace{0.17em}\left(L\left({\mu }_0,{\sigma }_0^2\mathrm{ },\mathrm{ }{\phi }_0;\underset{¯}{Y}\right)\right)\right)$$(14) ] for a given value of ${\widehat{\phi }}_0^{\mathrm{MLE}}.$

2. From the data generated in step 1, calculate $S_1$ sets of estimates of ${\tilde{\phi }}_{\mathrm{ }0,s_1}^{\mathrm{MLE}}$ using the MLE method.

3. The second level of bootstrapping for generating the empirical distributions of AR(1) Processes: Using each set of estimates obtained in step 2, generate $S_2$ samples as in step 1, denoted as $\underset{¯}{\tilde{Y}},$ and then for these data calculate $S_2$ sets of estimates ${\tilde{\tilde{\mu }}}_{0, s_1,s_2}^{\mathrm{MOM}}$ and ${\tilde{\tilde{\gamma }}}_{0, s_1,\mathrm{ }{s}_2}^{\mathrm{MOM}}$ using the MOM.

4. Estimate the probability of the NSE in Eq. [13(13) $$1-{\mathrm{FAP}}_0=1-\mathrm{FAP}=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)$$(13) ] by (15) $$\frac{1}{S_1}\sum _{s_1=1}^{S_1}\frac{1}{S_2}\left[\sum _{s_2=1}^{S_2}ID\left(\left(-c^{\left(U\right)}\right){\underset{¯}{1}}_m\le \underset{¯}{\tilde{X}}=\frac{\widehat{\underset{¯}{Y}}-{\tilde{\tilde{\mu }}}_{0,\mathrm{ }{s}_1,s_2}^{\mathrm{MOM}}}{\sqrt{{\tilde{\tilde{\gamma }}}_{0,\mathrm{ }{s}_1,\mathrm{ }{s}_2}^{\mathrm{MOM}}}}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)\right]\mathrm{ }$$(15)

where $ID(A)$ is an indicator function, with $ID\left(A\right)=1,$ when the event $A$ holds, otherwise, $ID\left(A\right)=0.$ Note that the required MLE is found using the R function Arima from the R package, forecast in Hyndman et al. (2020). The required charting constants can be found, on demand, using the PH1ARMA function in the accompanying R package (Yao et al. 2021). For ease of implementation, some are tabulated and provided in Table A1 of the Appendix.

In the next section, the charting constants using the proposed method are compared to those using the 3-sigma rule.

3.2. Proposed charting constants and comparison to the classical value 3

This section compares the charting constants using the proposed method to 3 in the classic 3-sigma rule. In the i.i.d. case of Phase I, Yao and Chakraborti (2021b) showed that the classic 3-sigma rule did not perform well. Thus, it is not surprising that the classic rule also does not have a good performance in the autocorrelated case of Phase I. An evidence can be found from the charting constants using the proposed method. It becomes clear that, generally speaking, the classic chart using the 3-sigma rule is not appropriate to set up the Phase I control chart. To explain this aspect in more detail, the relationship among the FAP, the charting constant, the autocorrelation and the sample size is shown and explored visually in the four panels of Figure 1, via the $\mathrm{FAP}$ contours obtained using Eq. [13(13) $$1-{\mathrm{FAP}}_0=1-\mathrm{FAP}=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)$$(13) ], carried out using the accompanying R package.

Figure 1. $\mathrm{FAP}$ contours calculated for different values of $c^{(U)},$ ${\phi }_0$ and $m.$

The horizontal and vertical dashed lines in all panels of Figure 1 represent the i.i.d. case (${\phi }_0=0$) when the charting constant equals 3, respectively. Hence, the intersection of these two lines represents the situation for the classic individual control chart using the 3-sigma rule. Many interesting observations can be made from Figure 1. First, it is seen that the Phase I charting constants for a given ${\mathrm{FAP}}_0$ do not converge to 3 or some other fixed value, as $m$ increases. This can be observed from the $\mathrm{FAP}$ contours in all the panels. For a specific ${\mathrm{FAP}}_0$ value, the contours gradually move to the right (meaning that the charting constant increases) as $m$ increases. Second, and as noted earlier, the 3-sigma chart is not appropriate in Phase I. For this, observe that from panels (a) and (b), it is seen that for $m$ = 10 and 20, the 3-sigma rule corresponds to a much larger ${\mathit{FAP}}_0$ than those from $0.05$ to $0.2$ (which are common choices for ${\mathrm{FAP}}_0$ in practice). From panels (c) and (d), for $m$ = 50 and 100, the 3-sigma rule is seen to work in a narrow range of ${\mathrm{FAP}}_0$ values, from 0.05 to 0.1 and from 0.1 to 0.2, respectively, especially when $\left|{\phi }_0\right|$ closes to zero or it is the i.i.d. case. By comparison, however, the proposed Phase I control chart always maintains the $\mathrm{FAP}$ at the nominal ${\mathrm{FAP}}_0.$ Third, the $\mathrm{FAP}$ contours become more symmetric around ${\phi }_0=0$ as $m$ increases. For example, in panel (d) for $m=100,$ the contours are almost symmetric and the curvature is much lower compared to those for $m=$20 or 50 as seen in panels (b) and (c), respectively. This means that the required charting constants are related to the magnitude (absolute value) of ${\phi }_0$ and the effect of the sign is small when $m\ge 20.$ In summary, Figure 1 shows that in addition to normality, the required charting constant depends on several other factors, such as the magnitude of the autocorrelation, the sample size and the nominal $\mathrm{FAP}$ value and therefore needs to be determined on a case-by-case basis.

The classical constant 3 may work in one or two special cases, where the data are independent, the sample size is medium-large and the nominal $\mathrm{FAP}$ is well-chosen (such as ${\phi }_0=0, m=50$ and ${\mathrm{FAP}}_0=0.1$), but not in general. Thus, the classic 3-sigma chart is not recommended for retrospective monitoring of individual, normally distributed, autocorrelated data from an AR(1) process. As noted earlier, the required charting constants for the proposed Phase I chart can be conveniently found using the accompanying R package in Yao et al. (2021).

Capizzi (2015) noted that the ability of Phase I charts to correctly identify the outliers is “unavoidably related to the correct specification of the underlying IC model. Many Phase I points may fall outside the control limits because the wrong model was assumed for the process variable…” Thus, since both the direct and the residual-based monitoring of autocorrelated data depend on fitting a suitable model to the data, which includes the assumption of a time series process, an important question is the effect of process miss-specification and hence of the model on the performance of the control chart. In the next section, we study the performance differences between the proposed chart and some charts available in the literature in this context using simulated data.

4. Model misspecification analysis: IC robustness and OOC performance

The charts used in the performance comparisons are the classic (I-MR) chart (Montgomery 2005), the residual chart (Alwan and Roberts 1988) and the corrected residual control chart (the AL chart, Apley and Lee 2003, 2008). As noted earlier, these charts are designed for Phase II applications and therefore they need to be adapted here for Phase I to be comparable to the proposed chart, that is, their charting constants need to be obtained for a specific time series process (discussed below) using a Bonferroni adjustment to maintain the FAP. Thus, the adjusted nominal false alarm rate for these charts is set equal to $\frac{{\mathrm{FAP}}_0}{2m}$ and this is used to find the charting constants for these charts as in Champ and Jones (2004). Vasilopoulos and Stamboulis (1978) considered a control chart for AR(1) data, but since their chart is for sub-grouped data, it is not included in the comparison. Another, more recent approach to control charting is the distribution-free approach (see, e.g., Capizzi and Masarotto 2013; Chakraborti and Graham 2019a, 2019b; Jones-Farmer, Jordan, and Champ, 2009; Kim et al. 2007; Li and Qiu 2016; Perry and Wang 2020). However, these distribution-free approaches mainly focus on Phase II monitoring, except for the CM chart which is the only available Phase I control chart for i.i.d. data with an accompanying software. Thus, the CM chart is included in the comparisons; both the IC and OOC cases are considered.

Jones-Farmer, Jordan, and Champ (2009) define two types of OOC conditions, namely an isolated shift and a sustained shift in the process mean. The isolated shift is a mean shift of size $\delta$ in units of the population standard deviation only in the first observation. The sustained shift is a shift in the mean, which is set to start two thirds of the way into the sample and continues until the end. In Phase I, the interest is to detect isolated shifts (like outliers) and examine the corresponding assignable causes, which present short-term unusual disturbance in the process mean. Following Graham, Human, and Chakraborti (2010) and Yao and Chakraborti (2021a, 2021b), we assume that ${\underset{¯}{Y}}_1={[Y_1,\dots ,Y_{t-1},Y_t, Y_{t+1},\dots ,Y_m]}^{\prime }$ is an ARMA process with the isolated shift at the $t$-th observation and, thus, the process is OOC. $Y_t$ is normally distributed with mean ${\mu }_1={\mu }_0+\delta \sqrt{{\gamma }_0},$ where $\delta \ne 0$ is the standard shift size and variance ${\gamma }_0;$ the other observations (those except for the $t$-th observation) follow a normal distribution with mean ${\mu }_0$ and variance ${\gamma }_0.$ The probability of at least one signal ($\mathrm{PAOS}$) is used as an OOC performance metric in Phase I. As in Eq. [13(13) $$1-{\mathrm{FAP}}_0=1-\mathrm{FAP}=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)$$(13) ], this is given by (16) $$\mathrm{PAOS}=\hspace{0.17em}1-P\left(\left({\widehat{\mu }}_0^{\mathrm{MOM}}-c\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}\right){\underset{¯}{1}}_m\le {\underset{¯}{Y}}_1\le \left({\widehat{\mu }}_0^{\mathrm{MOM}}+c\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}\right){\underset{¯}{1}}_m\right)=1-P\left(-c{\underset{¯}{1}}_m\le {\underset{¯}{X}}_1^{\left(U\right)}\le c{\underset{¯}{1}}_m\right)$$(16) where ${\underset{¯}{X}}_1^{\left(U\right)}=\left[X_{1,i}^{\left(U\right)}\right]=\frac{{\underset{¯}{Y}}_1-{\widehat{\mu }}_0^{\mathrm{MOM}}}{\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}}$ is the standardized OOC process. Note that when the process is IC, the $\mathrm{PAOS}$ reduces to the $\mathrm{FAP}$ in Eq. [13(13) $$1-{\mathrm{FAP}}_0=1-\mathrm{FAP}=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)$$(13) ]. Without loss of generality and simplicity, ${\mu }_0$ and ${\sigma }^2$ are set to be 0 and 1, respectively. $m$ is assigned to be 20 and 100 and, correspondingly, $t$ is assigned to be the first observation and the first observation after the first quarter observation of $m.$ That is, $t=1$ and $6$ for $m=20,$ and $t=1$ and $26$ for $m=100.$ Note that other $t$’s including the first observations after the second quarter, the third quarter and the last observations of $m$ have been studied. The results are fairly similar to those mentioned above, so they are ignored. In the performance comparison, various types of time series processes are studied. The AR(1) processes, written as Eq. [4(4) $$(Y_i-{\mu }_0)={\phi }_0(Y_{i-1}-{\mu }_0)+{\tau }_i,\mathrm{ }i=1,2,\mathrm{ }\dots ,\mathrm{ }m,\mathrm{ }\left|{\phi }_0\right|<1\text{ and }{Y}_0={\mu }_0$$(4) ], with $-0.9\le {\phi }_0\le 0.9,$ are considered for the IC situation. In order to study the effects of model miss-specification, the underlying process is assumed to be an MA(1) process, which is written as $\left(Y_i-{\mu }_0\right)=\theta {\tau }_{i-1}+{\tau }_i,$ and the charts are constructed under an AR(1) process assumption. Note that ${\gamma }_0$ is equal to $(1+{\theta }^2){\sigma }^2$ for the MA(1) process, different from that for the AR(1) process, where $\theta \in [-0.9, 0.9]$ is the coefficient for the first-order moving average term. Though some other common processes including an ARMA(1, 1), which is written as $\left(Y_i-{\mu }_0\right)=\phi \left(Y_{i-1}-{\mu }_0\right)+\theta {\tau }_{i-1}+{\tau }_i$ and an AR(2), which is written as $\left(Y_i-{\mu }_0\right)={\phi }_1\left(Y_{i-1}-{\mu }_0\right)+{\phi }_2\left(Y_{i-2}-{\mu }_0\right)+{\tau }_i$ have been studied, the results are provided in Figures A1–A3 of the Appendix. Also, note that the CM chart is used in the comparisons since as claimed by the authors, it can detect “single or multiple mean and/or scale shifts” (Capizzi and Masarotto 2013) and the results are obtained using the rsp function in the R package dfphase1 (Capizzi and Masarotto 2018). The computations for the CM chart are controlled by the maximum number of step shifts and the minimum length of each step, and these are set to 1 and 2, which are the smallest values for the two parameters in the rsp function, respectively, in order to better accommodate the one-isolated-shift situation considered in our simulations. Finally, each parameter combination for different processes generates 1000 realizations in the simulation.

From Figure 2, it can be seen that:

Figure 2. Performance comparison under the AR(1) process when the process goes OOC on the first observation.

1. The proposed charts are IC FAP robust for different ${\phi }_0$ and $m$ values. That is, the simulated $\mathrm{FAP}$ value is close to the ${\mathrm{FAP}}_0$ as shown in the first column of the graphs, in panels (a) and (e), in terms of the differences between the black curves and the dashed lines. The residual charts are also IC FAP robust as seen in these panels, from the difference between the yellow curves and the dashed lines. By comparison, as it might be expected, the Bonferroni adjusted classic I-MR Phase I control chart and the CM chart are seen (in all the panels) to be not IC FAP-robust and hence are not recommended in Phase I applications, when there is autocorrelation. These two charts are mainly affected by ${\phi }_0$ (instead of $\delta$ or $m$) and tend to have more (less) false alarms when ${\phi }_0$ is negative (positive). Note that this phenomenon has also been observed in Maragah and Woodall (1992) for the classic I-MR chart.

2. The AL chart performs worse than the other charts. When $\delta =0,$ the ${\mathrm{FAP}}_0$ of the AL chart equals almost zero, as seen in panels (a) and (e). This is unacceptable. When $\delta >0,$ the PAOS values for the AL chart is much lower than those for the other charts in all panels except for (a) and (e).

3. The proposed chart performs better when the sample size is small, such as $m\le 20.$ The proposed chart and the residual chart using the Bonferroni adjustment both have ${\mathrm{FAP}}_0$ values close to the nominal ${\mathrm{FAP}}_0=0.1$ when $\delta \le 1,$ as seen in panels (a), (b), (e) and (f). For $\delta >1$ and small sample sizes such as $m\le 20$ as seen in panel (c), and considering the $\mathrm{FAP}$ and the $\mathrm{PAOS}$ values, the proposed chart performs better than the residual chart, since the proposed chart (the purple curve) dominates the residual chart (the orange curve). For example, when $m=10$ and $\delta =3,$ the relative changes between the $\mathrm{FAP}$ values of the proposed and the residual chart are $\frac{0.711-0.687}{0.687}=0.0349, \frac{0.442-0.332}{0.332}=0.3313, \frac{0.430-0.310}{0.310}=0.3871,\hspace{0.17em}\frac{0.439-0.317}{0.317}=0.3849,\hspace{0.17em}\frac{0.538-0.404}{0.404}=0.3317, \frac{0.913-0.856}{0.856}=0.0666$ for ${\phi }_0=-0.9, -0.5, -0.1, 0.1, 0.5$ and $0.9,$ respectively. This shows that the proposed chart has higher signaling probabilities for the more common values of ${\phi }_0,$ the performance could be as much as 38% higher.

4. Performance of the proposed chart is similar to that of the residual chart when the sample size is large such as $m>20$ except in some extreme situations where the autocorrelation is very high such as $\left|{\phi }_0\right|=0.9,$ which is close to the non-stationarity bound. For example, when $m=100$ and $\delta =3,$ the relative changes in the $\mathrm{PAOS}$ values between the proposed and the residual charts are $\frac{0.518-0.992}{0.992}=-0.4778, \frac{0.390-0.403}{0.403}=-0.0323, \frac{0.390-0.399}{0.399}=-0.0226, \frac{0.395-0.400}{0.400}=-0.0125, \frac{0.404-0.377}{0.377}=0.0716, \frac{0.628-0.995}{0.995}=-0.3688$ for ${\phi }_0=-0.9, -0.5, -0.1, 0.1, 0.5$ and $0.9,$ respectively.

In summary, the proposed chart is recommended for monitoring individual data from an AR(1) process based on its IC FAP-robustness and OOC performance (in terms of $\mathrm{PAOS}),$ particularly for smaller sample sizes ($m\le 20$), except for some cases of extremely high autocorrelation, such as $\left|{\phi }_0\right|=0.9.$ However, note that for such high values of $\left|{\phi }_0\right|,$ the process is close to being non-stationary, which poses additional modeling and monitoring challenges.

Next, we study and compare the effects of model misspecification on the performance of the charts, designed for the AR(1) process, when the data are simulated using the “miss-specified” MA(1) process. For the MA(1) process, the range $-0.9\le \theta \le 0.9$ is considered for satisfying invertibility. The results are displayed in Figure 3.

Figure 3. Performance comparison under the MA(1) process when the process goes OOC on the first observation.

From Figure 3, it is seen that:

1. The classic, the AL, the residual and the CM charts are all problematic to monitoring the MA(1) data. Recall that in the IC case (when $\delta =0),$ the simulated $\mathrm{FAP}$ of the charts should maintain the nominal ${\mathrm{FAP}}_0=0.1.$ However, the simulated $\mathrm{FAP}$ of the classic, the residual and the CM charts are seen to be increasing as ${\phi }_0$ increases from $-0.9$ to $0.9.$ In other words, for this case, the chart gives fewer false signals for negative values of ${\phi }_0,$ and more false signals for positive values of ${\phi }_0,$ even if $\delta =0.$ For the AL chart, the performance is similar to that in Figure 2; the simulated $\mathrm{FAP}$ values are very close to zero and much less than ${\mathrm{FAP}}_0=0.1.$

2. By contrast, the proposed chart shows a much more reasonable IC performance for the MA(1) data. When $\delta =0,$ the simulated FAP of the proposed chart is close to the nominal ${\mathrm{FAP}}_0=0.1.$ However, the proposed chart is not sensitive to detecting small shifts such as $0<\delta \le 1,$ the simulated $\mathrm{PAOS}$ in the panels (b) and (f) is almost equal to the simulated $\mathrm{FAP}$ of the proposed chart in panels (a) and (e). When $\delta \ge 2,$ the difference between the simulated PAOS and the nominal ${\mathrm{FAP}}_0$ is much higher and the chart is able to detect the shift. Thus, the proposed chart is more capable of detecting larger shifts which is typically the goal in Phase I (Montgomery 2005).

Next, we study the performance of the charts when the process goes OOC on the first observation after the first quarter of observations in Figure 4.

Figure 4. Performance comparison under the AR(1) and MA(1) processes when the process goes OOCs on the first observation after the first quarter of observations.

From Figure 4, it can be seen that all charts perform in a manner similar to what is seen in Figures 2 and 3. In fact, similar performances are observed when the process goes OOC on the first observation after the second, the third quarter of observations and the last observation, respectively, and are therefore not shown here due to the space limitations. These results suggest that the chart performance is not affected by the locations of the OOC observations.

To conserve space, results for the OOC model miss-specification analysis using observations from the AR(2) and the ARMA(1, 1) processes are shown in Section A2 of the Appendix. Similar conclusions can be drawn as those for the MA(1) process discussed above; that is, the proposed chart is (1) sensitive to large shifts in means and (2) fairly FAP-robust for both the AR(2) and ARMA(1, 1) processes. By comparison, the residual charts are not IC FAP-robust. In summary, the classic, the AL and the CM charts are not recommended for monitoring individual autocorrelated data in Phase I. The proposed chart is a safer and better choice in practice from the perspective of IC FAP-robustness since in general there is no guarantee that the model is specified correctly.

5. Illustration

In this section, we discuss an application of the proposed methodology to a real-world dataset in the public health sector where retrospective monitoring is routinely done and yet much of the existing methodology is based primarily on a descriptive analysis, which raises the need for more sophisticated tools. It is well-known that the United States is (has been) experiencing an opioid epidemic. It is in this context of analyzing prescription opioid usage related data, the SPM methodology can add significant capabilities to fight the epidemic, by helping to objectively monitor the data, both in retrospective and prospective phases, which will allow taking appropriate, timely and statistically valid actions by the appropriate authorities. The data for this illustration are the prescription fentanyl (PF) (fentanyl is one of the most well-known opioids) transactional data from the DEA (Drug Enforcement Authority)’s ARCOS database. This is a comprehensive drug reporting system that monitors the supply chain of DEA controlled substances (more information can be found on https://www.deadiversion.usdoj.gov/arcos/index.html) and the selected dataset is a subset of the data from the ARCOS database. It is expected that the proposed Phase I (retrospective) monitoring methodology can be applied to other datasets involving other drugs in similar contexts, which can be a potent and beneficial strategy for effective public health surveillance systems.

Note that the unit of dosage for opioid drugs is reported in milligrams and the strength of the dosage of different opioid drugs is typically transformed to a unified measurement, called the morphine milligram equivalent (MME), using the conversion factor table provided by Centers for Medicare & Medicaid Services (2018). In many opioid-related studies, practitioners are more interested in the MME per capita which is an average of the MME over a population in a specific region such as a county (see, e.g., Guy et al., 2019). We work with the MME per capita transaction data and for illustration, we use the PF transactional data in Preston County, West Virginia (WV), between January 2009 and December 2013. Note that the state of West Virginia has been identified as the the state where the drug overdose is increasing considerably (West Virginia Department of Health and Human Resources, 2021) and the Preston County is selected, since the county is ranked 29 out of the 55 in the vulnerability ranking (West Virginia Department of Health and Human Resources 2020) and thus stands for a typical county in WV in terms of opioid transactions. The dataset is provided in Section A3 of the Appendix along with a time series plot, a histogram and a boxplot in Figure A4.

The solid curve with circles in Figure A4 of the Appendix shows the PF transactional data series (in the MME per capita). From the time series plot, we observe a clear peak around April 2009 and a drop around June 2009. This indicates that the underlying process may have shifted to a different level, which may raise some concerns. From the histogram it is seen that the data distribution is almost symmetric, which is also seen in the boxplot. The latter indicates two outliers, both occurring on the high side, in March and April 2009, respectively. Stationarity of the process is tested using the augmented Dickey–Fuller test provided by the function adf.test in the R package tseries (Trapletti, Hornik, and LeBaron 2019). The value of the test statistic is $-5.1816$ with a p-value less than 0.01, so that the stationarity is not rejected at a significance level 0.05. Continuing the analysis of the data, following the Box–Jenkins time series methodology, we next show the autocorrelation function (ACF) and the partial autocorrelation function (PACF) for the PF transactional data in Figure 5.

Figure 5. ACF and PACF of the PF transactional data in Preston County, WV, from January 2009 to December 2013.

From Figure 5, it is seen that both the ACF and PACF both have a “cut-off” at lag 1. Thus, several models such as the AR(1), MA(1), ARMA(1, 1) and AR(2) can be contemplated under the Box–Jenkins analysis. We check goodness-of-fit for all these models using the Akaike information criterion (AIC), the corrected Akaike information criterion (AICc) and the Bayesian information criterion (BIC). For AR(1), MA(1), ARMA(1, 1) and AR(2) models, the values using the criteria are, (-6.84, −6.41, −0.56), (-6.65, −6.22, −0.36), (-5.37, −4.65, 3.01) and (-5.8, −5.08, 2.57), respectively. Based on these figures (the criteria values are smaller) and the fact that these are real data, the AR(1) appears to be a reasonable first choice. Note that at some other lags, such as at lag 10 in panels (a) and (b), the ACF and the PACF values are close to but do not exceed the dashed lines which are the critical values at the 95% confidence level. This may raise some concerns with respect to the model identification. However, this is not unexpected for a real dataset such as this one and therefore we do not pursue this point further (and do not consider other models for these data here, which could be a topic for future investigations). After fitting the AR(1) model using MLE, the means and the standard deviations for the intercept and the slope for the auto-regressor are $(0.3878, 0.1179)$ and (1.0462, 0.0453), respectively. The normality of the residuals is tested by the Shapiro–Wilk test using the function shapiro.test in the R package stats. This yielded a test statistic 0.9803 with a p-value 0.4448. Hence, no significant evidence against normality is observed and we continue toward finding the proposed Phase I charting constant based on the fitted AR(1) model.

1. Proposed chart: The charting constant is calculated and our software. For ${\mathrm{FAP}}_0=0.05, 0.1, 0.2$ and $m=60,$ the charting constants are found to be 3.1710, 2. 9956, and 2.8082, respectively, so the pairs of lower and upper control limits, corresponding to each ${\mathrm{FAP}}_0$ value, are: (-3.1710, 3.1710), (-2.9956, 2.9956), and (-2.8082, 2.8082). It may be noted that coincidentally, for $m=60,$ the charting constants for the two charts are fairly close, but at other values of $m,$ smaller or larger, the charting constants are much different. Using the sample mean and standard deviation 1.0451 and 0.2381, respectively, process monitoring is done via the standardized variables in Eq. [11(11) $$P\left(({\widehat{\mu }}_0^{\mathrm{MOM}}-c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\le \underset{¯}{Y}\le ({\widehat{\mu }}_0^{\mathrm{MOM}}+c^{\left(U\right)}\sqrt{{\widehat{\gamma }}_0^{\mathrm{MOM}}}){\underset{¯}{1}}_m\right)=P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m\right)=E_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(P\left(-c^{\left(U\right)}{\underset{¯}{1}}_m\le {\underset{¯}{X}}^{\left(U\right)}\le c^{\left(U\right)}{\underset{¯}{1}}_m{\widehat{\phi }}_0^{\mathrm{MLE}}\right)\right)={\int }_{-1}^1\left[\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}\dots \underset{-c^{\left(U\right)}}{\overset{c^{\left(U\right)}}{\int }}{f}_{\mathit{mvn}\left(\underset{¯}{0},\widehat{\underset{¯}{\mathrm{\Sigma }}}\left({\widehat{\phi }}_0^{\mathrm{MLE}}\left(t\right)\right)\right)}\left(x_1,x_2,\dots ,x_m\right)d{x}_{1}d{x}_2\dots d{x}_m\right]f_{{\widehat{\phi }}_0^{\mathrm{MLE}}}\left(t\right)dt$$(11) ].

2. Residual chart: Using the Bonferroni adjustment, for ${\mathrm{FAP}}_0=0.05,\mathrm{ }0.1,\mathrm{ }0.2$ and $m=60,$ the charting constants are found to be 3.1440, 2.9352 and 2.7131, respectively. Also, the value of the mean of moving range with span 2 is 0.2370. Thus, the pairs of lower and upper charting limits with unbiasing constant $d_2=1.128$ are (-0.6606, 0.6606), (-0.6167, 0.6167) and (-0.5701, 0.5701) respectively.

The two control charts are shown in Figure 6 along with the plotted data.

Figure 6. The residual and the proposed control charts for AR(1) model. The solid, dashed and dotted lines are for ${\mathrm{FAP}}_0=0.05,$ 0.1 and 0.2, respectively. In both charts, the charting statistics are shown by open circles and the circles are connected by straight lines for better visualization.

It is seen that at ${\mathrm{FAP}}_0=0.05,$ both charts support the IC status of the process, since no charting statistic plots outside the control limits and no systematic pattern seems to be present. However, the two other pairs of the control limits, for ${\mathrm{FAP}}_0=0.1$ and $0.2,$ respectively, show a different result. Both of these proposed charts signal in April 2009. By comparison, the residual chart shows no signal. To the best of our knowledge, no one has reported this signal in the literature, and this may be an interesting point of investigation in a future study.

When a chart signals, that is a point plots outside the control limits, typically a careful examination is conducted as to any assignable causes and once that is confirmed, the point may be removed and the chart is redone (limits recalculated, remaining points plotted). This process is repeated until the process is declared to be IC (all points plot within the control limits with a random pattern). This is a common strategy for i.i.d. data (Montgomery 2005). To this end, Dasdemir et al. (2016) suggested an iterative way to handle the Phase I OOC observations from an ARMA process. However, the autocorrelation value may change after eliminating the observation(s). In Dasdemir et al., the initial value is 0.387 compared to 0.535 of the final filtered value. This change may lead to a different time series model and hence, a different chart (charting constant) using the proposed method. Also, the effect of the iterative elimination and the imputation of the resulting missing values remain unknown at this point. This will be studied elsewhere.

In summary, the proposed control chart can help in an objective retrospective monitoring (surveillance) of the PF data, which, in turn, can help researchers identify anomalies, make decisions, take corrective actions if any, and generate reference data, to be subsequently used (for model building, parameter estimation, developing control charts, etc.) in prospective (on-line/Phase II) monitoring.

6. Summary and recommendations

In this article, we present a Shewhart-type Phase I control chart for retrospective monitoring of individual autocorrelated data. Such data frequently arise in modern monitoring environments, for example, when data are collected at a high rate of speed by automatic sensors. An AR(1) process is used to model the data with the assumption of normality in the unknown parameter case. The proposed chart accounts for the effects of parameter estimation and autocorrelation and is designed for a nominal IC FAP, a metric typically recommended for Phase I applications to control the overall false alarm rate. We derive the required charting constants, and some tabulation is presented in Table A1 of the Appendix. The accompanying R package PH1XBAR will facilitate the on-demand deployment of the proposed control chart in practice. A simulation study is conducted for monitoring a number of different time series processes in the ARMA family. Findings strongly suggest that it is ill-advised to use a Phase I Shewhart-type chart for i.i.d. data in case of autocorrelated data because of the risk of increased or decreased false alarms. This may cause a disruption in process monitoring, which may result in an inefficient use of resources to situations that are actually in control. The autocorrelation, if it exists in the data, cannot be ignored and must be accounted for in the determination of the appropriate control limits along with the issue of estimated parameters. Also, the proposed chart is more robust to the ARMA processes compared to those in the literature.

Given the normality assumption which needs to be thoroughly examined before deployment, some adaptations may be necessary for practical applications. For example, the errors may be non-normal in which case a “quick and dirty” solution will be to transform the data using say the Wilson–Hilferty (Wilson and Hilferty 1931) transformation or the Box–Cox transformations (Box and Cox 1964) and then applying the proposed charts on the transformed data. In general, one could also develop a model for non-normal errors such as some skewed parametric distribution, like the lognormal, or, use a distribution-free method to build a control chart. As previously noted, although some distribution-free methodology, including control charts, have been considered in the literature, the autocorrelated data case has not been studied. We believe these are all interesting topics to be studied in the future.

Also, note that other common attributes of time series data that may be helpful in the modeling process, such as seasonality, non-stationarity and additional regressors (co-variates) are not accounted for in the proposed methodology. Thus, if these are present, users may need suitable adaptations such as appropriate differencing (to deal with non-stationarity and seasonality), which will be incorporated in a future update of the R package. Finally, the development of the online counterpart for this methodology, useful for Phase II surveillance, is in progress and will be reported elsewhere.

Acknowledgments

The authors thank two anonymous reviewers and the editor for their encouragement and comments that have improved the article.

Disclosure statement

No potential conflict of interest was reported by the authors.