Email updates

Keep up to date with the latest news and content from BMC Medical Research Methodology and BioMed Central.

Open Access Research article

Statistical process control of mortality series in the Australian and New Zealand Intensive Care Society (ANZICS) adult patient database: implications of the data generating process

John L Moran1*, Patricia J Solomon2 and and for the ANZICS Centre for Outcome and Resource Evaluation (CORE) of the Australian and New Zealand Intensive Care Society (ANZICS)

Author Affiliations

1 Department of Intensive Care Medicine, The Queen Elizabeth Hospital, 28 Woodville Road, Woodville, SA 5011, Australia

2 School of Mathematical Sciences, University of Adelaide, Adelaide, SA 5000, Australia

For all author emails, please log on.

BMC Medical Research Methodology 2013, 13:66  doi:10.1186/1471-2288-13-66


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2288/13/66


Received:7 January 2013
Accepted:17 May 2013
Published:24 May 2013

© 2013 Moran and Solomon; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Statistical process control (SPC), an industrial sphere initiative, has recently been applied in health care and public health surveillance. SPC methods assume independent observations and process autocorrelation has been associated with increase in false alarm frequency.

Methods

Monthly mean raw mortality (at hospital discharge) time series, 1995–2009, at the individual Intensive Care unit (ICU) level, were generated from the Australia and New Zealand Intensive Care Society adult patient database. Evidence for series (i) autocorrelation and seasonality was demonstrated using (partial)-autocorrelation ((P)ACF) function displays and classical series decomposition and (ii) “in-control” status was sought using risk-adjusted (RA) exponentially weighted moving average (EWMA) control limits (3 sigma). Risk adjustment was achieved using a random coefficient (intercept as ICU site and slope as APACHE III score) logistic regression model, generating an expected mortality series. Application of time-series to an exemplar complete ICU series (1995-(end)2009) was via Box-Jenkins methodology: autoregressive moving average (ARMA) and (G)ARCH ((Generalised) Autoregressive Conditional Heteroscedasticity) models, the latter addressing volatility of the series variance.

Results

The overall data set, 1995-2009, consisted of 491324 records from 137 ICU sites; average raw mortality was 14.07%; average(SD) raw and expected mortalities ranged from 0.012(0.113) and 0.013(0.045) to 0.296(0.457) and 0.278(0.247) respectively. For the raw mortality series: 71 sites had continuous data for assessment up to or beyond lag40 and 35% had autocorrelation through to lag40; and of 36 sites with continuous data for ≥ 72 months, all demonstrated marked seasonality. Similar numbers and percentages were seen with the expected series. Out-of-control signalling was evident for the raw mortality series with respect to RA-EWMA control limits; a seasonal ARMA model, with GARCH effects, displayed white-noise residuals which were in-control with respect to EWMA control limits and one-step prediction error limits (3SE). The expected series was modelled with a multiplicative seasonal autoregressive model.

Conclusions

The data generating process of monthly raw mortality series at the ICU level displayed autocorrelation, seasonality and volatility. False-positive signalling of the raw mortality series was evident with respect to RA-EWMA control limits. A time series approach using residual control charts resolved these issues.

Keywords:
Statistical process control; Time series; Autocorrelation; Seasonality; Volatility; Exponentially weighted moving average smoothing; Autoregressive moving average models; GARCH models

Background

Statistical process control (SPC), deriving from Shewart’s work in 1920-30 and in the 1950’s with Deming’s refinements [1], has been more recently applied in health care and public health surveillance [2], generating considerable interest in the general [3-5] and specialist medical literature [6-10]; and has been subject to a detailed exposition from a “quality-in-medicine” perspective [11]. Important statistical principles underlying SPC or control-chart methodology are those of the monitored process being “in control” and subject to the independence of observations [12]. The presence and impact (possible increase in frequency of false alarms) of process autocorrelation in industrial/engineering series have long been documented [13-16]. Somewhat surprisingly, little formal attention has been directed to this problem in the bio-medical literature [17,18], one review suggesting that there was “…limited advice on how to manage [autocorrelation]…” [5].

We have previously drawn attention to the data-generating mechanisms of overall monthly mortality series, at the aggregate level, from a bi-national intensive-care (ICU) database, where persistent autocorrelation (to lag24) was evident in a seasonal ARIMA (auto-regressive integrated moving average) model of the mortality series [19]. We now extend this study to further characterise the data generating process of mortality series at the individual ICU level and the impact of autocorrelation upon (i) mortality monitoring using EWMA (exponentially weighted moving average) control charts and (ii) time-series modelling of the data process using residual control charts.

Methods

As previously described [19,20], the ANZICS (Australian and New Zealand Intensive Care Society) adult patient database [21] was utilised to define an appropriate patient set, 1995-(end)2009. Physiological variables collected in accordance with the requirements of the APACHE III algorithm [22,23] were the worst in the first 24 hours after ICU (intensive care unit) admission, and all first ICU admissions to a particular hospital for the period 1995-2009 were selected. Records were used only when all three components of the Glasgow Coma Score (GCS) were provided; records for which all physiologic variables were missing were excluded, and for the remaining records, missing variables were replaced with the normal range and weighted accordingly. The mortality endpoint was at hospital discharge. Exclusions: unknown hospital outcome; patients with an ICU length of stay ≤ 4 hours, and patients aged < 16 years of age. Access to the data was granted by the ANZICS Database Management Committee in accordance with standing protocols; local hospital (The Queen Elizabeth Hospital) Ethics of Research Committee approval was waived.

Statistical analysis

(i) Monthly raw (risk-unadjusted) and risk-adjusted (RA) mortality time series at the individual ICU were generated. Risk adjustment was undertaken, generating the “expected” series, using a random coefficient logistic model (intercept as ICU site and “slope” as (centred) APACHE III score; unstructured covariance using adaptive quadrature, estimated via the Stata™ module “xtmelogit” [24]), as previously described in detail [20], and extended to both ventilated and non-ventilated patients. No formal adjustment for potential seasonality (trigonometric seasonality using sine/cosine functions or monthly dummy variables) was undertaken. Individual ICUs were allocated an identifier based upon a random number sequence.

(ii) Graphical inspection of the mortality series and formal testing of normality to confirm that the “…distributions of… (observed) and… (predicted) [series] … were sufficiently similar and are robustly normal and symmetrical…” [25]. Classical seasonal decomposition [26] was undertaken using the “decompose” module in R statistical software (Version 15.2 [27]). Autocorrelation plots (scatterplot grid of series versus lagged values) were performed via the R user-written module “lag1.plot” [28].

(iii) Generation of EWMA charts with confidence limits.

a. assuming iid (independent and identically distributed) observations, the EWMA statistic (zi) is defined as: λxi + (1 − λ)zi−1 and the variance (σ2) as <a onClick="popup('http://www.biomedcentral.com/1471-2288/13/66/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/13/66/mathml/M1">View MathML</a>, where 0 < λ ≤ 1 is a constant (smoothing parameter) [29].

a. For the variance of (non-stationary) auto-correlated series, we followed Montgomery & Mastrangelo [15]: division of the sum of squared (prediction) errors for optimal λ by n; leading to the plotting of a moving centre-line EWMA control chart [12].

a. Default values (“optimal”) in Stata™ statistical software for λ were chosen to minimize the in-sample sum-of-squares forecast errors [30], a method also recommended by Montgomery and Mastrangelo [15]; albeit small values of λ may inhibit the detection of large sudden process shifts; the “inertia” phenomenon [31].

a. Average run length (ARL): that is, the average number of “points”, when the data-generating process is in fact in-control, plotted before out-of-control is declared (ARL0). For instance, with iid observations and a Shewhart control-chart with three sigma limits, ARL0= 1/p=1/0.0027=370 (where p is the probability that any point exceeds the control limits [32,33], when the data-generating process is in fact in-control). Under the iid assumption, for various mortality series and values of λ, scenario based increments of the (mean of the) underlying series were computed using Statgraphics® Centurion XVI statistical software [34].

a. Using conventional SPC methods, EWMA control limits (at 3 sigma) were applied to the raw mortality series using the expected series as reference process; that is, RA control limits were generated.

(iv) Establishment of time-series models at the individual ICU level was based upon classic Box-Jenkins methodology (autoregressive moving average (ARMA) models) with investigation of (G)ARCH ((Generalised) Autoregressive Conditional Heteroscedasticity) effects [35,36], as previously described [19].

a. A stationary time series {xt; t = 0, ± 1, ± 2, …} has an autoregressive moving average (ARMA(p,q)) structure: xt = ϕ1xt − 1 + … ϕpxt − p + ωt + θ1ωt − 1 + … θqωt − q where ϕ1, ϕ2, …, ϕp are the “autoregressive” (AR) coefficients relating the value of x at time t to its past p values, and θ1, θ2, …, θq are the “moving average” (MA) coefficients, relating the current “white-noise”,ωt, to its past q values and <a onClick="popup('http://www.biomedcentral.com/1471-2288/13/66/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/13/66/mathml/M2">View MathML</a>. If xt has a non-zero mean (μ), then a constant α = μ(1 − ϕ1 − … − ϕp) is introduced into the structure. An integrated series accumulates (some) past effects and is therefore non-stationary. A series is integrated, say, of order 1 (I(1)) if the changes (or differences: Δxt = xt − xt − 1) of the series generate stationarity (I(0)), leading to the expanded ARIMA model (ARIMA(p,d,q)), where d is the degree of differencing [37]. This being said, careful attention was directed to the question of trend versus difference stationarity [38], especially in medical series where, as opposed to stochastic random walks, “deterministic” trends may be present. [39,40].

a. Model diagnostics: the use of auto- (ACF) and partial-autocorrelation (PACF) function displays, testing for the presence of a unit-root (ADF (augmented Dickey-Fuller) and DF-GLS (modified Dickey–Fuller t test) tests [30] and variants), residual white-noise (Bartlett’s periodogram-based- and Portmanteau (Q)-test) and seasonality were undertaken after Shumway & Stoffer [41] and as previously described [19].

a. Volatility of the (squared) residuals (ϵ) of the mean equation (conditional heteroscedasticity [42]) was checked using the PAC of the squared residuals and the user-written Stata™ “armadiag” module [43]; that is, ARCH and GARCH effects ((Generalised) Autoregressive Conditional Heteroscedasticity of the error variance process). For an ARCH model, the mean equation is yt = xtβ + ϵt and the variance equation <a onClick="popup('http://www.biomedcentral.com/1471-2288/13/66/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/13/66/mathml/M3">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2288/13/66/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/13/66/mathml/M4">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2288/13/66/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/13/66/mathml/M5">View MathML</a> are the squared residuals (innovations) and γi are the ARCH parameters; the conditional variance is thus modelled as an AR process. A GARCH(m,k) model includes lagged values of the conditional variance

<a onClick="popup('http://www.biomedcentral.com/1471-2288/13/66/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2288/13/66/mathml/M6">View MathML</a>, where δi are the GARCH parameters (an ARMA process) [19,44]. Exploration of different error term distributions (normal, t and generalised error) was also undertaken [30].

a. Under the conditions of an appropriately specified time-series model, the behaviour of the residuals was investigated, after Alwan and Roberts [45], on the basis that a shift in the mean of a time series is transmitted to the residuals [46].

b. As residuals are assumed to be independent (white-noise: a sequence of iid random variables with finite mean and variance, all ACFs being [close to] zero [47]), standard control chart methods were used to generate residual-EWMA charts [33]. Thus, determination of the residual-EWMA smoothing parameter (λ) was based upon methods for independent observations.

b. Control limits were also determined using standard errors (3×) of the one-step-ahead forecasts [45].

a. Model selection was guided by penalized information criteria (Akaike (AIC) and Bayesian (BIC) information criteria) [48].

a. Formal exegesis proceeded using a single exemplar complete ICU series (1995-(end)2009).

(v) Graphical displays: line-graphs of series were produced for appropriate illustration of relevant stages of analysis

a. Line graph(s) of the raw series were produced with 3*SE control limits of the expected series.

a. EWMA control limits (including residual control charts) were generated using default values of “optimal exponential coefficient” in Stata™ statistical software [49].

a. Values of λ for scenario based increments (say, 5% or 10%) of target mean were calculated using the SPC module of Statgraphics® statistical software [34] and appropriate 3*SE control limits of the expected series as in (a) above or EWMA line graphs were produced as in (b) above.

Results

The overall data set, 1995-2009, consisted of 491324 records from 137 ICU sites; mean (hospital) mortality was 14.07%. The random coefficient logistic regression model (Hosmer-Lemehsow statistic 62.97, ROC area under the curve 0.89) generated an overall predicted mortality probability of 0.1407 (SD 0.0202, range 0.00004-0.993). Over the 137 sites mean raw and expected (RA-) mortalities ranged from 0.012(0.113) and 0.013(0.045) to 0.296(0.457) and 0.278(0.247) respectively.

Of the raw mortality series from the 137 ICUs, 71 had continuous monthly data (excluding missing values or zero monthly mortality) for assessment up to or beyond lag40. For 25 of these series (35%), there was a significant Q test (null hypothesis being that the series is white noise) and autocorrelation through to lag40. Thirty six had continuous monthly data (excluding missing values) for ≥ 72 months; all series demonstrated marked seasonality and 30 demonstrated an obvious trend decline in mortality. Of the expected mortality series, 72 had appropriately assessable data to lag40 and in 46 (64%) there was a significant Q test and autocorrelation through to lag40. Similarly, in the same 36 series with continuous (raw) monthly data for ≥ 72 months, all expected mortality series demonstrated marked seasonality and 30 demonstrated an obvious trend decline in mortality.

Data from site “4” over 1995-2009 was used to generate an exemplar mortality time series. The mean raw mortality was 0.139(0.047) with skewness 0.216 and kurtosis 2.53; and the expected mortality was 0.138(0.028) with skewness 0.361 and kurtosis 3.47. The Shapiro-Wilk normality test was not rejected for either series (P =0.23 for both series). Kernel density estimates of raw and expected mortality are seen in Figure 1 (upper panel), with obvious difference in the degree of kurtosis between the two series. Time series plots, 1995-2009, for raw and expected (RA-)mortality are seen in the lower panel; a gradual time-decline in mortality for both series is evident. Additive seasonal decomposition of both series is seen in Figure 2, revealing marked seasonality and a trend decline in mortality. Autocorrelation plots are seen in Figure 3, showing correlation (positive and negative) decreasing variably with increase in lag in both series.

thumbnailFigure 1. Upper panel: Kernel density estimates of raw and expected mortality. Lower panel: Time series (1995-2009) of mean monthly raw and expected mortality.

thumbnailFigure 2. Decomposition plot of raw (upper panel) and expected (lower panel) mortality series.

thumbnailFigure 3. Lagplot of series versus lagged values (to lag24); upper panel, raw mortality; lower panel, expected mortality.

Figure 4 displays a plot of raw mortality series with control limits as 3SE of expected mortality (upper panel) and a scenario based mortality increment of 5% (5% false positive rate and desired ARL= 6 months) with control limits as 3SE of expected mortality. Frequent signalling is seen in both panel-plots. Figure 5 shows a plot of the raw mortality series with (fixed) EWMA 3 SE control limits derived from a projected 5% (upper panel) and 10% (lower panel) increment in expected mortality, assuming: an in-control ARL of 370, mean (expected) mortality 0.1381(0.0276) and target mean (expected) mortality of 0.145 (5% increment) and 0.152 (10% increment), for an EWMA λ = 0.02 and 0.05, respectively (calculations preformed in Stagraphics®). For both 5% and 10% projected increments of expected mortality, the raw mortality series signalled frequently, mainly in the early periods. Figure 6 shows the same scenarios with a time-varying variance EWMA control chart; again, there was frequent signalling of the raw mortality series.

thumbnailFigure 4. Upper panel: raw mortality series with 3SE control limits of expected mortality; lower panel: EWMA (λ=0.51) of raw mortality series for anticipated 5% increase in raw mortality, 5% fale positive rate and desired ARL=6 months with 3SE control limits of expected mortality.

thumbnailFigure 5. Raw mortality series with EWMA fixed 3SE control limits: upper control limit (red line), lower control limit (green line), signal (navy line); for anticipated 5% (upper panel) and 10% (lower panel) increments in expected mortality.

thumbnailFigure 6. Raw mortality series with EWMA time-varying 3SE control limits: upper control limit (red line), lower control limit (green line), signal (navy line); for anticipated 5% (upper panel) and 10% (lower panel) increments in expected mortality.

The autocorrelation evident in the raw and expected mortality series suggested a formal time series approach to SPC:

(i) Raw mortality: both the DF-GLS and ADF tests (with trend) rejected the null-hypothesis of presence of a unit-root and the series was de-trended using linear regression (raw mortality against time) and the residuals (also not evidencing a unit-root) of the linear regression model were used for subsequent formal analysis. The de-trended series from the raw mortality displayed seasonality but, not surprisingly, no trend decline (graphics not shown). An initial additive seasonal ARMA model satisfied conventional diagnostic requirements, but displayed ARCH effects. Of the (G)ARCH models assessed, the most parsimonious was a simple [ARCH-lag1, GARCH-lag1] model (Table 1). Although the individual GARCH term was nominally non-significant, there was a highly significant (P=0.0001) test of joint significance of the ARCH and GARCH parameters. There was no advantage of either t or general error distribution in the development of the (G)ARCH models.

(ii) Expected mortality: trend stationarity was demonstrated by rejection of the null-hypothesis of existence of a unit-root by the DF-GLS and ADF tests (with the trend option) and de-trending (linear regression of expected mortality against time) yielded residuals (also not evidencing a unit-root) for subsequent formal analysis. A simple (multiplicative) seasonal autoregressive model was generated with no evidence of ARCH effects (Table 1). Although an ARMA(1,1) model satisfied model diagnostic tests, the multiplicative seasonal AR model was favoured on clinical grounds.

Table 1. Parameters for GARCH (estimated from raw mortality series; de-trended linear model residuals) and ARMA (estimated from the expected mortality series; de-trended linear model residuals) models

Both the GARCH and ARMA models were considered parsimonious and the de-trended signals for each model were within 3SE limits of respective model predictions (Figure 7). The residuals from both the formal GARCH and ARMA models (mean: 0(0.0423) and 0(0.0257) respectively) satisfied multiple criteria of Gaussian white-noise and were within residual-EWMA control limits (default values of “optimal exponential coefficient” in Stata™ statistical software; 3SE control limits; λ = 0.0001 for both series; Figure 8). To address any potential inertial problems consequent upon the small λ, control limits were also established for projected 1 (λ=0.16), 2 (λ=0.42) and 3 (λ=0.71) SD increments of the mean of the GARCH residuals; the latter were within these control limits (Figure 9).

thumbnailFigure 7. De-trended series (navy line) generating the GARCH (upper panel) and ARMA (lower panel) models with one-step-ahead forecast control limits (3SE); upper control limit (red line), lower control limit (green line).

thumbnailFigure 8. “Optimal” residual-EWMA control chart (3 SE control limits): upper control limit (red line), lower control limit (green line), residuals (navy line); for GARCH (upper panel) and ARMA (lower panel) model residuals, respectively.

thumbnailFigure 9. Residual-EWMA control charts (3 SE control limits) for projected 1 (upper left panel), 2 (upper right pane) and 3 (lower left panel) SD increase of residual mean of the GARCH model (for raw mortality series); model residuals (navy line), upper control limit (red line), lower control limit (green line).

Discussion

The current analysis of monthly mortality series confirms the existence of autocorrelation and seasonality in both the raw and expected series at the individual ICU level, avoiding any potential confounding at the aggregate level due to Simpson’s paradox. We thus concur with the findings of Alwan [13,50] and Bisgaard and Kulahci [51], who documented the pervasiveness of autocorrelation in a variety of series, industrial and non-industrial. We also established that out-of-control signalling of the raw mortality series with respect to both 3 standard error risk-adjusted and RA-EWMA control limits was not evident with analysis of the residuals from the GARCH time series model. Thus the identification of (G)ARCH processes is an important issue for SPC [35].

As our focus was directed to an understanding of the underlying data-generating process [45] and the performance of the RA-EWMA control limits under conditions of autocorrelation, we deemed it appropriate to also subject the expected series, from which the control limits for the raw mortality series were established, to formal time-series estimation. Not surprisingly, as the underlying mortality estimates from a random coefficient model are obligatorily “smoothed” (see also Figure 1), no ARCH effects, representing “volatility”, were demonstrated and a relatively simple seasonal autoregressive model was established (Table 1). As the EWMA is based upon an ARIMA(0,1,1), that is an integrated moving average process [45,52], it has been applied to autocorrelated data [15], although the majority of studies have used relatively simple non-seasonal autoregressive models (AR(1) or AR(2)) with fixed λ (usually 0.2, which is the default for the SPC model of Statgraphics software). Residual-EWMA charts, in the context of time series modelling, would appear to be more robust than EWMA applied to the original (autocorrelated) data [53,54]. Reynolds and Lu have recommended that under autocorrelation “…traditional control chart methodology should not be applied without modification…” [55] and Human et al. have recently sounded a cautionary note about the robustness of the conventional EWMA [56].

For the classical SPC model, a process is in control if the mean and standard deviation estimate remain within prescribed control limits [57], usually three-sigma; that is, for a normally distributed series, 99.7% of observations should lie within the limits [58] and there a probability of 0.0027 that any point exceeds the control limit [32,59]. However this definition does not necessarily entail the formal time-series notion of stationarity (strict or weak), where the requirement for stationarity is that the first two moments (mean and variance [45,50]) and the autocorrelation function are time-invariant, albeit a stationary processes may be auto-correlated [60]. In the industrial/engineering sphere, practitioner response to process autocorrelation [61] was to embrace a time-series paradigm and apply SPC methods to the residuals of a formal time-series model [45,62], albeit there were different tactical approaches [15,63]; or to develop modified control limit schemes [64-66]. It is instructive to note that the non-model based EWMAST chart (EWMA chart for stationary processes [66]), recommended by Winkel and Zhang [11], pre-supposes a stationary (not “in-control”) process. In a systematic review of the application of statistical process control in healthcare, Thor et al. [5] adduced only one literature reference [67] and a calendar year 2003 monograph which discussed autocorrelation in medical series. As argued by Alwan and Roberts [45], systematic non-random patterns in series make separation of the classic common and special causes difficult, as departures from control, nominally traceable to special causes, are confounded by autocorrelation and, in the current series, seasonality. Two further concerns were raised by the authors; first, the undue emphasis placed upon normality and the (erroneous) assumption that “approximate normality” implies a state of statistical control; and second, in the presence of a well-fitting time series model with residuals consistent with white-noise (“randomness”), it is “…futile to search for departures from statistical control and their corresponding special causes…”. The latter caution resonates with the current finding of frequent signalling of the raw mortality series compared with in-control residuals from an apposite time series model; with respect to the error process, such signalling represents false positivity [13].

Cook and co-workers, “…explicitly compare[d] EWMA(observed) and EWMA(predicted) …[with] thresholds around the EMWA(predicted)…”, employing the EWMA (λ = 0.005-0.020) to “…effectively attenuate noise in the data and smooth an erratic but unbiased risk model” [25], although no criteria of “erratic” were provided. Smoothed control limits for the expected series were also utilised in a review paper by Cook et al. ([68]) and Pilcher et al. ([69], λ = 0.005), albeit the data structure differed; sequential plotting of each patient admission versus monthly mortality rates in the current paper. Our focus and methodology were different, in that we were concerned to both understand and formally model the “noise in the data”. This being said, in the current series, the smoothed EWMA (λ = 0.51) raw series (Figure 4) was demonstrated to signal using 3 standard error expected mortality control limits.

The sophistication of time-series modelling in standard statistical software packages makes the formal analyses of the current study feasible; in particular, automated routines for application of time series models [70]. However, for the application of appropriate SPC to mortality series from multiple ICUs in a data-base, there are unresolved statistical issues [71,72]. From the perspectives of this study, a multivariate approach may be established using more conventional estimators (multivariate GARCH [73] and vector autoregression models [74]) or by newly described hierarchical/functional time series [75,76].

Conclusions

The underlying data generating process of monthly mortality series at the ICU level displayed autocorrelation and seasonality, with volatility evident in the raw mortality series. Failure to accommodate these characteristics by SPC measures resulted in false-positive signalling. A time series approach to SPC, using residual control charts, would appear to resolve such issues.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

The study was conceived, designed, (data)-analysed, written and critically revised jointly by both authors (JLM, PJS). Both authors read and approved the final manuscript.

Acknowledgements

ANZICS Centre for Outcome and Resource Evaluation (CORE) of the Australian and New Zealand Intensive Care Society (ANZICS):

Australian and New Zealand Intensive Care Society, Carlton, Victoria 3053, Australia.

References

  1. Montgomery DC: Quality Improvement in the Modern Business Environment. In Introduction ot Statistical Quality Control. 7th edition. Edited by Montgomery DC. Hoboken, NJ: John Wiley & Sons, Inc; 2013:3-47. OpenURL

  2. Woodall WH: The Use of control charts in health-care and public-health surveillance.

    J Qual Technol 2006, 38:89-104. OpenURL

  3. Benneyan JC, Lloyd RC, Plsek PE: Statistical process control as a tool for research and healthcare improvement.

    Qual Saf Health Care 2003, 12:458-464. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Mohammed MA, Worthington P, Woodall WH: Plotting basic control charts: tutorial notes for healthcare practitioners.

    Qual Saf Health Care 2008, 17:137-145. PubMed Abstract | Publisher Full Text OpenURL

  5. Thor J, Lundberg J, Ask J, Olsson J, Carli C, Haerenstam KP, Brommels M: Application of statistical process control in healthcare improvement: systematic review.

    Qual Saf Health Care 2007, 16:387-399. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  6. Collins G, Jibawi A, McCulloch P: Control chart methods for monitoring surgical performance: a case study from gastro-oesophageal surgery.

    Ejso 2011, 37:473-480. PubMed Abstract | Publisher Full Text OpenURL

  7. Cook DA, Steiner SH, Cook RJ, Farewell VT, Morton AP: Monitoring the evolutionary process of quality: risk-adjusted charting to track outcomes in intensive care.

    Crit Care Med 2003, 31:1676-1682. PubMed Abstract | Publisher Full Text OpenURL

  8. Duclos A, Voirin N, Touzet S, Soardo P, Schott AM, Colin C, Peix JL, Lifante JC: Crude versus case-mix-adjusted control charts for safety monitoring in thyroid surgery.

    Qual Saf Health Care 2010, 19:1-4. OpenURL

  9. Kirkham JJ, Bouamra O: The use of statistical process control for monitoring institutional performance in trauma care.

    J Trauma 2008, 65:1494-1501. PubMed Abstract | Publisher Full Text OpenURL

  10. Rogers CA, Reeves BC, Caputo M, Ganesh JS, Bonser RS, Angelini GD: Control chart methods for monitoring cardiac surgical performance and their interpretation.

    J Thorac Cardiovasc Surg 2004, 128:811-819. PubMed Abstract | Publisher Full Text OpenURL

  11. Winkel P, Zhang NF: Statistical development of Quality in Medicine. Chichester, West Sussex: John Wiley & Sons Ltd; 2007. OpenURL

  12. Montgomery DC: Other univariate statistical process monitoring and control techniques. In Introduction ot Statistical Quality Control. 7th edition. Edited by Montgomery DC. Hoboken, NJ: Wiley; 2013:448-508. OpenURL

  13. Alwan LC: Effects of autocorrelation on control chart performance.

    Commun Stat Theory Methods 1992, 21:1025-1049. Publisher Full Text OpenURL

  14. Berthouex PM, Hunter WG, Pallesen L: Monitoring sewage-treatment plants - some quality-control aspects.

    J Qual Technol 1978, 10:139-149. OpenURL

  15. Montgomery DC, Mastrangelo CM: Some statistical process-control methods for autocorrelated data.

    J Qual Technol 1991, 23:179-193. OpenURL

  16. Wardell DG, Moskowitz H, Plante RD: Control charts in the presence of data correlation.

    Manag Sci 1992, 38:1084-1105. Publisher Full Text OpenURL

  17. Winkel P, Zhang NF: Serial correlation of quality control data - on the use of proper control charts.

    Scand J Clin Lab Invest 2004, 64:195-203. PubMed Abstract | Publisher Full Text OpenURL

  18. Winkel P, Zhang NF: Control charts for autocorrelated data. In Statistical development of Quality in Medicine. Edited by Winkel P, Zhang NF. Chichester, West Sussex: Wiley; 2007:92-110. OpenURL

  19. Moran JL, Solomon PJ, Adult Database Management Committee (ADMC) of the Australian and New Zealand Intensive Care Society (ANZICS): Conventional and advanced time series estimation: application to the Australian and New Zealand Intensive Care Society (ANZICS) adult patient database, 1993-2006.

    J Eval Clin Pract 2011, 17:45-60. PubMed Abstract | Publisher Full Text OpenURL

  20. Moran J, Solomon P: Mortality and Intensive Care volume in ventilated patients, 1995-2009, in the Australian and New Zealand bi-national adult patient intensive care database.

    Crit Care Med 2012, 40:800-812. PubMed Abstract | Publisher Full Text OpenURL

  21. Stow PJ, Hart GK, Higlett T, George C, Herkes R, McWilliam D, Bellomo R: Development and implementation of a high-quality clinical database: the Australian and New Zealand Intensive Care Society adult patient database.

    J Crit Care 2006, 21:133-141. PubMed Abstract | Publisher Full Text OpenURL

  22. ANZICS Centre for Outcome and Resource Evaluation (CORE) of the Australian and New Zealand Intensive Care Society (ANZICS):

    APD Data Dictionary: Version 3.2.1 Updated February 2012.

    [http://www.anzics.com.au/core/data-collection-tools webcite], Accessed June 30th 2012

    OpenURL

  23. Knaus WA, Wagner DP, Draper EA, Zimmerman JE, Bergner M, Bastos PG, Sirio CA, Murphy DJ, Lotring T, Damiano A: The APACHE III prognostic system. Risk prediction of hospital mortality for critically ill hospitalized adults.

    Chest 1991, 100:1619-1636. PubMed Abstract | Publisher Full Text OpenURL

  24. Rabe-Hesketh S, Skrondal A: Multilevel and Longitudinal Modeling Using Stata. 2nd edition. College Station, TX: Stata Press; 2008. OpenURL

  25. Cook DA, Coory M, Webster RA: Exponentially weighted moving average charts to compare observed and expected values for monitoring risk-adjusted hospital indicators.

    BMJ Qual Saf 2011, 20:469-474. PubMed Abstract | Publisher Full Text OpenURL

  26. Shiskin J: Decomposition of economic time series.

    Science 1958, 128:1539-1546. PubMed Abstract | Publisher Full Text OpenURL

  27. R Development Core Team: R : A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2012.

    [http://www.R-project.org webcite]

    OpenURL

  28. Stoffer D:

    Applied Statistical Time Series Analysis (“astsa”): R package (V 1.1).

    [http://www.stat.pitt.edu/stoffer/tsa3/ webcite], Accessed 12th August 2012

    PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  29. Montgomery DC: Cumulative sum and exponentially weighted moving averaged control charts. In Introduction ot Statistical Quality Control. 7th edition. Edited by Montgomery DC. Hoboken, NJ: Wiley; 2013:413-447. OpenURL

  30. StataCorp: Time Series Manual: Release 12. College Station, TX: StataCorp LP; 2011. OpenURL

  31. Reynolds MR, Stoumbos ZG: Comparisons of some exponentially weighted moving average control charts for monitoring the process mean and variance.

    Technometrics 2006, 48:550-567. Publisher Full Text OpenURL

  32. Montgomery DC: Methods and Philosophy of Statistical Process Control. In Introduction ot Statistical Quality Control. 7th edition. Edited by Montgomery DC. Hoboken, NJ: Wiley; 2013:187-233. OpenURL

  33. Zhang NF: Statistical control for autocorrelated data.

    Proc Soc Photo Opt Instrum Eng 1999, 3742:65-70. OpenURL

  34. StatPoint Technologies Inc:

    STATGRAPHICS Centurion XVI.1. 2012.

    Warrenton, Virginia USA

    OpenURL

  35. Fang Y, Zhang J: Performance of control charts for autoregressive conditional heteroscedastic processes.

    J Appl Stat 1999, 26:701-714. Publisher Full Text OpenURL

  36. Schipper S, Schmid W: Control charts for GARCH processes.

    Nonlinear Anal-Theor 2001, 47:2049-2060. Publisher Full Text OpenURL

  37. Shumway RH, Stoffer DS: ARIMA models. In Time Series Analaysis and Its Applications With R Examples. third edition. Edited by Shumway RH, Stoffer DS. New York, NY: Springer Science+Business Media, LLC; 2011:83-172. OpenURL

  38. Wang S-H, Hafner C: Estimating autocorrelations in the presence of deterministic trends.

    J Time Ser Econom 2011., 3

    Article4; http://www.degruyter.com/view/j/jtse.2011.3.2/jtse.2011.3.2.1022/jtse.2011.3.2.1022.xml webcite

    OpenURL

  39. Nelson CR, Plosser CR: Trends and random walks in macroeconmic time series : some evidence and implications.

    J Monet Econ 1982, 10:139-162. Publisher Full Text OpenURL

  40. Pierce DA: Trend and autocorrelation.

    Commun Stat 1975, 4:163-175. OpenURL

  41. Shumway RH, Stoffer DS: Additional time domain topics. In Time Series Analaysis and Its Applications With R Examples. third edition. Edited by Shumway RH, Stoffer DS. New York, NY: Springer Science+Business Media, LLC; 2011:267-318. OpenURL

  42. Tsay RS: Asset volatility and Volatility models. In An Introduction to Analysis of Financial Data with R. Edited by Tsay RS. Hoboken, NJ: Wiley; 2013:176-241. OpenURL

  43. Karlsson S:

    ARMADIAG: Stata module to compute post-estimation residual diagnostics for time series. 2009.

    [http://econpapers.repec.org/scripts/search.asp?ft=armadiag webcite], Accessed May 2009

    OpenURL

  44. Engle R: GARCH 101: the use of ARCH/GARCH models in applied econometrics.

    J Econ Perspect 2001, 15:157-168. Publisher Full Text OpenURL

  45. Alwan LC, Roberts HV: Time-series modeling for statistical process-control.

    J Bus Econ Stat 1988, 6:87-95. OpenURL

  46. Koehler AB, Marks NB, O’Connell RT: EWMA control charts for autoregressive processes.

    J Oper Res Soc 2001, 52:699-707. Publisher Full Text OpenURL

  47. Tsay RS: An Introduction to Analysis of Financial Data with R. John Wiley & Sons, Inc: Hoboken, NJ; 2013. OpenURL

  48. Kuha J: AIC and BIC: comparisons of assumptions and performance.

    Sociol Method Res 2004, 33:188-229. Publisher Full Text OpenURL

  49. StataCorp: tssmooth exponential -Single-exponential smoothing. In Time Series Manual: Release 12. College Station, TX; 2011:477-484. OpenURL

  50. Alwan LC, Roberts HV: The problem of misplaced control limits.

    Appl Stat-J Roy St C 1995, 44:269-278. Publisher Full Text OpenURL

  51. Bisgaard S, Kulahci M: Quality quandaries: the effect of autocorrelation on statistical process procedures.

    Qual Eng 2005, 17:481-489. Publisher Full Text OpenURL

  52. Box G, Narasimhan S: Rethinking statistics for quality control.

    Qual Eng 2010, 22:60-72. Publisher Full Text OpenURL

  53. Apley DW, Lee HC: Robustness comparison of exponentially weighted moving-average charts on autocorrelated data and on residuals.

    J Qual Technol 2008, 40:428-447. OpenURL

  54. Lu CW, Reynolds MR: EWMA control charts for monitoring the mean of autocorrelated processes.

    J Qual Technol 1999, 31:166-188. OpenURL

  55. Reynolds MR, Lu CW: Control charts for monitoring processes with autocorrelated data.

    Nonlinear Anal-Theor 1997, 30:4059-4067. Publisher Full Text OpenURL

  56. Human SW, Kritzinger P, Chakraborti S: Robustness of the EWMA control chart for individual observations.

    J Appl Stat 2011, 38:2071-2087. Publisher Full Text OpenURL

  57. Vasilopoulis AV, Stamboulis AP: Modification of control chart limits in the presence of data correlation.

    J Qual Technol 1978, 10:20-30. OpenURL

  58. Shahian DM, Williamson WA, Svensson LG, Restuccia JD, DAgostino RS: Applications of statistical quality control to cardiac surgery.

    Ann Thorac Surg 1996, 62:1351-1358. PubMed Abstract | Publisher Full Text OpenURL

  59. Tennant R, Mohammed MA, Coleman JJ, Martin U: Monitoring patients using control charts: a systematic review.

    Int J Qual Health Care 2007, 19:187-194. PubMed Abstract | Publisher Full Text OpenURL

  60. Kirchgassner G, Wolters J: Introduction to Modern Times Series Analysis. Berlin: Springer; 2008. OpenURL

  61. Lwin T: Parameter estimation in first-order autoregressive model for statistical process monitoring in the presence of data autocorrelation.

    J Stat Plan Infer 2011, 141:2556-2575. Publisher Full Text OpenURL

  62. Roberts HV, Tsay RS: Making control charts more effective by time series analysis: three illustrative applications.

    Commun Stat-Theory 1996, 25:2767-2796. Publisher Full Text OpenURL

  63. Runger GC: Assignable causes and auto correlation: control charts for observations or residuals?

    J Qual Technol 2002, 34:165-170. OpenURL

  64. Apley DW: Time series control charts in the presence of model uncertainty.

    J Manuf Sci E-T Asme 2002, 124:891-898. Publisher Full Text OpenURL

  65. Lee HC, Apley DW: Improved design of robust exponentially weighted moving average control charts for autocorrelated processes.

    Qual Reliab Engng Int 2011, 27:337-352. Publisher Full Text OpenURL

  66. Zhang NF: A statistical control chart for stationary process data.

    Technometrics 1998, 40:24-38. Publisher Full Text OpenURL

  67. Solodky C, Chen HG, Jones PK, Katcher W, Neuhauser D: Patients as partners in clinical research - a proposal for applying quality improvement methods to patient care.

    Med Care 1998, 36:AS13-AS20. PubMed Abstract | Publisher Full Text OpenURL

  68. Cook DA, Duke G, Hart GK, Pilcher D, Mullany D: Review of the application of risk-adjusted charts to analyse mortality outcomes in critical care.

    Crit Care Resusc 2008, 10:239-251. PubMed Abstract OpenURL

  69. Pilcher DV, Hoffman T, Thomas C, Ernest D, Hart GK: Risk-adjusted continuous outcome monitoring with an EWMA chart: could it have detected excess mortality among intensive care patients at Bundaberg Base Hospital?

    Crit Care Resusc 2010, 12:36-41. PubMed Abstract OpenURL

  70. Hyndman RJ, Khandakar Y: Automatic time series forecasting: the forecast package for R.

    J Stat Softw 2008, 27:1-22. OpenURL

  71. Bottle A, Aylin P: Predicting the false alarm rate in multi-institution mortality monitoring.

    J Oper Res Soc 2011, 62:1711-1718. Publisher Full Text OpenURL

  72. Marshall T, Mohammed MA, Rouse A: A randomized controlled trial of league tables and control charts as aids to health service decision-making.

    Int J Qual Health Care 2004, 16:309-315. PubMed Abstract | Publisher Full Text OpenURL

  73. Bauwens L, Laurent S, Rombouts JV: Multivariate GARCH models: a survey.

    J Appl Econ 2006, 21:79-109. Publisher Full Text OpenURL

  74. Pan X, Jarrett JE: Why and how to use vector autoregressive models for quality control: the guideline and procedures.

    Qual Quant 2012, 46:935-948. Publisher Full Text OpenURL

  75. de Silva A, Hyndman RJ, Snyder R: The vector innovations structural time series framework: a simple approach to multivariate forecasting.

    Stat Model 2010, 10:353-374. Publisher Full Text OpenURL

  76. Hyndman RJ, Ahmed RA, Athanasopoulos G, Shang HL: Optimal combination forecasts for hierarchical time series.

    Comput Stat Data An 2011, 55:2579-2589. Publisher Full Text OpenURL

Pre-publication history

The pre-publication history for this paper can be accessed here:

http://www.biomedcentral.com/1471-2288/13/66/prepub