1 Introduction
Interest in monitoring and continuous improvement of medical outcomes has grown steadily over the last three decades. Results of monitoring exercises are used by regulatory organizations, hospital administrators, medical practitioners and other stakeholders to detect changes in performance and to compare practitioners or different hospitals. Current methods of monitoring medical outcomes focus on final outcomes of a given procedure, such as surgical failure rate. The methods used are based on those developed for monitoring in industrial contexts (montg2009). To take into account heterogeneity of expected outcomes among patients, riskadjusted monitoring was proposed (lovegrove1997; poloniecki1998; steiner2000). In riskadjusted monitoring, the expected outcome is not uniform among patients, but rather depends on key risk factors of the patient. Some recent developments include simultaneous monitoring of multiple outcomes (waterhouse2012) and monitoring that allows for underlying changes in performance over time (steiner2014).
In addition to patientspecific risks, there is heterogeneity among practitioners in the decisions they make and the processes they follow during the procedure. This issue has not received much attention in the literature. Furthermore, the final outcome of a procedure may be related to outcomes of preceding stages within the procedure. To address both issues, sibanda2016 proposed monitoring outcomes of intermediate stages in addition to the final outcome, adjusting for heterogeneity in patient risk and processes used by practitioners. A model for the procedure in its entirety was developed incorporating patient risks, processes and outcomes at all stages. A test statistic based on the ratio of observed () to expected () number of failures in consecutive patient subgroups was used for each outcome. While the test statistic has a clear clinical interpretation, it is slow in detecting small shifts and is unsuitable for rare outcomes. In this article we propose monitoring the regression coefficients of the model that represents the relationships among the outcomes, processes and patient risk factors relevant to the procedure. To ensure quick detection of small gradual shifts we use a sequential monitoring approach where the chart statistic is updated at each successive observation.
In the statistics and econometrics literature, considerable attention has been paid to problems of monitoring constancy of parameters in statistical models. nyblom1989 derived a Lagrange multiplier (LM) test based on likelihood score equations for a timeseries model and showed this to be the locally most powerful test for the alternative that the parameters follow a martingale process. hansen1992
extended the LM test to linear regression models, and
hjort2002 suggested a general class of likelihood scorebased structural change tests. In this article, we develop a monitoring procedure for model coefficients based on likelihood score equations.In the statistical process control literature, profile monitoring is used for model based processes where the quality of the process or a stage of the process is represented by a functional relationship between a response variable and one or more explanatory variables.
noorossana2011 describe profile monitoring in a number of settings. The main objective of profile monitoring is to check the stability of a functional relationship over time. mandel1969 and hawkins1993give examples of model based profile monitoring procedures. In these references, the response variable was assumed to follow a Normal distribution and regression coefficients were estimated using ordinary least squares (OLS).
skinner2003 and jeark2003incorporated use of the generalized linear model in profile monitoring of Poisson and Gamma distributed response variables.
jeark2005 and Shang2013 developed a procedure based on monitoring deviance residuals using Shewhart charts for a multistage process with nonGaussian response variables. Shang2013 pointed out that an exponentially weighted moving average (EWMA) scheme would be required to detect small shifts in the multistage process. yeh2009developed a profile monitoring procedure using logistic regression for a single binary response variable. In this paper, we develop a multivariate EWMA (MEWMA) procedure based on likelihood score equations for simultaneous monitoring of the vector of model coefficients of a multistage process. In this way we monitor stability of the functional relationship between multiple response variables and associated predictors.
There are a number of issues to consider when monitoring a multivariate vector of parameters in a model. For example, monitoring model parameters is more complex than the traditional SPC approach of monitoring the response variables themselves. Also, a medical procedure can have categorical outcome variables. Our approach is flexible in that it can handle such processes, since the likelihood function and score equations can be obtained in a straightforward manner once the assumed distributions of the response variables are known. In the next section, we describe a motivating example and represent the process using a graphical model presentation. Then in Section 3, we present the general formulation of our proposed monitoring procedure, together with derivation of the associated chart statistic. Performance of the proposed chart is assessed using simulation studies in Section 4 and the results are discussed in Section 5.
2 Motivating example
The motivating example is the delivery process in a maternity unit. The example is described in detail by sibanda2016 and is revisited briefly here. The medical procedure of interest is infant delivery in the labour ward of a maternity unit. Four maternal outcomes, represented by , are identified. These outcomes are temporally ordered as they occur at various ordered stages of the process. Temporal ordering allows us to represent correlations among the outcome variables using a graphical model framework. In addition to being influenced by outcomes at upstream stages, each depends on a number of process factors, , and patient risks, that can also be included in the graphical model. The graphical representation of a procedure extends the literature on the regression adjustment approach and multistage control charts developed in industrial applications (hawkins1991; hawkins1993; zantek2002). The graphical representation of the delivery procedure is shown in Figure 1.
2.1 Model specification for the delivery procedure
There are three major stages in the delivery process with a total of four maternal outcomes defined as:
Patients who do not undergo labour due to an elective or emergency Cesarean section have been omitted from this study since the only outcome observed for them is . Each outcome variable is influenced by one or more process variables given by
and risk factors given by
A combination of literature (thorngren2001; cameron2006; rouse2009; hirayama2012) and empirical evidence was used to determine the structure of the graphical presentation as shown in Figure 1.
The graphical model in Figure 1
can be viewed as a Bayesian network
(pearl1986; lauritzen1988). Note that for clarity of presentation we have deviated from standard graphical model presentation by not including separate edges for model coefficients. A Bayesian network is a directed acyclic graph (DAG) with node setrepresenting a set of random variables,
having a joint probability distribution function that can be written as
(1) 
The term represents the set of parent nodes of the node . The power of a DAG representation is that once the structure is known, the joint probability distribution of can be written in the form of equation (1) using the conditional independence axioms introduced by dawid1979. In equation (1), each node is conditionally independent of all other nodes, given its parent nodes. Based on the graphical structure in Figure 1 we can write:
where is a vector of parameters characterizing the relationship between variables.
2.2 Process model and monitoring
For each patient, , we observe the outcome vector , the process vector, and the riskfactors . Using generalized linear model (GLM) notation we write
(2) 
where is an appropriate link function and , for and . In general, the process model may comprise a mix of continuous, discrete, ordinal or nominal variables. In such cases, it is possible to have a different link function for each outcome variable
. In our case all outcome variables are binary in nature and here we choose to use the logit link for all
. We therefore have:(3) 
for all and . Our proposed procedure is based on monitoring the model coefficient vector where for all .
3 Score based Multivariate EWMA chart
The scorebased multivariate EWMA (MEWMA) chart is the monitoring scheme we introduce to overcome the problems associated with rare outcomes and detection of small shift (sibanda2016). The MEWMA (lowry1992), accumulates information from past observations for quick detection of small shifts. Given a vector, , of quality indicators for patient , the MEWMA chart is based on the statistic
where the MEWMA vector is calculated using
and is the covariance matrix for given by
is a matrix with diagonal element for , and zero elements on the offdiagonal.
is the identity matrix and
is the covariance matrix for . The parameter is called the smoothing parameter, and it determines the relative weight of current and past observations for the quality characteristic. If there is no a priori reason to use different weights for the quality characteristics, the diagonal elements of are all set equal to a constant . When , montg2009 shows that the covariance matrix simplifies toThe MEWMA chart signals at the first observation at which , where is a control limit chosen to achieve a specified incontrol average run length (ARL).
In our case, the multivariate quality indicator of interest is the vector of regression coefficients for the functional relationships between the response variables and the process variables and risk factors. To detect a change in the functional relationship, we propose monitoring the score vector of the regression coefficients. For a graphical model with response variables , process variables , patientrisk variables and parameter vector of length , the likelihood function is given by
where is a vector of observed values of for patients and is the density or probability distribution function of . The loglikelihood function is given by
(4) 
Let be a length vector of first derivatives of ,
Then has an asymptotic variate normal distribution with mean zero and covariance matrix given by the Fisher Information matrix, , where is a matrix of second derivatives of with respect to (rao1973).
A common approach for testing the stability of a dimensional vector of regression coefficients is based on a cumulative sum of standardised scores at observation of the form
where is the specified value of
under the null hypothesis of no change in the regression coefficients. The components of
tend in distribution to independent Brownian motions under the null hypothesis (hjort2002) and are used for separate tests of the parameters. When is unknown, a maximum likelihood estimate can be determined from data from a previous time period (Phase 1 data) and the estimated cumulative standardized score becomesOur approach is based on the cumulative sum of standardised scores, but with different weights given to current and past observations. This gives a scorebased multivariate exponentially weighted moving average (MEWMA) chart. This scorebased MEWMA chart gives greater weight to more recent observations and will therefore be quicker at detecting smaller gradual changes in the process parameters.
When the MEWMA is applied to the score vector , we calculate
at time , for i.i.d with and if the process is in control. This gives and . We would then chart
with control limit identified through simulation to achieve a desired incontrol .
3.1 Score equations and information matrix for the delivery process
For each , let be a dimensional vector and , where is the number of regression coefficients for and is the corresponding design matrix. For the model in Figure 1 and equation 3, the loglikelihood function is given by
(5) 
with dimensional score vector with elements
(6) 
for . The observed Fisher information matrix is a block diagonal matrix given by
where
The block diagonal structure of the information matrix suggests a natural partition of the regression coefficient vector into independent sets that can be used to monitor subsets of the coefficients.
4 Simulation study
In this section, we present a simulation procedure and results for how the scorebased MEWMA chart reacts to different types of shifts or changes at different stages of the multistage process for a large tertiary hospital. We check the efficiency of this score based chart in detecting small and large shifts in the process overall quality due to changes in process variables, risk factors or in upstream stage outcomes. The aim of the simulation is to obtain empirical estimates of average run length (ARL) as a measure of chart performance. The upper control limit used in this case is the one obtained when we set the chart to achieve in control ARL of 200 under the assumption that Phase I parameters are known without error. We used data from the Southmead maternity unit in Bristol, updated from that published in a previous paper (sibanda2007) to obtain model parameter estimates and these are shown in Table 1. These model parameter estimates are used as the basis for the simulation procedure for testing the proposed control chart.
Parameter  Value 

1.724  
0.730  
1.682  
1.262  
0.597  
0.342  
0.467  
0.758  
0.316  
1.140  
0.482  
1.267  
0.374 
There is an average of about 6400 deliveries a year at Southmead Hospital, which gives an average of 533 deliveries a month. Hence, for each hypothesised shift, a single data set of size 500 patients was generated for model development to correspond to the size of dataset available for monitoring at monthly intervals. This size of 500 fulfills the asymptotic assumption of the scorebased MEWMA chart. We consider four types of shifts.

An additive shift in a given regression coefficient(s): , where is the linear predictor excluding the term.

Simultaneous additive shifts in pairs of regression coefficients.

An additive shift in the mean response: . Values of were restricted to where .

A more appropriate approach for binary outcomes is to consider shifts in the odds ratio
(steiner2000). Under current conditions, we assume an odds ratio, of 1. If a shift occurs in one of the outcome variables the odds ratio is so thatThis means
For each shift type, a range of values of are used and for each value of , the chart is constructed 5,000 times and the run length recorded for each run. The arithmetic mean of the run lengths is used to estimate the outofcontrol ARL at which the scorebased MEWMA chart is likely to detect the proposed shift.
First, we study performance of the chart when the shift in the functional relationship occurs in , the effect of process variable on when remains unchanged. These kinds of changes are likely to occur with changes in staff or equipment. For illustration, changes of various sizes are introduced to and , the effects of (use of mechanical instruments) on (occurrence of 3rd or 4th degree tear) and on (occurrence of postpartum haemorrhage), respectively. The empirical outofcontrol ARLs under various shifts are shown in Table 2, indicating that detecting a change in the functional relationship of the process depends on the size of change and the size of the parameter affected. In this case, , and when an equivalent change is applied to both parameters, with all other parameters kept constant, detection of changes occurs more quickly for a shift in than it does for a shift in .
Change factor  Shift in  Shift in 

0.2  197.9  195.7 
0.4  194.4  183.5 
0.6  190.2  169.5 
0.8  187.9  153.3 
1.0  181.6  137.6 
1.2  174.8  123.8 
1.4  166.8  112.8 
1.6  157.9  102.4 
1.8  148.4  95.1 
2.0  140.2  88.9 
2.2  133.4  84.7 
2.4  123.6  80.2 
2.6  117.4  77.1 
2.8  111.0  73.7 
3.0  102.8  72.2 
3.2  98.2  70.3 
3.4  91.4  69.1 
3.6  86.9  68.3 
3.8  82.3  67.2 
4.0  77.3  66.3 
In our next investigation we consider performance of the chart when there is a shift in the effect of an upstream outcome variable in a downstream outcome. Such changes are likely to occur when changes in management of upstream adverse events occur that could be due to policy or staff changes. For illustration, changes of various sizes are introduced to and . The results are shown in Table 3 and show that detection of changes occurs more quickly for equivalent shifts in coefficients that are larger. In conjunction with the results in Table 2, we see that it is the size of the coefficient and size of shift that determines how quickly a change is detected. The type of coefficient, that is whether it is associated with a process variable, upstream outcome or risk factor does not seem to have an impact on how quickly shifts are detected. This is to be expected since the risk factors and process variables play the same role in the model.
Change factor  Shift in  Shift in  Shift in 

0.2  198.3  197.6  196.8 
0.4  195.9  195.4  172.8 
0.6  193.1  190.2  143.5 
0.8  193.8  179.2  129.4 
1.0  189.3  172.8  119.1 
1.2  187.5  155.8  88.4 
1.4  180.3  147.3  113.3 
1.6  173.6  137.0  62.8 
1.8  169.1  129.1  87.0 
2.0  163.7  114.5  22.0 
2.5  147.8  89.2  57.5 
3.0  129.7  78.1  97.5 
3.5  114.7  68.3  88.7 
4.0  100.0  53.4  60.0 
In our next investigation, we consider simultaneous shifts in more than one coefficient, indicating shifts in multiple parts on the process. For multistage health care procedures where there is no tight control of the process unlike in manufacturing, it is realistic to expect that changes can occur simultaneously at various stages. The chart we propose presents an approach for detection of any such changes at a global level. For illustration, we consider simultaneous shifts in the pairs and in . The results are shown in Figures 2 and 3 and indicate that shifts in pairs of coefficients are detected more quickly than shifts in a single coefficient (dotted line).
Next, we investigate chart performance when a shift occurs in the mean of an upstream outcome variable that affects a downstream outcome. Such a change represents an outcome shift whose source is unknown. A typical case where such a shift may occur is where small shifts occur in various parts of the functional relationship that collectively result in an overall mean shift. These kinds of changes are likely to occur in most health procedures that are complex in nature with multiple contributing components. For illustration, we investigate performance of the chart when there is a shift in , the probability of a 3rd or 4th degree tear occurring for patient . We consider both an additive shift and a shift in the odds ratio. The results are shown in Table 4. As expected the larger the shift, the lower the outofcontrol average run length.
Change factor  ARL for:  

Additive shift  Odds ratio shift  
0.2  188.2  188.2 
0.4  169.7  138.1 
0.6  153.8  137.3 
0.8  138.5  135.3 
1.0  124.9  130.6 
1.2  114.0  123.2 
1.4  103.8  116.0 
1.6  95.1  109.5 
1.8  88.3  102.6 
2.0  82.0  96.5 
2.2  75.5  91.1 
2.4  71.1  86.4 
2.6  67.7  81.8 
2.8  63.8  78.1 
3.0  60.0  74.6 
3.2  57.1  71.3 
3.4  54.3  68.6 
3.6  52.3  65.7 
3.8  49.8  63.2 
4.0  46.9  61.1 
Therefore the chart we propose has the ability to detect changes of various types in the process. The greater the shift, whether it is in a single coefficient, an outcome mean or in multiple coefficients, the quicker it is to detect. The chart can therefore be used as a tool for ongoing monitoring of a multistage procedure at the global level. An outofcontrol signal would then be followed by detailed investigations of various parts of the process.
5 Discussion
In this article, we proposed a new multistage multivariate control chart based on likelihood score equations to monitor the outcomes of health care procedures. A multistage approach is used to track outcomes at all stages through monitoring the coefficients of a model that represents the relationships among the outcomes, processes and patient risk factors relevant to the procedure. To ensure quick detection of small gradual shifts we used a sequential monitoring approach where the chart statistic is updated at each successive observation. We use likelihood score equations to construct a multivariate EWMA chart statistic. The advantage of our scorebased MEWMA proposed chart is that it can be designed to detect small or large shifts. Moreover, our approach is flexible in that it can handle processes with different types of variables since the likelihood function and score equations can be obtained in a straightforward manner once the assumed distributions of the response variables are known. Using simulations, we demonstrated the sensitivity of our proposed approach in reacting to various shift in the process and the effect of upstream outcomes on monitoring end stage outcomes. The chart showed efficient performance in detecting small shift sizes. Inclusion of risk factors in the model mean that factors outside of the practitioner’s control are accounted for. Signals that occur should be followed up with a review of the cases with adverse outcomes to determine the appropriate action to take.
Conflict of Interest: None declared.
Acknowledgements
We thank Professers Stefan Steiner and Jock Mackay for useful discussions at the initial phase of this work.
Comments
There are no comments yet.