Spatial MCUSUM Control Chart

This paper proposes a spatial multivariate CUSUM control chart in order to monitor the mean of a single characteristic of a product or process, when the measurements are taken in different locations on each sampled item. To estimate the variance and covariance matrix some tools from the geostatistics are used, taking into account the spatial correlation between the measurements. The performance of this control chart is explored by simulation and its use is illustrated with an example.


Introduction
Statistical process control (SPC) is a set of tools for monitoring a production process, useful to achieve the stability of the process and reduce its variability, a Statistician. E-mail: judrojasgo@unal.edu.co b PhD. E-mail: rdguevarag@unal.edu.co among which are included the control charts. The processes always have variation, which can be inherent or natural to the characteristics of the process, called common cause variation; or obey special circumstances that are not permanent in the process, called assignable (or special) cause variation. When a process is operating in the presence of special causes, it is considered to be unstable or out of statistical control (OC); otherwise, the process is considered to be stable or in statistical control (IC). A fundamental objective of a control chart is to analyze the variability of the process and quickly detect assignable causes when they occur. Readers interested in this topic can consult the following textbooks: Ryan (2011), Qiu (2013) and Montgomery (2019), among others.
Generally, the statistical process control can be divided in two phases called phase I or retrospective phase and phase II or prospective phase. In phase I a set historical process data is used to learn about the process. Among the main goals of phase I are: to understand the sources of process variability, to evaluate the process stability searching unusual results and removing assignable causes of variation, and to estimate the in-control process parameters. The major goal of the phase II is to detect any change from the assumed in-control model. Technical details, comments on some of the important aspects of phase I analysis and review about this area can be found in Chakraborti, Human & Graham (2008) and Jones-Farmer, Woodall, Steiner & Champ (2014).
In industry there are many processes characterized for two or more quality characteristics, generally correlated. Several multivariate control charts have been proposed for monitoring these processes, among which are the proposals elaborated by Hotteling (1947), Alt & Smith (1988), Crosier (1988), Pignatiello Jr & Runger (1990), Lowry, Woodall, Champ & Rigdon (1992), Khoo, Wu, Castagliola & Lee (2013), (Lee, Khoo & Xie 2014) and Nezhad & Niaki (2013). The traditional multivariate control charts ignore the correlation between adjacent sample locations, a situation common in some modern manufacturing processes. For example, the problem of monitoring at different locations bottle-wall thickness used for soft drinks and other beverages (Grimshaw, Blades & Miles 2013), or the identification of local distortions surrounding the cylinder bores in a deck face on an automotive engine head, which impact the sealing performance of engine assemblies and reflect cutting force dynamics during machining processes, Suriano, Wang, Shao, Hu & Sekhar (2015). Some proposals that consider the correlation of adjacentes sites have been recently developed by Wang, Jiang & Li (2016) and Zhang, Liu & Jung (2019). Some recent review papers about multivariate control charts are Woodall & Montgomery (2014), Capizzi (2015), Bersimis, Sgora & Psarakis (2018), and Peres & Fogliatto (2018).
Most multivariate control charts seek to monitor multiple quality characteristics simultaneously, however Grimshaw et al. (2013) propose a spatial T 2 control chart and a multivariate EWMA chart (MEWMA) in order to monitor the mean of spatial data (a single quality characteristic at different spatial locations), that is, instead of measuring multiple quality characteristics, only one quality characteristic is measured in multiple sites of the product or process. These control charts are based on multivariate control charts where the sample covariance matrix is replaced with a parametrized covariance matrix based on the semivariogram that models the spatial correlation, under the assumption that the process is isotropic stationary Gaussian. While the proposals presented by Grimshaw et al. (2013) monitor the mean of a spatial process, Garthoff & Otto (2015) recently considered the problem of simultaneously monitoring the mean and covariance of multivariate spatial processes using conventional MEWMA and multivariate cumulative sums (MCUSUM) control charts. It is important to highlight that the proposals presented by Grimshaw et al. (2013) deal with the detection of changes in the mean of the spatial process over time, while the proposal presented by Garthoff & Otto (2015) searches to detect changes in the mean and covariance of the spatial process in terms of the distance from a known centre of the spatial process; therefore, the run length is not anymore a point of time but a distance in space.
In many manufacturing processes the quality of a process or product is best described by a function, called a profile, this is, a functional relationship between a variable and a set of explanatory variables. In these cases, the point of interest are the changes in the profile over time. Profile monitoring is an example of multivariate monitoring and it also includes the monitoring of two-dimensional shapes and three-dimensional surfaces, Woodall & Montgomery (2014). Some control charts proposed for monitoring processes characterized for profiles are based on multivariate control charts, among which are the approaches proposed by Kang & Albin (2000) and Williams, Woodall & Birch (2007). Among the control charts proposed for monitoring surface that consider the spatial correlation are the proposals made by Wang, Wang & Tsung (2014) and Colosimo, Cicorella, Pacella & Blaco (2014). Wang et al. (2014) proposed a control chart in phase II, that describes the spatial correlations among variables in 2-D surface data. This chart is based on the Gaussian-Kriging model, in which the spatial correlations within the 2-D surface profile are represented by a parametric function. The proposal presented by Colosimo et al. (2014) models the manufactured surface via Gaussian processes models and monitoring the deviations of the actual surface from the target pattern estimated in phase I. In this sense, the control charts proposed by Grimshaw et al. (2013) may be extendend for monitoring profiles with spatial dependence. Readers interested in reading about profile monitoring can consult Woodall, Spitzner, Montgomery & Gupta (2004), Woodall (2007), and Noorossana, Saghaei & Amiri (2011), among others.
The control charts proposed by Grimshaw et al. (2013), based on the T 2 and multivariate EWMA control charts, allow to monitor the mean of a spatial process. The first is designed to detect great shifts in the mean of the spatial process while the last is designed to detect smaller shifts. Although Grimshaw et al. (2013) proposed these two new control charts taking account the spatial dependence, only the performance of the control chart based on the T 2 statistic is evaluated in the simulations and illustrated in the example. Other well-known and effective alternative to the multivariate EWMA control chart for detecting small sustained shifts in the mean vector is the multivariate cumulative sum control chart, Montgomery (2019). According to our knowledge, there is no MCUSUM control chart in the literature to monitor the mean of a spatial process. In this paper, we propose an extension of a multivariate CUSUM chart that takes into account the spatial dependence between the measurements of a single quality characteristic in different product locations, where the covariance matrix is determined by the semivariogram that models the spatial correlation. This paper is organized as follows: Section 2 presents some basic concepts of spatial statistics. Section 3 describes the control charts proposed by Grimshaw et al. (2013) for monitoring the mean of a spatial process. In Section 4, the spatial multivariate CUSUM control chart is proposed. Section 5 presents a simulation study to evaluate the performance of the spatial T 2 , EWMA charts and the new spatial MCUSUM control chart, to detect changes in the mean of a spatial process in phase II. An example that illustrates the use of the control chart proposed is given in Section 6. Finally, a summary of results, conclusions and some recommendations are presented in Section 7.

Spatial Data Analysis
Spatial data analysis is the set of statistical techniques for the proper analysis of spatial data, which makes use of the spatial referencing and Geographical Information Systems (GIS) techniques and methods. In this context, the spatial dependence is considered a fundamental property (Haining 2003), where spatially near observations tend to be more similar than distant observations, property known as Tobler's law of geography (Tobler 1970). In this sense, the strength of the relationships between spatially distributed quantities is a function of their spatial separation, (Schabenberger & Pierce 2001).
Let s ∈ D be a generic location data in d-dimensional Euclidean space and suppose the potential datum Z(s) at spatial location s is a random quantity, where s varies over index set D ⊂ R d generating the multivariate random field {Z(s) : s ∈ D}, with Z(s) = [Z(s) 1 , . . . , Z(s p )] ′ . The characteristics of the index set D determine the type of spatial data: Geostatistical data, if D is continuous and fixed; lattice data, if D is discrete and fixed; and point patterns, if the set of locations is itself random, (Cressie 2015).
Stationarity means that the random fields looks similar in different parts of the domain D, (Schabenberger & Pierce 2001). There are three types of stationarity: strict, second-order and intrinsic. The process is strict (or strong) stationarity if the spatial distribution is invariant under translation of the coordinates by the vector h, that is, a. The expected value (mean of the random field) exists and does not depend The covariance function is finite and only depends on the distance between the locations, that is, Cov(Z(s), Z(s + h)) = C(h). The function C(·) is called covariogram or covariance function.
c. The variance function is defined by the literal b when h = 0, therefore V ar(Z(s)) = C(0) = σ 2 < ∞ does not depend on the location.

The process {Z(s) : s ∈ D ⊂ R d } is called intrinsically stationary if E[Z(s)] = µ and V ar(Z(s) − Z(s + h)) = 2γ(h). The function γ(h)
is called the semivariogram of the spatial process.
If a random field is strictly stationary it is also second-order stationary, but the reverse is not necessarily true. A Gaussian random field is defined as a random function whose finite-dimensional distributions are multivariate Gaussian. If a Gaussian random field is second-order stationary it is also strictly stationary, (Schabenberger & Pierce 2001).
A second-order stationary process is also intrinsically stationary, but intrinsic stationarity does not imply second-order stationarity, (Schabenberger & Gotway 2017). Second-order stationarity implies the following relationship between the semivariogram function and the covariogram function: (1) Whereas in a stationary random field the absolute coordinate differences are immaterial, the orientation (angle) of the lag vector h matters. If the semivariogram or covariance function do not depend on the direction of the lag vector h, but these depends only on the absolute distance between points, the random field is termed isotropic, (Schabenberger & Gotway 2017).
In a second-order isotropy stationary random field with isotropic covariogram the covariance between any two points Z(s) and Z(s + h) is only a function of the Euclidean distance between the two points, Cov(Z(s), Z(s + h)) = C(∥h∥), Schabenberger & Pierce (2001).
For a second-order stationary random field, the isotropic semivariogram has a typical form showed in the Figure 1, where one can observe the following elements: nugget effect, sill and range. The nugget effect, denoted by τ 2 , represents the discontinuity at the origin of the semivariogram, which may be due to the fact that part of the spatial structure is concentrated at distances less than observed. The sill, denoted by σ 2 , corresponds to the upper asymptote of the semivariogram and it represents the V ar[Z(s)] = C(0). The range, denoted by ϕ, refers to the distance from which two observations are spatially uncorrelated.
The semivariogram is a function of the spatial process and as such satisfies certain properties, among these, the semivariogram must be conditionally negativedefinite, this is  The goal of semivariogram estimation is to estimate the unknown parameters of a theoretical semivariogram model γ(h, θ). Estimation of the semivariogram is not easy because it must be conditionally negative definite. Estimating a semivariogram is usually a two-step process: (i) to derive an empirical estimate of the semivariogram from the data and (ii) to fit a theoretical semivariogram model to the empirical estimate, Schabenberger & Pierce (2001).
The semivariogram conveys information of the spatial dependence of a random field. An unbiased estimator of the semivariogram, called the classical semivariogram estimator or Matheron estimator, is given bŷ where N (h) is the set of locations pairs that are separated by the lag vector h and |N (h)| denotes the cardinality of this set.
Estimation of the semivariogram by parametric statistical methods requires the selection of a semivariogram model γ(h, θ), where θ is a vector of parameters that is estimated from the data by direct or indirect methods (methods that process the data in some form generating summary amounts, and then fit the semivariogram model to these summaries). Functions that serve as semivariogram models must be conditionally negative definite, Schabenberger & Pierce (2001). Among the semivariogram models more common we found the linear, spherical, exponential, Matern and Gaussian model. In particular, the exponential model is given by For the exponential semivariogram to be valid we need to have τ 2 ≥ 0, σ 2 ≥ 0, and ϕ ≥ 0. The covariance function of the exponential model without nugget effect is The more flexible Matern semivariogram model is given by where Γ(·) is the gamma function and K k (·) is the modified Bessel function of order k. This function is valid for ϕ > 0 and k > 0. The case when k = 0.5 reducing to the exponential model.
The readers interested in explore more about spatial data analysis can consult (Schabenberger & Pierce 2001), (Haining 2003), (Cressie 2015), and (Schabenberger & Gotway 2017) among other textbooks. Grimshaw et al. (2013) proposed two control charts for monitoring the mean of a spatial process, based on the classical multivariate T 2 control chart and MEWMA control chart. In these control charts instead of measuring multiple quality characteristics, only one quality characteristic is measured in multiple sites of the product or process, where the spatial dependendency sctructure of the process is represented by the parametrized covariance matrix based on the semivariogram that models the spatial correlation, under the assumption that the process is isotropic stationary Gaussian.

Control Charts for the Mean of a Spatial Process
Let Z j,t = [Z j,t (s 1 ), . . . , Z j,t (s p )] ′ be a p × 1 vector with the measurements of the quality characteristic on the jth product, j = 1, . . . , n, at the time t, from the spatial process {Z(s) : s ∈ D}. When the process is in control, µ(s) = µ 0 , for all s, but when the process is out of control, µ(s) ̸ = µ 0 for at least one location s. Let Z t be mean vector of the n sampled products at each of the p locations, at the time t. Under the assumption that Z j,t is an isotropic stationary Gaussian process with mean E(Z(s)) = µ(s) = µ 0 , the T 2 control chart produces an out-control signal when exceeds the upper control limit χ 2 (1/ARL 0 ; p), where ARL 0 denotes the incontrol ARL, µ 0 = [µ 0 , . . . , µ 0 ] ′ is the in-control mean vector and Σ(θ 0 ) is the covariance matrix based on the semivariogram that models the spatial correlation with parameters θ 0 .
The spatial multivariate EWMA control chart proposed by Grimshaw et al. (2013) detects small and medium shifts, sustained over time, in the mean of a spatial process, based on the classical MEWMA control chart. This new chart signals a process is out of control at time t when exceeds the upper control limit, where with 0 < λ ≤ 1 and where Σ(θ 0 ) is constructed from the semivariogram. Grimshaw et al. (2013) proposed two control charts for monitoring the mean of a spatial process based on the classical T 2 and MEWMA control charts. The spatial control chart is a straighforward application of the T 2 control chart, therefore this control chart is relatively insensitive to small and moderate shifts in the mean vector because it uses information only from the current sample. The spatial control chart can be used in both phase I and phase II situations. The spatial MEWMA control chart is a modified version of the classical MEWMA chart, was designed to detect small and medium shifts that persist over time, and as in the MEWMA chart, this version is a phase II procedure.

Spatial multivariate CUSUM Control Chart
Other well-known and effective alternative to the multivariate EWMA chart for detecting small sustained shifts in the mean is the multivariate cumulative sum chart, Montgomery (2019). In this paper, we propose an extension of a multivariate CUSUM chart that takes into account the spatial dependence between the measurements of a single quality characteristic in different product locations for monitoring the mean of a spatial process. In this new chart, under the assumption that the process is isotropic stationary Gaussian, the spatial dependendency sctructure of the process is represented by the parametrized covariance matrix based on the semivariogram that models the spatial correlation. We will call this new chart as spatial MCUSUM, which is based on the classical MCUSUM proposed by Crosier (1988).
Let Z j,t (s) = [Z j,t (s 1 ), . . . , Z j,t (s p )] ′ , j = 1, . . . , n, be an vector in phase II, at a time t, with the measurements on the jth product of the quality characteristic of interest taken at the locations s 1 , . . . , s p . We assume that Z j,t is an isotropic stationary Gaussian process with mean vector µ 0 = [µ 0 (s 1 ), . . . , µ 0 (s p )] ′ , E(Z j (s i )) = µ 0 (s i ) = µ 0 , and covariance matrix determined by the semivariogram that models the spatial correlation, Σ(θ 0 ). The in-control mean vector µ 0 describes the mean of a single quality characteristic measured at different locations. The statistic plotted on the control chart at that point of time is where S 0 = 0,Z t = 1 n n ∑ j=1 Z j,t (s) is the mean vector of the n sample products at each of the p locations at the time t, k = 0.5 is an allowance constant, and The chart gives a signal of spatial process mean shift at the t-th time point when U t > h, where h > 0 is the control limit chosen to reach a pre-specified ARL 0 value, Σ(θ 0 ) is constructed from the semivariogram, and θ 0 contains the parameters of the semivariogram of a second-order stationary random field, for example, θ 0 = (σ 2 , ϕ) for an exponential semivariogram model without nugget effect. When C t is less than or equal to k, S t is set to be 0, because there is little evidence of process mean shift in such cases, Qiu (2013).
The spatial MCUSUM control chart is based on the assumption that the process is isotropic stationary Gaussian, therefore the assumptions of normality, constant variance at all s, and isotropy must be guaranteed. Qiu (2008) and Qiu (2013) showed that the actual IC ARL values of the traditional MCUSUM chart are far away from their nominal IC ARL value for some non-normal distributions. If the number of locations is larger than the number of observations, (p ≫ n), the central limit theorem is not applicable, therefore the distribution of the sample mean vector may be multivariate non-normal when the observations from a multivariate non-normal process. Therefore, it is important to ensure that the process is Gaussian. Initially, we would test the multivariate normality for the finite-dimensional distribution at the specific locations where the quality characteristic is measured, using some tests such as Mardia's, Henze-Zirkler's, or Royston's test or some graphical approaches such as chi-square Q-Q plot; but this does not guarantee that the spatial process is Gaussian. Some tests have been designed to determine whether a spatial process is Gaussian. Yuan (2000) derived tests of Gaussianity for stationary random fields on Z d (d ≥ 1) based on the estimated kth-order spectrum (k > 2) at all the discrete frequencies. He used the likelihood-ratio principle. Di Bernardino, Estrade, León et al. (2017) dealt with a real-valued stationary isotropic random field and used the information given by level functionals of a single realization of the process to build a test of Gaussianity. The level functional is the Euler characteristic (EC) of the excursion sets above some levels. More research is required to study the robustness of the proposed spatial MCUSUM chart, as well as the development of new robust control charts to monitor the mean of a single characteristic of a product or process, when the measurements are taken in different locations on each sampled item.
An user may use graphical diagnostics, such as directional sample variograms, rose diagrams, and empirical semivariogram contour plot to decide whether the assumption of isotropy is reasonable. These graphical techniques can be difficult to assess, open to subjective interpretations, and misleading. Hypothesis tests of the assumption of isotropy may be more objective. To this end, we found among others, the proposals formulated by Cabaña (1987), Maity & Sherman (2012), and Weller & Hoeting (2020).
There are only few papers that discuss hypothesis tests of stationarity, among these, Fuentes (2005) presented a test of stationarity for spatial random fields based on spatial spectral analysis (spectral functions which are space dependent), under the assumption that the spatial random field is on a regular grid. The proposed method consists essentially in testing the homogeneity of a set of spatial spectra evaluated at different locations.
Because the assumption of isotropic stationary spatial process, when the process is in control, the mean of the quality characteristic that describes the spatial process is constant for all location. If the spatial MCUSUM chart exhibits an out-of-control signal, at least one of the physical locations departs from the specifications. Because the spatial correlation, using kriging the users can predict the quality charateristic at locations where measurements were not taken. Therefore, the sample unit can be diagnosed with a graphic representing measurements at all locations, and thus detect defective regions.
Because the assumption of stationarity of the spatial process and to the fact that the in-control covariance matrix is a function of the parameters of the semivariogram and the distances for each pair of points, the spatial MCUSUM control chart can be computed even if the spatial measurement locations and/or the number of spatial measurements p change.
The efficient estimation of semivariogram parameters is affected by the selected sampling design, this is, the number of sites and the spatial configuration where measurements are taken. Grimshaw et al. (2013) pointed out, the optimal allocation scheme depends on the type of defects that are of concern: the best performance is obtained when the measurement grid is chosen to cover the locations and types of anticipated defects and then the sample size maximized subject to resource constraints. In this paper, we do not study the effect that the implementation of some sample design criteria would have on the performance of the proposed control chart. It would be interesting to investigate this issue. Some design criteria, among others, have been proposed by Bogaert & Russo (1999), Müller & Zimmerman (1999), Marchant & Lark (2007). An interesting review about this issue is presented by Zimmerman & Buckland (2019).
In the spatial MCUSUM chart the covariance matrix is determined by the semivariogram that models the spatial correlation. The semivariogram is a parameter of the spatial process that generally must be estimated from the data in phase I, usually in a procedure of two steps: (1) derive an empirical estimate of the semivariogram from the data and (2) fit a theoretical semivariogram model to the empirical estimate. There are various methods of fitting semivariogram models, such as least squares, maximum likelihood, and robust methods (Cressie 2015). These techniques require a large number of variogram points. Another alternative is a visual fit of the variogram points to a few standard models. A visual check against a fitted theoretical model is appropriate (Hohn 1998). Studying both the impact of the misspecifying the semivariogram model and the effect of estimating the parameters of the semivariogram on the performance of the spatial MCUSUM chart in phase II are issues of great importance that are not addressed in this paper, which deserve further investigation.

Simulations
In this section, we present some simulation results to evaluate, in phase II, the numerical performance of the proposed spatial MCUSUM chart. We have organized these results in three subsections: in the first, the spatial T 2 , MCUSUM and MEWMA control charts are compared and the effects of sample sizes and magnitude of the non-centrality parameter are analyzed; in the second, the effects of the range of the semivariogram model in the performance of the spatial MCUSUM chart are evaluated, and finally the third subsection evaluates the robustness of the performance of the spatial MCUSUM control chart under deviations in the normality assumption.
The simulations are carried out under the hypothesis that the spatial correlations come from a exponential model without a nugget effect with parameters σ 2 = 0.03 and ϕ = 1. A spatial measurement grid density of 12 points is considered, with 1.5 units between points.

Performance Comparison Between the Spatial T 2 , MEWMA and MCUSUM Control Charts
We compare the numerical performance of the proposed spatial MCUSUM chart with the performance of the spatial T 2 and spatial MEWMA (considering λ = 0.1, 0.2, 0.4, 0.8) control charts. In order to measure the numerical performance of the proposed control chart in phase II, we calculate the average run length for different changes in the mean vector, (ARL 1 ). The level of shifts in the mean vector is described by the non-centrality parameter (ncp), which is defined as where µ 0 represents a vector of means under control and µ 1 the new vector of means, (Montgomery 2019). Large values of δ represent large shifts in the mean vector. When δ = 0, the mean process is under control and the measure of performance is denoted as ARL 0 .
Two sample sizes are explored: n 1 = 1 and n 6 = 6. For the three spatial control charts considerated, the control limits and ARL are calculated using the next algorithm: 1. An initial control limit (h) is established.
2. Under the assumption that the spatial process is isotropic, stationary and Gaussian, random vectors from a multivariate normal distribution are simulated with mean vector zero (µ 0 = 0) and covariance matrix defined by an exponential semivariogram model (Σ 0 ) until the control statistic exceeds the value of the control limit h.
3. The run length (RL) is calculated.
4. Steps 2 and 3 are repeated 100.000 times and the ARL is estimated.
, h is decreased) and go to step 2.
7. 100.000 random vectors are generated from a multivariate normal distribution with mean vector µ 1 and covariance matrix Σ 0 .
The R code for implementing the proposed scheme, including the procedure for finding the control limit, is available from the authors upon request.  Table 1 shows the performance of the spatial T 2 , MEWMA (with λ = 0.1, 0.2, 0.4, 0.8) and MCUSUM control charts for processes with spatial correlation, where the spatial correlation from an exponential model without nugget effect with σ 2 = 0.03 and ϕ = 1. The performace of these control charts is calculated for 9 values of δ. For small shifts in the mean vector (δ=0.06 to δ=0.312) the spatial MCUSUM chart has better performance than the spatial T 2 and MEWMA control charts. For moderate shifts in the mean vector (δ = 0.696, 1 and 2.02) the spatial MEWMA(λ = 0.1) presents better performance. For large shifts (δ = 3), the spatial T 2 and MEWMA control charts show a better performance when n = 6, while the MEWMA(λ = 0.1) is better for n = 1. Table 2 presents a comparison in the performance of the spatial MCUSUM chart for n = 6 and different values of range (ϕ) in the exponential model considered. According to the results, the chart presents a similar behavior for different values of the range, that is, the spatial MCUSUM chart is not sensitive to changes presented by the parameter ϕ. Table 2: Performance comparison of the spatial MCUSUM chart for n = 6 and processes where the spatial correlation from an exponential model without a nugget effect with parameters σ 2 = 0.03 and different values of ϕ.  Table 3 shows the influence of the sample size in the performance of the spatial MCUSUM chart under two range values and two sample sizes. For the two range values, the performance of the spatial MCUSUM chart is better when the sample size is larger. Table 3: Influence of the sample size in the performance of the spatial MCUSUM chart for n = 1 and 6 from an exponential model without a nugget effect with parameters σ 2 = 0.03 and ϕ = 1 and 4.

Robustness of the Performance of the Spatial MCUSUM Control Chart, Under Deviations in the Normality Assumption
Both the multivariate t and the multivariate normal are members of the general family of elliptically symmetric distributions. Multivariate t distributions are generalizations of the classical univariate Student t distribution. Multivariate t-distributions offer a more viable alternative with respect to real-world data, particularly because its tails are more realistic than the multivariate normal distribution. This distribution has been used in different fields such as physics, engineering and finance. Its degrees of freedom parameter, υ, is also referred to as the shape parameter, because the peakedness of the density function may be diminished, preserved, or increased by varying υ (Nadarajah & Kotz 2005). The covariance matrix of this distribution is only defined if υ > 2 (Frey 2010). Other family of distributions with applications in the real world is the Multivariate generalized hyperbolic distributions (MGH distributions). The MGH distributions were introduced by Barndorff-Nielsen (1978). These distributions in general are skewed and represent an attractive family of distributions with exponentially decreasing tails for multivariate data modelling (Frey 2010). Under the exponential model considered and n = 6, Table 4 shows the performance of the spatial MCUSUM chart under three different distributions: multivariate normal, multivariate t-student with υ = 3 and 10, and M GH(0.5, 0.5, 2, µ 0 , Σ(θ 0 ), 0). Table 4: Effect of the distributional assumption in the performance of the spatial MCUSUM chart δ N (µ 0 , Σ(θ 0 ) t 3 (µ 0 , Σ(θ 0 )) t 10 (µ 0 , Σ(θ 0 )) M GH(0.5, 0.5, 2, µ 0 , Σ(θ 0 )), 0) 0 According to the obtained results, the performance of the spatial MCUSUM control chart is sensitive to the distributional assumption. However, because the control limits are found by simulation to guarantee a predetermined ARL incontrol, no distributional assumptions are required for the data.

Example
To illustrate the application of the proposed spatial MCUSUM control chart, we consider the problem of monitoring bottle thickness presented by (Grimshaw et al. 2013). The production of bottles is done by a blown model, which can vary the thickness of the bottle-wall when the air pressure is not uniform in the time, causing the surface of the bottle to curve inward. Bottle-wall thickness is measured using a nondestructive ultrasound micrometer. Bottle is a three-dimensional object, which is represented more simply in two dimensions "unwrapping" the bottle.
In this example, according to the requirements of the manufacturer, 8 locations of the bottles were selected (p), which belongs the phase I sample, with a uniform thickness of 0.055 inch. The sample is made up 26 bottles, includes four locations around the cylinder at two heights, p = 8. There were m = 4 days where n = 6 bottles were measured in a thick spatial measurement grid, to improve production of semivariogram for distances close to zero were collected in two other denser spatial measurement grids, a 5 × 5 grid with dimensions of 0.25 inches away one day and a 5 × 5 grid with dimensions of 0.125 inches away the other day. These measurements are skewed, but the log transformation produces data that are approximately distributed as multivariate normal, according to multivariate version of the statistic Shapiro-Wilk (p-value= 0.6832).
The in-control mean vector is µ 0 = ln(0.055)1 p = −2.91 p , where 1 p is an p × 1 column vector of ones. The in-control covariance matrix is estimated based on the semivariogram that models the spatial correlation. Figure 2 shows the plot of (Z(si)−Z(sj )) 2 2 against the arc distance measured in inches, d(s i , s j ), for each pair of points on the bottles from the phase I sample (first 4 days), represented by circles, which indicate that there is no evidence of a nugget effect. The figure also includes the estimated semivariogram based on the exponential model without nugget effect (solid line), the Matern model without a nugget efect (black dash-dot line) and the Spherical model (gray dotted line). The parameters of the exponential semivariogram model without a nugget effect are estimated using n-weighted least squares and a numerical minimization algorithm in the geoR package of software R, (Ribeiro Jr. & Diggle 2001). Notice in Figure 2 that small d have smaller error bars than larger d. Cressie (2015) describes how a simple nonlinear least-squares estimator would ignore the fact that the variance is a function of distance and a weighted least-squares approach could improve model fit by allowing different variances. Of the three spatial correlation models the exponential model better fit. For the exponential modelσ 2 = 0.034 andφ = 3.194; for the Matern model σ 2 = 0.028,φ = 0.167 and k = 52.791, and for the spherical modelσ 2 = 0.028 and ϕ = 5.
There is no evidence of a nugget effect because the empirical semivariogram cloud approaches zero for small distances. Based on these phase I sample bottles, the spatial control chart will use the covariance based on the exponential model without a nugget-effect. Figure 3 shows the spatial T 2 , MEWMA and MCUSUM control charts using the exponential, Matern and spherical models to estimate the semivariogram, where the upper control limit for the spatial T 2 is calculated as the percentile χ 2 0.005,8 = 21.95, the control limits for the spatial MCUSUM chart and spatial MEWMA (λ = 0.2) chart are h = 12.7 and h = 27.2, respectively, which are chosen by simulations to obtain an ARL 0 = 200. The spatial MCUSUM chart behaviour is similar in the three models; the same situation occurs for the spatial T 2 and MEWMA(λ = 0.2) chart. These control charts show 2 out-of-control signals, corresponding to the days 7 and 8.  In this example, the three control charts show the same performance, this may be because only 4 days are used for the estimation of parameters in phase I.

Conclusions
This paper presents a new type of MCUSUM control chart for monitoring, in phase II, the mean of a single characteristic with spatial dependence. Under the assumption that the observations from an isotropic stationary Gaussian process, the control chart proposed taking account the spatial correlation that exists between the observations of the characteristic that are measured at different sites. This new control chart is called spatial MCUSUM. To ensure a good performance of the proposed control chart, the assumptions of multivaraite normality, constant variance at all location, and isotropy must be verified.
Although the spatial MCUSUM chart constitutes an extension of the classical multivariate cusum control chart, the control chart proposed differs in three fundamental points: only a single quality characteristic is measured in different sites, the covariance matrix is estimated based on the semivariogram that takes into account the spatial dependence, and the observations from an isotropic Gaussian process.
According to the simulation study, which considered only the exponential semivariogram, for small shifts in the mean vector the spatial MCUSUM chart has better performance than the spatial T 2 and MEWMA control charts. For moderate shifts in the mean vector the spatial MEWMA(λ = 0.1) presents better performance. For large shifts the spatial T 2 (with n > 1) and MEWMA control charts show a better performance. It would be interesting to explore the performance of these three control charts using other semivariogram models. The performance of the spatial MCUSUM control chart is sensitive to the distributional assumption, although the control limit can be found by simulation to guarantee a predetermined ARL in-control. More research is required to study the robustness of the proposed spatial MCUSUM chart.
The spatial MCUSUM control chart is based on the assumption that the process is isotropic stationary Gaussian, therefore is important verify, in phase I, that the spatial process satisfies the assumptions of isotropy, stationarity and normality. An user may use graphical diagnostics, such as directional sample variograms, rose diagrams, and empirical semivariogram contour plot to decide whether the assumption of isotropy is reasonable, or use the proposals presented by Cabaña (1987), Maity & Sherman (2012), or Weller & Hoeting (2020 to test this hypothesis. The stationarity can be verified using the proposal presented by Fuentes (2005), for example; while the assumption of normality can be verified using the proposals formulated by Yuan (2000) or Di Bernardino et al. (2017).
Another interesting research topic could be the extending of this procedure for health surveillance. The timely detection of various types of adverse health events is very important for the agencies responsible for health surveillance. For example, detect disease outbreaks such as asthma or influenza. The delay in the detection of and response to an adverse health event can result in the loss of lives and high economic costs. The early detection of these events allows preparing contingency plans to contain them. Using the spatial information that is given by the location can potentially enable localized outbreaks of a disease to be detected, or variations in regional patterns to be identified (Unkel, Farrington, Garthwaite, Robertson & Andrews 2012). The proposed spatial MCUSUM control chart could help detect these outbreaks or changes. This extension must consider, among others, the following observations: 1) risk adjust the outcome data before constructing a control chart, 2) find measures appropriate to evaluate the performance of the control proposed chart in the health care context, 3) adjust the population at risk, which it varies over time, 4) capture the space-time dependence presents in the data. Some interesting references on this topic are Sonesson & Bock (2003), Woodall (2006), and Tsui, Chiu, Gierlich, Goldsman, Liu & Maschek (2008).
It would be interesting to investigate a version of the proposed control chart for phase I, as well as the effect of the estimation in phase I on the performance of the chart in phase II.