Variable parameters memory-type control charts for simultaneous monitoring of the mean and variability of multivariate multiple linear regression profiles

Variable parameters (VP) schemes are the most effective adaptive schemes in increasing control charts' sensitivity to detect small to moderate shift sizes. In this paper, we develop four VP adaptive memory-type control charts to monitor multivariate multiple linear regression profiles. All the proposed control charts are single-chart (single-statistic) control charts, two use a Max operator and two use an SS (squared sum) operator to create the final statistic. Moreover, two of the charts monitor the regression parameters, and the other two monitor the residuals. After developing the VP control charts, we developed a computer algorithm with which the charts' time-to-signal and run-length-based performances can be measured. Then, we perform extensive numerical analysis and simulation studies to evaluate the charts’ performance and the result shows significant improvements by using the VP schemes. Finally, we use real data from the national quality register for stroke care in Sweden, Riksstroke, to illustrate how the proposed control charts can be implemented in practice.


Multivariate multiple linear profiles
For the kth sample of size n, with p response variables (profiles), Y k is a n × p matrix.Y k is a linear function of some independent variables x, so that: where X is a n × q + 1 matrix of explanatory (independent) variables,q is the number of independent variables, B is a q + 1 × p matrix of regression parameters, and E k is a n × p matrix of correlated error terms ( ε) , which follows a multivariate normal distribution (0, ) , where =    σ 11 σ 12 • • • σ 1p . . . . . . . . .
, and σ gh denotes the covariance between the error vector terms of gth and the hth response variables at each observation.Therefore, we can write Eq. ( 1) as:

Max-type and SS-type memory-type control charts
In this section, we describe two memory-type Max-type control charts, namely Max-MEWMA, and Max-MCU-SUM charts, which form a single statistic by taking the maximum value between the absolute values of two statistics (one for the process mean vector and the other one for the process variability).We also describe two SS-type memory type control charts, namely SS-EWMAe and SS-CUSUMe charts, which form a single statistic by adding the squared values of two statistics (again, one for the process mean vector and the other one for the process variability).As previously mentioned in the Introduction section, these Max-type and SS-type control charts are proposed by Ghashghaei & Amiri 20 and Ghashghaei & Amiri 21 , respectively.

Memory-type control charts using a Max operator
In this section, we describe two memory type Max-type control charts for regression parameters, namely Max-MEWMA and Max-MCUSUM charts, which were introduced by Ghashghaei & Amiri 20 .

Max-MEWMA control chart
In this section, we detail a single control chart for simultaneously monitoring of the process mean vector and variability ( B and matrices), by assuming that their values are known.First, we need to develop a statistic to represent the process mean vector.To monitor the process mean, the Hotelling's T 2 k statistic can be used to monitor the changes in the B matrix.The sample estimate of the B matrix ( B k ) is computed as: By changing B k into a p(q+1)×1 vector, we have: Next, we need to compute its average and variance-covariance matrix.For an in-control process, and since we assume that the parameters' values are known, the expected value of B k is equal to B .Therefore, we have: For its variance-covariance matrix, we have: gh is a q + 1 × q + 1 matrix equal to X T X −1 σ gh .
Next, we define: (1) where z 0 is the starting point and is equal to zero, is the smoothing parameter and its value can vary between 0 and 1.However, most commonly, = 0.2 is used.
Then, the following statistic is defined for monitoring the process mean vector: where H (q+1)p (.) is the chi-square cumulative distribution function with q + 1 p degrees of freedom, and �(.) is the standard normal cumulative distribution function.
To construct the statistic for monitoring the process variability, first, we define W k as: So that W k has a chi-square distribution with np degrees of freedom.Then, we define: where g 0 is the starting point and is equal to zero.The statistic for monitoring the variability is defined as: Finally, the ME k statistic is formed by combining C k and S k : Since this statistic only generates positive values, we only need an upper control limit (UCL) for this control chart, and its value is obtained using simulation to achieve any desired ARL performance.

Max-MCUSUM control chart
The statistic for monitoring the process mean vector for the Max-MCUSUM control chart is defined as: where The statistic for monitoring the process variability is: where v = log(τ ) τ τ −1 , τ is the multiplier with which the variance-covariance matrix shifts ( → τ ) , and similar to Ghashghaei and Amiri 20 we assumed τ =1.1 and L 0 = 0.
Finally, MC k is formed by combining U k and L k , as: Again, since this statistic only generates positive values, we only need a UCL for this control chart, and its value is obtained using simulation to achieve any desired ARL performance.

Memory-type control charts using a SS operator
In this section, we describe the SS-type control charts for the residuals introduced by Ghashghaei & Amiri 21 .They developed these kinds of control charts to monitor both regression parameters and residuals (four control charts).However, we only use the ones for monitoring the residuals, mostly because we already have our Maxtype ones for monitoring the regression parameters, and also, they concluded that in most situations their residual monitoring charts perform better than the other ones.

SS-EWMAe control chart
To compute the statistic for monitoring the mean vector ( P k ) , we first define: where H p (.) is the chi-square cumulative distribution function with p degrees of freedom, �(.) is the standard normal cumulative distribution function, z k = e k + (1 − )z k−1 , z 0 = 0, e k = e 1k , e 2k , . . ., e pk is the average residual vector in the sample k, and is the smoothing parameter.
Then, we have: where f k = n i=1 (e ik ) −1 (e ik ) T , and is the variance-covariance matrix of the error terms.Then we have: where V 0 = 0.
Finally, the SS-type statistic in the EWMAe scheme is defined as: The same as in the case of Max-type control charts, since this statistic only generates positive values, we only need a UCL for this control chart, and its value is obtained using simulation to achieve any desired ARL performance.

SS-CUSUMe control chart
The statistic for monitoring the mean vector in this scheme is: is the reference value, and T k was defined in the previous chart.
Similarly, the statistic for monitoring the process variability is: where is the reference value, and F k was defined in the previous chart.
Finally, the SS-type statistic in the CUSUMe scheme is defined as: Note that following Ghashghaei & Amiri 21 , in this paper we choose k 1 = 1 and k 2 = 1.5.

Design parameters in a variable parameters scheme
The adaptive scheme in which all the control chart (design) parameters are allowed to vary from sample to sample, is called a VP (Variable Parameters) scheme.In this paper, we consider two types of sample sizes with n 1 <n 2 , two types of sampling intervals with t 2 < t 1 , and two types of Type-I error probabilities with α 1 <α 2 .In addition to these parameters, we have to define two upper control limits UCL 1 and UCL 2 with UCL 2 < UCL 1 as well as two upper warning limits UWL 1 and UWL 2 , satisfying UWL 1 < UCL 1 and UWL 2 < UCL 2 .
In a VP scheme, we should have the following three constraints (each one is related to one design parameter) satisfied: By solving Eqs. ( 20)-( 22) together, P 0 , t 1 and α 2 are obtained as: Note that P 0 is the conditional probability of being in the safe zone while the process is in-control.After deter- mining the values of the UCLs and UWLs, which we will later show how to determine by introducing algorithms 1 and 2, we use the following sampling strategy in a VP scheme: • If at sample k, the statistic's output ∈ 0, UWL (k) , then the process is declared as being in-control and the parameters for the next sample must be n 1 , t 1 , UCL 1 , UWL 1 .• If at sample k, the statistic's output ∈ UWL (k) , UCL (k) , then the process is also declared as being in-control, but the parameters for the next sample must be n 2 , t 2 , UCL 2 , UWL 2 .• If at sample k, the statistic's output ∈ UCL (k) , ∞ , then the process is declared as being out-of-control and the corrective actions might be required, where UCL (k) ∈ {UCL 1 , UCL 2 } and UWL (k) ∈ {UWL 1 , UWL 2 } are the upper control and warning limits used for sample k=1,2,..., respectively.
For determining the UCL values, we assume that we have an FP scheme and set the values of each UCL separately (UCL for an FP control chart, and UCL 1 & UCL 2 for a VP control chart).We use the following algorithm for determining the values of UCLs.First note that in all the following algorithms, it is assumed that the average sampling interval is equal to one time unit, which results in ARL=ATS.Otherwise, it would have been ATS=t × ARL.Also, all the algorithms in this section should be run in an in-control state.In fact, obtaining the values of all the control chart parameters should be performed while the process is in-control.

Algorithm 1: Adjusting the UCL values
Step 1-Choose a value for α (the probability of Type-I error) and the sample size n.
Note that for the FP schemes there is only one α and one n.However, for the VP schemes we have two α s and two ns, therefore we set UCL 1 by using α 1 and n 1 , and UCL 2 by using α 2 and n 2 .
Step 3-Obtain the initial value for the UCL by generating and increasingly sorting 10000 in-control samples by using the statistic and choosing the [10000(1-α)] th value in the range.
Step 4-Run 10000 simulations and adjust the UCL so that you get ARL= 1 α .After determining the value of the UCLs, we should obtain the values of UWLs for the VP scheme (remember that the FP scheme has no UWL).For obtaining the UWL values, we can assume that we only have a variable sampling interval (VSI) scheme and compute each UWL separately with its corresponding parameters using the following algorithm (as is supposed to be done in a VSI scheme).

Algorithm 2: Adjusting the UWL values
Step 1-Choose the values of α , t 2 , E(t), the corresponding UCL value obtained via the previous algorithm, and the corresponding sample size ( Step 2-Compute P 0 using Eq. ( 23) and t 1 using Eq. ( 24).
Step 3-Choose the same statistic used for determining the corresponding UCL value.
Step 4-Run 10000 simulations and adjust each UWL with its corresponding UCL so that you get P 0 equal to the value you obtained in Step 3 and at the same time get ARL= 1 α and ATS= E(t)ARL=E(t) 1 α (note that if E(t)=1, as in our case, then ARL= ATS= 1 α ).

Performance measures
The Run length and time to signal based measures are the two most important control charts' performance measures.In an FP scheme, computing the average run length based measures are enough, considering one can multiply the average run length by the sampling interval to obtain the average time to signal.However, in a VP scheme, the average time to signal should be computed separately.Both average and standard deviation of run length (ARL and SDRL) as well as time to signal (ATS and SDTS), are important to consider.Although the SDRL and SDTS are always expected to be low, the ARL and ATS are expected to be as high as possible when the process is in-control and as low as possible when the process is out-of-control.
To compute the performance measures for an FP scheme, algorithm 1 can still be used with the only difference that here we have obtained the UCL value, and we are now only interested in computing the values of the ARL and SDTS in an out-of-control situation.
To compute the performance measures for a VP control chart, the following computer algorithm is developed and can be used.Note that, to reduce the paper size, we only present the algorithm for the case of the Max-MEWMA control chart, but it can easily be modified for the other proposed control charts as well.
Algorithm 3 Computing the performance measures in a memory-type VP scheme.

Simulation studies
In this section, we perform numerical analyses and simulation studies to evaluate and compare our developed adaptive control charts to one another in adaptive and one-adaptive conditions.We evaluate the performance of the proposed control charts under different shift scenarios and in different dimensions.
Although we report the values of the ATS, SDTS, ARL, and SDRL in the following tables, the comparisons are mainly made using the ATS values.All the simulation environments and the chosen values for the process and Table 1.ATS=ARL, SDTS=SDRL for FP and ATS, SDTS (ARL, SDRL) for VP schemes for shifts in the intercept vector ( β 0 ) and the error variation ( τ ), when p=2 and q=2.chart parameters as well as the shift sizes are the same as in Sabahno & Amiri 27 .The in-control ARL and ATS for all the considered control charts are set to 200 runs and 200 hrs (α = 0.005), respectively.The analysis is conducted for the case of two and six response variables (p=2 and 6), i.e. two and six multiple linear profiles.The following multiple regression models are used for the case of p=2: The error's variance-covariance matrix for this case is assumed to have the following elements: , where σ 1 and σ 2 are the standard deviations of the first and second profiles, respectively, and ρ is their correlation.For its in-control value, we have: 0 = 1 0.5 0.5 1 .
The sample size for the non-adaptive (FP) scheme is 4 and the sampling interval is 1hr.For the adaptive VP control chart, in which we need two values for each chart parameter, the following parameters' values are chosen: www.nature.com/scientificreports/ 004,E(t) = 1 hr,and t 2 = 0.1 hr.t 1 is computed by using Eq. ( 24) as 1.9 hrs, and the upper control and warning limits are computed using the algorithms outlined in Section "Design parameters in a variable parameters scheme".www.nature.com/scientificreports/For each sample size, a set of explanatory variables is required.The X matrix for n = 4 is assumed to be  and for the sample size 8 (only needed for the VP chart) it is considered ATS=ARL, SDTS=SDRL for FP and ATS, SDTS (ARL, SDRL) for VP schemes for shifts in the intercept vector ( β 0 ) and the error variation ( τ ), when p=6 and q=6.www.nature.com/scientificreports/ We also assume that each element in the error's variance-covariance matrix shifts with the same multiplier (τ ).In addition, the correlations between the responses in this section are assumed to be equal (regarding the p=6 case) and are fixed at ρ = 0.5.
The results for the p=2 case are presented in Table 1 (which contains separate and simultaneous shifts in the variability and the intercepts), Table 2 (which contains separate and simultaneous shifts in the variability and the first slopes), and Table 3 (which contains separate and simultaneous shifts in the variability and the second slopes).
The results in all these three tables show that in all the FP and VP control charts, as the intercept/slope or variability shift increases, the charts signal faster.Moreover, if the number of profiles whose intercepts/slopes shift increases from one to two, i.e. from (0.2, 0) to (0.2, 0.2) in Table 1, only the Max-MCUSUM chart shows a significant increase in the performance (decrease in the ATS value), and the performance of the other charts remains more or less the same.In addition, by comparing the charts' FP and VP schemes, we realize that all the charts show significant performance improvements if the VP scheme is used (with more than a 70% performance increase in some cases), and by comparing different control charts, it is clear that the Max-type control charts www.nature.com/scientificreports/mostly perform better than the SS-type control charts (only if one profile shifts, the SS-type control charts perform better and that also only in some cases of no or low variability shifts).As for the Max-type control charts, the Max-MEWMA chart mostly performs better as the variability shift increases and the Max-MCUSUM chart mostly performs better as the mean shift increases.
For the case of p=6, the real case adopted in Sabahno & Amiri 27 with the following model is used:  Since there are two sampling strategies in our adaptive scheme, we use n 1 = 8 and n 2 = 16 , and E(n) = 12 .Therefore, since we have two sets of sample sizes, we need two value sets for the explanatory variables as well.Again, we use the same value sets used by Sabahno & Amiri 27 , which we don't include in this paper to save space.
The results of this case are presented in Tables 4, 5 and 6, for separate and simultaneous shifts in the variability and the intercepts, separate and simultaneous shifts in the variability and the first slopes, and separate and simultaneous shifts in the variability and the second slopes, respectively.
The results in the p=6 problem show that while in the cases of slope shifts (Tables 5 and 6) the conclusions are almost the same as in the previous case (p = 2), the same does not completely apply to the case of the intercept shift (Table 4).The Max-type control charts' performance mostly gets worse (or their performance remains rather unchanged) as the shift in the intercept increases (with the Max-MCUSUM chart being the worst between the two).On the contrary, The SS-type charts' performance mostly gets better (or their performance remains rather unchanged) under a similar situation.In addition, except for some cases of no/low variation shifts, the y 4 = 0.04 + 0x 1 + 0x 2 + 0x 3 + 10.53x 4 − 0.47x 5 + 0.21x 6 + ε 4 , www.nature.com/scientificreports/Max-MEWMA control charts perform better than the other control charts.Moreover, the VP adaptive control charts are still mostly much faster than the FP charts.By comparing the p = 6 case with the p = 2 case, we realize that in the case of no variation shift ( τ = 1), all the charts perform worse when the mean shift is in the intercept.However, if the mean shift is in the slopes, and also in the case of low variation shift ( τ = 1.1) when there is no mean (intercept/slops) shift, all the charts perform better, but as the mean shift increases, the charts mostly tend to perform worse when the process dimension increases.Furthermore, in the cases of moderate/large variation shifts ( τ = 1.3 and 2), all the charts perform better in the case of p=6 compared to the case of p = 2.

A real case
To illustrate how one can implement the proposed control charts in practice, we study a real healthcare-related case.Stroke is one of the most common causes of death and disability in the world.In addition to severe consequences for individuals, stroke causes a high financial burden on societies.Intravenous thrombolysis within 4.5 hours from the stroke onset is an established treatment for ischemic stroke.The benefit of treatment reduces for every minute's delay (Darehed et al. 40 ), and thrombolysis delay times are key quality indicators of stroke care, and essential to monitor and maintain a good quality of stroke care.We study a dataset containing all the stroke patients who received thrombolysis in Sweden from 2016 until the year 2020 and were registered in the national quality register for stroke care in Sweden (Riksstroke; Asplund et al. 41 ).The study has been performed in accordance with the relevant guidelines/regulations.The Declaration of Helsinki was followed.Patients were informed about registration in Riksstroke, that their registered data could be used for research purposes, and their right to remove themselves from the registry at any time (opt-out consent).According to the Swedish Patient Data Act, data from national quality registers may be processed for research purposes without additional individual consent, if processing has been approved by an Ethics Review Board in accordance with the Ethical Review Act.The use of data from Riksstroke for this study was approved by the Swedish Ethical Review Authority (reference no.2021-06152-01).
The objective of this study is to monitor the efficiency of the stroke care process in terms of thrombolysis treatment delays.Therefore, we investigate whether the relationship between two correlated responses, y 1 = the time from stroke onset until getting the treatment (onset-to-needle time, ONT) and y 2 = the time from admission to the hospital until getting the treatment (door-to-needle time, DNT) in relation to three crucial covariates (patient characteristics) of x 1= Age, x 2= Sex, and x 3= Stroke Severity (as measured by NIH stroke score, i.e.NIHSS) is being kept constant over time or not.x 1 and x 3 are modeled as continuous variables and x 2 as a binary variable (0 for male and 1 for female).
We used the data from two recent years 2016 and 2017 (which were considered to be stable years in the stroke system), from all the hospitals in Sweden to estimate this relationship (profile).
After cleaning the dataset by removing the missing data (the proportion of missing data was relatively low compared to the overall dataset) and the erroneous entries (times more than 4.5 hours, which showed they had been entered wrongly), our first analysis was to see if ONT and DNT were normally distributed or not.The analysis revealed that their distributions were skewed, indicating that they deviated from a normal distribution and appeared more like Lognormal distributions.This outcome is commonly expected for time-related variables with many factors influencing them.Consequently, we applied a logarithmic transformation (base 10) to the data in order to approximate a normal distribution.The variance-covariance matrix of these response variables was estimated as: = 0.0399 0.0207 0.0207 0.0743 .
Then, we estimated the multivariate multiple regression model using R.The results were as follows.
We used the VP adaptive control charts (Max-MEWMA, Max-MCUSUM, SS-EWMAe, and SS-CUSUMe) to monitor these regression models over time.The design parameters for our control charts are the same as in our simulation study section (except for the sampling intervals) are set to n 1 = 4 , n 2 = 8 , E(n) = 6 , E(α) = 0.005 , α 1 = 0.004 , E(t) = 2 months and, t 2 = 1 month.The other chart parameters are computed using Equations ( 23)- ( 25) and the proposed algorithms in Sect."Design parameters in a variable parameters scheme" and can be seen in Tables 7, 8, 9 and 10 for each control chart.
After designing the control charts, we first checked whether the process was really under control in Phase-I or not (during the years 2016 and 2017).Otherwise, the estimated regression models and consequently the developed control charts would not be valid to be used in Phase -II.Researchers usually use non-adaptive (FP) www.nature.com/scientificreports/schemes to do so.The result showed that the process was in-control according to all the control charts (the details of this analysis are not included in this paper but can be requested from the corresponding author).For Phase-II, we employed the developed control charts for the years 2018 and 2019 to see if any assignable causes could be detected in those years.The control charts for the Max-MEWMA, Max-MCUSUM, SS-EWMAe, and SS-CUSUMe schemes can be seen in Figs. 1, 2, 3 and 4, respectively.More details regarding each sample can be seen in Tables 7, 8, 9 and 10.
Tables 7, 8, 9 and 10, from left to right, show the sample number (k), the sample size (the number of patients investigated in the sample), the cumulative number of samples up to the current sample, the sampling interval (in months) used to reach the current sample, the cumulative number of the sampling intervals up to the current sample (the time from the start of the process monitoring up to the current sample), the mean of the current sample's first and second response variables, the values of the sample statistics for the mean and variability, the value of the final statistic, the UWL and UCL values used for the current sample, and the status of the process based on the current sample, respectively.Remember that each control chart has two statistics, one for monitoring the mean vector and the other one for monitoring the variability.However, by using a Max or SS operator, a single final statistic will be formed and plotted in each control chart, i.e. the statistic in column ten.
As can be seen in Figs. 1, 2, 3 and 4 as well as Tables 7, 8, 9 and 10, the Max-MEWMA control chart was able to detect the out-of-control situations at samples eleventh (November 2019) and twelfth (December 2019).This control chart was able to signal first after 22 months and after observing 60 patients.The Max-MCUSUM control chart was able to detect the out-of-control situations in samples ninth (September 2019) and 11th (December 2019).This control chart was able to signal first after 20 months and after observing 44 patients.The SS-EWMAe control chart signaled at samples 14th (October 2019), 15th (November 2019), and 16th (December 2019).This chart was able to first signal after 21 months and after investigating 90 patients.Finally, the SS-CUSUMe control chart signaled at samples 13th (November 2019) and 14th (December 2019).This chart was able to first signal after 22 months and after investigating 80 patients.
One thing that is clear based on the obtained result is that in all the control charts, the statistic responsible for the mean vector monitoring (column eight) has caused the signal, meaning the shift has happened in the coefficient values and not the responses' variability.
Further investigation is required to discover the reasons for these signals and to see which ones are real assignable causes and which ones are outliers that can be ignored.Then, these assignable causes should be removed if they are undesirable.We might even need to update the profile's parameters if we discover that those assignable causes are desirable and should be kept.

Conclusions
In this paper, we improved the performance of four memory-type control charts for monitoring multivariate multiple linear profiles.These control charts are Max-MEWMA, and Max-MCUSUM control charts that use a single Max-type statistic and monitor the regression parameters, and SS-EWMAe and Max-CUSUMe control charts that use an SS-type statistic and monitor the residuals.We designed a VP adaptive scheme for all these control charts, in which all the design parameters can be varied throughout the process monitoring to increase their capability in detecting shifts.After that, we developed an algorithm with which the time to signal and run length-based performance measures of these charts could be measured.Then, we performed extensive simulations to evaluate these charts' performance under different shift sizes and types as well as different process dimensions.Two different cases of two profiles (p = 2)-two covariates (q = 2), and also, six profiles (p = 6)-six covariates (q = 6) were studied in this paper.
The results in the p = 2 and q = 2 case showed that: (i) as the intercept/slope or variability shift increases, all the FP and VP charts signal faster, (ii) if the number of profiles whose intercepts/slopes shift increases, only the Max-MCUSUM chart shows a significant performance improvement, (iii) all the charts show significant performance improvements if the VP scheme is utilized, (iv) the Max-type control charts mostly perform better than the SS-type control, and (v) the Max-MEWMA chart mostly performs better as the variability shift increases and the Max-MCUSUM chart mostly performs better as the mean shift increases.The results in the p = 6 and q = 6 case show that, in the case of slope shifts, the conclusions are more or less the same as in the case of p = 2 and q = 2.However, in the case of intercept shift, the main difference is that the Max-type control charts' performance mostly gets worse (or their performance remains rather unchanged) as the shift in the intercept increases (with the Max-MCUSUM chart being the worst between the two) and the SS-type charts' performance mostly gets better (or their performance remain rather unchanged).
Finally, we used a real dataset to estimate two profiles in a stroke care process and then developed and utilized the VP control charts to monitor those profiles over time to show how these charts can be implemented in real practice.
For future studies, implementing similar adaptive strategies for more advanced profiles such as non-parametric and nonlinear profiles can be suggested.Furthermore, developing and implementing other control charts to improve the healthcare-related processes in general, and the stroke care process in particular, could be a great contribution considering the availability of very few studies in this regard.

Figure 4 .
Figure 4. Max-CUSUMe VP control chart in the illustrative example.

Table 7 .
Details of process monitoring from Jan 2018 to the end of 2020 in the real case for the Max-MEWMA control chart.Significant values are in bold.
k n

Table 8 .
Details of process monitoring from Jan 2018 to the end of 2020 in the real case for the Max-MCUSUM control chart.Significant values are in bold.
k n

Table 9 .
Details of process monitoring from Jan 2018 to the end of 2020 in the real case for the SS-EWMAe control chart.Significant values are in bold.
k n

Table 10 .
Details of process monitoring from Jan 2018 to the end of 2020 in the real case for the SS-CUSUMe control chart.Significant values are in bold.