Density Estimation and Wavelet Thresholding via Bayesian Methods: A Wavelet Probability Band and Related Metrics Approach to Assess Agitation and Sedation in ICU Patients

1. Introduction

A wave is usually defined as an oscillating function that is localized in both time and frequency. A wavelet is a “small wave”, which has its energy concentrated in time providing a tool for the analysis of transient, non-stationary, or time-varying phenomena [1, 2]. Wavelets have the ability to allow simultaneous time and frequency analysis via a flexible mathematical foundation. Wavelets are well suited to the analysis of transient signals in particular. The localizing property of wavelets allows a wavelet expansion of a transient component on an orthogonal basis to be modelled using a small number of wavelet coefficients using a low pass filter [3]. This wavelet paradigm has been applied in a wide range of fields, such as signal processing, data compression and image analysis [4 -10].

Typically agitation-sedation cycling in critically ill patients involves oscillations between states of agitation and over-sedation, which is detrimental to patient health, and increases hospital length of stay [11-14]. The goal of the research specifically in reference [14] was to develop a physiologically representative model that captures the fundamental dynamics of the agitation-sedation system. The resulting model can serve as a platform to develop and test semi-automated sedation management controllers that offer the potential of improved agitation management and reduce length of stay in the intensive care unit (ICU). A minimal differential equation model to predict or simulate each patient’s agitation-sedation status over time was presented in [14] for 37 ICU patients, and was shown to capture patient A-S dynamics. Current agitation management methods rely on subjective agitation assessment and an appropriate sedation input response from recorded at bedside agitation scales [15, - 19]. The carers then select an appropriate infusion rate based upon their evaluation of these scales, their experience and intuition [20]. This process is depicted in Figure 1 (see [14]). Recently a more refined A-S model, which utilised kernel regression with an Epanechnikov kernel and better captured the fundamental agitation-sedation (A-S) dynamics was formulated [12, 13].

A secondary aim of this chapter is to test the feasibility of wavelet statistics to help distinguish between patients whose simulated A-S profiles were “close” to their mean profile versus those for whom this was not the case (i.e. their simulated profiles are not “close” to their actual recorded profiles). This chapter builds on a preliminary study [21] to assess wavelet signatures for modelling ICU agitation-sedation profiles, so as to, as in this chapter, evaluate “closeness” or “discrimination” of simulated versus actual A-S profiles with respect to wavelet scales - as recently analysed using DWT and wavelet correlation methods in [29] (see also [22]-[24]). The recent work of Kang et al. [29] investigated the use of DWT signatures and statistics on the simulated profiles derived in [12] and [13], to test for commonality across patients, in terms of wavelet (cross) correlations. Another earlier application of this approach was the study of historical Australian flowering time series [22], where it was established that wavelets add credibility to the use of phenological records to detect climate change. This study was also recently expanded and reported by Hudson et al. [23, 24] (see also references [25-28]).

Figure 1.

Diagram of the feedback loop employing nursing staff’s feedback of subjectively assessed patient agitation through the infusion controller (diagram is sourced from [14]).

The density function is very important in statistics and data analysis. A variety of approaches to density estimation exist. Indeed the density estimation problem has a long history and many solutions [30, 31, 32]. A large body of existing literature on nonparametric statistics is devoted to the theory and practice of density estimation [32-36]. The local character of wavelet functions is the basis for their inherent advantage over projection estimators – specifically that wavelets are straightforward and well localized in both space and frequency. The relevant estimation methods belong to the class of so-called projection estimators, as introduced by [36] or their non-linear modifications. Section 3 traces, in brief, the development of some basic methods used in density estimation. We then link these and apply wavelet methods for (density) function estimation to the ICU data of [29].

In this chapter the density is estimated using wavelet shrinkage methods, as based on Bayesian methods. Specifically the minimax estimator is used to obtain a patient specific wavelet tracking coverage index (WTCI). All values of the WTCI are obtained using Bayesian wavelet thresholding, and are shown to differentiate between poor versus good tracking. A Bayesian approach is also suggested in this chapter by which to assess a parametric A-S model – this by constructing a wavelet probability band (WPB) for the proposed model and then checking how much the nonparametric regression curve lies within the band. The wavelet probability band (WPB) is shown to provide a useful tool to measure the comparability between the patient’s simulated and recorded profiles. Moreover, the density profile is then successfully used to define and compute two numerical measures, namely the average normalized wavelet density (ANWD) and relative average normalized wavelet density (RANWD) – both measures of agreement between the recorded infusion rate and simulated infusion rate. Our WPB method is shown to be a good tool for detecting regions where the simulated infusion rate (model) performs poorly, thus providing ways to help improve and distil the deterministic A-S model. The so-called Wavelet Time Coverage Index (WTCI) developed is analogous to the metrics based on a kernel based probability band of [13, 14]. The research in [29] and that formulated in this chapter have successfully developed novel quantitative measures based on wavelets for the analysis of A-S dynamics.

2. Density estimation using wavelet smoothing

In order to apply wavelets to various function estimation problems, it is useful to examine some of the existing techniques in use. This provides a useful lead in to a discussion of wavelet methods for (density) function estimation, since typical techniques can be modified in a straightforward manner for use in a wavelets approach and for a subsequent wavelets based analysis.

In exploratory data analysis it is important to have an idea of the shape of the data distribution, whereby interesting features become evident. For example, in describing the shape of the data, by a histogram, we easily obtain an overall feel for the data. Specialised versions of histograms that can be constructed using the Haar wavelet basis are now discussed in brief. Important theoretical properties of this estimator are discussed further in [37]. The Haar wavelet approach and histogram leads naturally to density estimators with smoother wavelet bases and lend themselves to histogram estimators, as we require.

Given the Haar scaling function, as on the left hand side of Equation (1), and then applying the usual dilation and translation gives,

$\varphi \left(t\right)=\left\{\begin{array}{c}1,\text{if}0\le t<1\\ 0,\text{otherwise}\end{array}\right\}\Rightarrow {\varphi }_{j,k}\left(t\right)=\left\{\begin{array}{l}{2}^{j/2},2^{-j}k\le t<2^{-j}\left(k+1\right)\\ 0,\text{otherwise}\end{array}\right\}.$E1

We can then count the number of data points that lie within a particular interval, say, $\left[2^{-j}k,2^{-j}\left(k+1\right)\right)$ using the quantity $2^{-j/2}{\displaystyle \sum _{i=1}^n{\varphi }_{j,k}\left(t_i\right)}.$

Now for any t∈ R and j∈Z,

$2^{-j}\left[2^{j}t\right]\le t<2^{-j}\left(2^{j}t+1\right),$E2

where [t] denotes the greatest integer function of t, the number of data points that lie in the same interval as any real number t can be computed by

$2^{-j}{\displaystyle \sum _{i=1}^n{\varphi }_{j,\left[2^{j}t\right]}\left(t_i\right)}={\displaystyle \sum _{i=1}^n\varphi \left(2^j{t}_i-2^{j}t\right)}.$E3

The histogram density estimator with origin 0 and bins of width $2^{-J}$ is given by

${\tilde{f}}_J\left(t\right)=\frac{1}{n}{2}^{J/2}{\displaystyle \sum _{i=1}^n{\varphi }_{J,\left[2^{J}t\right]}\left(t_i\right)}.$E4

This estimator can be regarded as being the best estimator of the density f on the approximation space V_J, where V_J is defined as length N/2^J vector scaling coefficients associated with averages on a scale of length 2^J = 2λ_J. Construction of histograms using the Haar basis, then leads to more general wavelet density estimators. The decomposition algorithm can be applied to Equation (1) and the resultant histogram can be written in terms of the Haar wavelets as follows:

${\tilde{f}}_J\left(x\right)={\displaystyle \sum _k{\tilde{c}}_{j_0,k}{\Phi }_{j_0,k}}\left(x\right)+{\displaystyle \sum _{j=j_0+1}^{J-1}{\displaystyle \sum _k{\tilde{d}}_{j^{\text{'}},k}{\Psi }_{j^{\text{'}},k}}}\left(x\right).$E5

The Haar-based histograms are given in Figure 2 (for level 1, 2, 3, and 4) for Patient 12 and Patient 18, with the simulated infusion rate (light) and the recorded infusion rate (dark) shown. Figure 2 shows a similar distribution between each patient’s recorded and simulated infusion rates - skewed right for both patients (P12 and P18). Correspondingly Figure 3 presents the simulated data and recorded A-S data of Patient 2 and Patient 27. Each patient’s simulated and recorded series are clearly from a differing distribution type to each other for these poor trackers (P2 and P27) (Figure 3). Figures 2-3 are clearly more informative than the histogram where the former density estimates are based on the Haar wavelet basis. These graphical comparisons allow us to visualize differences in the distribution between the poor and good tracking patients in ICU.

Estimating density functions using smooth wavelets can be performed in the same way for any orthogonal series. This estimation procedure, which is a natural application of wavelets, results from a straightforward extension of the Haar-based histogram approach [30]. The same approach used to estimate a density in terms of the Haar basis above can thus also be used with smooth wavelet bases, as we now illustrate.

Let ϕ and ψ be an orthogonal scaling function and mother wavelet pair that generates a series of approximating spaces ${\left\{V_j\right\}}_{j\in \mathbb{Z}}.$, then $f\left(x\right)$, which is a square integrable density function, is

$f\left(x\right)={\displaystyle \sum _k{c}_{j_0,k}{\varphi }_{j_0,k}\left(x\right)}+{\displaystyle \sum _{j>j_0}{\displaystyle \sum _k{d}_{j,k}{\psi }_{j,k}\left(x\right)}},$E6

where j₀ represents a coarse level of approximation. Haar coefficients are estimated using

${\widehat{c}}_{j,k}=⟨f^{~},{\varphi }_{j,k}⟩=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\varphi }_{j,k}\left(X_i\right)}$E7

${\widehat{d}}_{j,k}=⟨f^{~},{\psi }_{j,k}⟩=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\psi }_{j,k}\left(X_i\right)}.$E8

From Equations (7) and (8) above, the wavelet estimator for f at level J ≥ j₀ is given by

$\begin{array}{l}{\widehat{f}}_J\left(x\right)={\displaystyle \sum _k{\widehat{c}}_{j_0,k}{\varphi }_{j_0,k}\left(x\right)}+{\displaystyle \sum _{j>j_0}{\displaystyle \sum _k{\widehat{d}}_{j,k}{\psi }_{j,k}\left(x\right)}}\\ \text{ = }{\displaystyle \sum _k{\widehat{c}}_{J,k}{\varphi }_{J,k}\left(x\right).}\end{array}$E9

The smoothing parameter in Equation (9) is the index J of the highest level considered. Smooth wavelet-based density estimates are plotted in Figure 4 for levels 4, 6, and 8, using Patient 4’s simulated infusion and recorded infusion rate data (via Daub (4) which denotes the Daubechies wavelet filter of length 4). We sampled 2048 (=2¹⁰) data points without loss of any generality from the original data of Patient 4.

Figure 5 shows Patient 29’s smooth wavelet-based density estimates (using Daub (4)). Figures 4 and 5 indicate that Patent 4 (P4) is potentially a poor tracker and Patient 29 (P29) a good tracker, given that the original and wavelet smooth densities are similar in P29, but not for P4. Note that when level J is increased, abrupt jumps disappear but can also lead to over-smoothing and loss of information.

Figure 2.

Haar-based histogram of the simulated series (light) and recorded empirical A-S series (dark) for varying resolution levels for two “good trackers”: Patient 12 (left side, 4 plots) and Patient 18 (right side, 4 plots).

Figure 3.

Haar-based histogram of the simulated series (light) and recorded series (dark) for varying resolution levels for two “poor trackers”: Patients 27 (right side, 4 plots) and 2 (left side, 4 plots).

Figure 4.

Smooth wavelet-based density estimates for P4’s recorded data (light) and simulated data (dark) using the Daubechies wavelet (Daub4) with sub- sample N=2048 and for different choices of J.

To quantify the relationship between the two variables $(x_i,y_i)$, we can employ the standard regression model as follows,

$y_i=f(x_i)+{\epsilon }_i,i=1,\mathrm{...},n,$E10

where the ${\epsilon }_i$’s are independent and identically distributed $\text{N}\left(0,{\sigma }^2\right)$ random variables. It will be assumed that the design points $\left(x_1,\mathrm{...},x_n\right)$ are equally spaced, and, without further loss of generality, that they lie on the unit interval: $x_i=1/n,i=1,\mathrm{...},n.$ Our approach constitutes projecting the raw estimator f onto the approximating space V_J, for any choice of the smoothing parameter J, which represents a linear estimation approach. In contrast to this, the approach reference [38] offers a non-linear wavelet based approach to this problem of nonparametric regression. The approach in [38] begins with computing the DWT of the data $y_i,$ by generating a new data set of empirical wavelet coefficients with which to represent the underlying regression function f.

The estimation procedure in [38] has three main steps as follows: First, transform the data $y_i$ to the wavelet domain by applying a DWT. If d is the DWT of f, and $d^{\prime }={\left[{c^{\prime }}_{0,0},{d^{\prime }}_{0,0},\mathrm{...},{d^{\prime }}_{J-1,2^J-1}\right]}^T$ the vector of empirical coefficients, we then have a sequence of wavelet coefficients $d^{\prime }=\text{d+}{\epsilon }^{\prime }\text{,}$ where ${\epsilon }^{\prime }$ is a vector of n independent $\text{N}\left(0,{\sigma }^2\right).$ In the second step, the true coefficients d are estimated by applying the thresholding rule to the empirical coefficients $d^{\prime }$ to obtain estimates, $\tilde{d}$. Finally, the sampled function values f are estimated by applying the inverse DWT (IDWT) to obtain $\tilde{f}=W^T\tilde{d}$, where $W^T$ is the transpose of an orthonormal n×n matrix. We can then represent the DWT as the sum

$\tilde{f}\left(x\right)={\tilde{c}}_{0,0}\varphi \left(x\right)+{\displaystyle \sum _{j=0}^J{\displaystyle \sum _{k=0}^{2^j-1}{\tilde{d}}_{j,k}{\psi }_{j,k}\left(x\right)}}$E11

Figure 5.

Smooth wavelet-based density estimates for P29’s recorded (light) and simulated data (dark) using the Daubechies wavelet (Daub4) with N=2048 and for different choices of J.

This procedure, as formulated in [38], is schematised below:

Scheme 1.

Schema of DWT procedural steps.

A technique for selective wavelet reconstruction similar to this general approach was proposed in [39] in a study to remove random noise from magnetic resonance images. The technique is further developed from a statistical point of view in [38] by framing selective wavelet reconstruction as a problem in multivariate normal decision theory. From [38] estimating an unknown function involves including only coefficients larger than some specific threshold value. A large coefficient is taken to mean that it is large in absolute value. Choosing an excessively large threshold will make it difficult for a coefficient to be judged significant and be included in the reconstruction, resulting in over smoothing. On the other hand, a very small threshold allows many coefficients to be included in the reconstruction, resulting in under-smoothed estimates.

Two methods of global wavelet thresholding were proposed in [38], namely the universal threshold and the minimax threshold method. The wavelet assumption of a dyadic length of the time series is not always true. A natural approach would then be to pre-condition the original data set, so as to obtain a set of values of length 2^J for some positive integer J. The resulting pre-conditioned data is then plugged directly into any standard DWT. We observed that most of the 37 ICU patient data is not to the power of two. One obvious remedy was to pad the series with values and increase its length to the next power of two. There are several choices for the value of these padded coefficients. The approach adopted in this chapter was to pad with zeros so as to increase the size of the data set to the next larger power of two, or some other higher composite number, and then apply the DWT. The minimax estimator approach with soft thresholding, as applied to the simulated infusion profile of Patient 2, for example, yielded the profile in Figure 6.

3. New non Bayesian wavelet based metrics for tracking (WTCI, ANWD, RANWD)

Based on the development of density estimation via wavelet smoothing discussed in section 2, specifically equations (9) to (11), we now derive three new wavelet based, but non Bayesian metrics for tracking, namely the Wavelet Time Coverage Index (WTCI) (section 3.1), the Average Normalized Wavelet Density (ANWD) and the Relative Average Normalized Wavelet Density (RANWD) (section 3.2).

3.1. Numerical approach 1: Wavelet Time Coverage Index (WTCI)

The most commonly used criterion to obtain a successful wavelet estimator of the signal $\tilde{y}$ in estimating $y$ is the mean square error (MSE) [40]. In this chapter we devise a variant based on the development of the smoothed recorded infusion. This then lays the foundation for the development of our Wavelet Time Coverage Index (WTCI). The WTCI is a quantitative parameter indicating how well the patient’s simulated infusion represents their average recorded infusion profile over the entire time series. Our approach uses wavelet coefficients on a scale by scale basis.

The WTCI is defined as follows:

$\text{WTCI}=\left\{1-\frac{{\displaystyle {\sum }_{j,k}\left|{\tilde{d}}_{j,k}-d_{j,k}\right|}}{{\displaystyle {\sum }_{j,k}{\tilde{d}}_{j,k}}}\right\}\times 100$E12

where ${\tilde{d}}_{j,k}$ is given in Equation (11) and $d_{j,k}$ is the DWT of f in Equation (10). A WTCI of 100% represents perfect tracking, which arises when the DE simulated infusion profile is identical to that of the wavelet smoothed infusion profile.

Figure 6.

Minimax estimator applied to Patient 2’s simulated profile. The thick line represents the wavelet threshold estimator of the simulated infusion rate and the thin line that of the recorded infusion data. A soft thresholding rule was used to obtain all estimates.

Figure 7 presents the box and whisker plot for values of the WTCI from bootstrap [56] realizations (per patient). Each box and whisker [57] in Figure 7 displays two main components of information. First, the median represents a measure of how well the agitation-sedation simulation models the recorded infusion profile on average. Second, the spread of the box and whisper provides an indication of how reliable that particular WTCI median is per patient. Further details regarding interpretation of the WTCI are given in section 5.1.

3.2. Numerical approach 2: Average Normalized Wavelet Density (ANWD) and the Relative Average Normalized Wavelet Density (RANWD)

We now propose wavelet analogues of the AND and RAND diagnostics developed in [12] and [13]. The density profile is used to compute the numerical measures of ANWD and RANWD, so that objective comparisons of model performance can be made across different patients. The ANWD value for the simulated infusion rates is the average of these normalized density values over all time points for a given patient. Similarly, the RANWD value for the smoothed infusion rate is obtained by superimposing the smoothed values by using the first cumulant from a normal posterior distribution onto the same density profile, after which RANWD can be readily computed.

Figure 7.

Box and whisker plot of the WTCI index for each of the 37 patients.

Let $y_t=\left\{y_{1,}{y}_2,\mathrm{...},y_n\right\}$ be the output data produced by a proposed model. These are often called the simulated data. Define the average normalized wavelet density (ANWD) of $y_t$ as

$\text{ANWD}\left(y_t\right)=\frac{1}{n}{\displaystyle \sum _{t=1}^n\frac{{\tilde{f}}_t\left(y_t\right)}{\text{max}\left({\tilde{f}}_t\right)}}$E13

where $\text{max}\left({\tilde{f}}_t\right)$denotes the maximum value of the wavelet density function ${\tilde{f}}_t$, which is estimated by wavelet smoothing via Equation (9) at time t. Thus, ANWD is an average of normalized densities, where each component in the sum is the value of ${\tilde{f}}_t$ at $y_t$ normalized by $\text{max}\left({\tilde{f}}_t\right)$. At time t the normalized wavelet density equals 1 when $X_t$ coincides with the point where ${\tilde{f}}_t$ is maximum. An infusion profile that coincides with the maximum wavelet density at every time point would therefore have ANWD equal to 1. Whereas the value of ANWD for an infusion profile distant from the high-density regions would approach 0. Finally, ANWD ($y$) is calibrated using the ANWD from the wavelet smoothed recorded infusion data, denoted by $\tilde{y}$ giving the relative average normalized wavelet density (RANWD):

$$\text{RANWD}=\frac{\text{ANWD}\left(y\right)}{\text{ANWD}\left(\tilde{y}\right)}$$E14

RANWD indicates the value of ANWD ($y$) relative to a typical realisation in the form of $\tilde{y}$ from the density profile. Therefore, the RANWD statistic estimates how probabilistically alike the model outputs are to the smoothed data, and hence the degree of comparability between the model (simulated) and the actual (recorded) data. A RANWD of 0.6 implies that the model outputs are 60% similar, on average, to the wavelet smoothed data. Greater similarity means higher values of RANWD for the given patient under investigation.

4. Wavelet thresholding via Bayesian methods

Bayesian wavelet shrinkage methods are discussed in section 4.1, and are subsequently used to develop a novel 90% wavelet probability band (WPB 90%) per patient (section 4.2). A 90% wavelet probability band is constructed for each of the 37 patient profiles, and the time and duration of any deviations from the wavelet probability band is recorded. A WPB 90% value of 70% implies that for at least 70% (time under ICU observation), the estimated mean value of the recorded infusion rate, for a given patient, lies within the 90% confidence interval of its wavelet probability band. For illustration, we refer the reader to Figure 8 which shows the WPB 90% curves for 4 patients: two good trackers (Patients 8 and 25) and two poor trackers (Patients 9 and 34).

4.1. Brief mathematical background

Recall the regression equation (Equation (10)) for an observed data vector $y_1,y_{2,\mathrm{...},}{y}_n$ satisfying

$y_i=f\left(x_i\right)+{\epsilon }_i,i=1,\mathrm{...},n,$

where the ${\epsilon }_i$’s are independent and identically distributed $\text{N}\left(0,{\sigma }^2\right)$ random variables, assuming that $\left(x_1,\mathrm{...},x_n\right)$ are fixed points.

We now consider a method to approximate the posterior distribution of each f (x_i), using the same prior utilised by the BayesThresh method of [41] and [42]. Posterior probability intervals of any nominal coverage probability can be calculated accordingly. For Haar wavelets, the scaling function and mother wavelet are $\varphi \left(t\right)=I\left(0\le t<1\right)$ and $\psi \left(t\right)=I\left(0\le t<1/2\right)-I\left(1/2\le t<1\right)$, respectively, where $I(\cdot )$ is the indicator function. Clearly the square of the Haar wavelet is just the Haar scaling function, ${\psi }_{j,k}^2\left(t\right)=2^{j/2}{\varphi }_{j,k}\left(t\right),{\psi }_{j,k}^3\left(t\right)=2^j{\psi }_{j,k}\left(t\right)$; ${\psi }^3\left(t\right)=\psi \left(t\right)\text{ and }{\psi }^2\left(t\right)={\psi }^4\left(t\right)=\varphi \left(t\right)$ and ${\psi }_{j,k}^4\left(t\right)=2^{3j/2}{\varphi }_{j,k}\left(t\right).$ All these terms can be included in a modified version of the IDWT algorithm which incorporates scaling function coefficients. By the development in [43] we approximate a general wavelet ${\psi }_{j_0,0}^r,\left(0\le j_0\le J-m\right)$, by

${\psi }_{j_0,0}^r~{\displaystyle \sum _t^{}{e}_{j_0-m,l}}{\varphi }_{j_0-m,l}\left(t\right)$E15

for r = 2, 3, 4, where m is a positive integer. The choice of m follows below, since scaling functions (instead of wavelets), as the span of the set of scaling functions at a given level j, are the same as that of the sum of $\varphi \left(t\right)$ and the wavelets at levels 0, 1,…, j-1. Moreover, if scaling functions ${\varphi }_{j,k}\left(t\right)$are used to approximate some function h(t), and both ϕ and h have at least ν derivatives, then the mean squared error in the approximation is bounded by $C{2}^{-\nu j}$, where C is some positive constant, (see, for example reference [44]).

To approximate ${\psi }_{j_0,k}^r\left(t\right)$ for some fixed $j_0$, we simply compute $y_{j_0,0}^r\left(t\right)$ using the pyramid algorithm [45], then take the DWT and set the coefficients $e_{m_0,l}$ to be equal to the scaling function coefficients $e_{m_0,k}$ at level $m_0$, where $m_0=j_0+m$. Recall that the wavelets at level j are simply shifts of each other $y_{j,k}\left(t\right)=y_{j,0}\left(t-2^{-j}k\right),$ hence

$y_{j_0,k}^r\left(t\right)y_{j_0,0}^r\approx {\displaystyle \sum _t^{}{e}_{m_0,l-2^{m}k}}{\varphi }_{m_0,l}\left(t\right)$E16

Figure 8.

Wavelet Probability Bands (WPB 90%’s) (thin lines) with simulated infusion profile (thick line) for Patients 8, 25 (P8, P25: good trackers, LHS) and Patients 9, 34 (P9, P34: poor trackers, RHS).

As we are assuming periodic boundary conditions, the $e_{m_0,l}$ can be cycled periodically. Given the localised nature of wavelets, the coefficients $e_{m_0+1,l}$ employed to approximate ${\psi }_{j_0-1,0}^r\left(t\right)$ can be found by inserting $2^{m_0}$ zeros into the vector of $e_{m_0,l}$

$$e_{m+1,l}=\{\begin{array}{ccc}\sqrt{2}{e}_{m_0,l}& & l=0,\mathrm{...}\begin{array}{cc}{2}^{m_0-1}-1& \begin{array}{ccc}& & \end{array}\end{array}\\ 0& & l=2^{m_0-1},\mathrm{...},2^{m_0+1}-2^{m_0-1}-1\\ \sqrt{2}{e}_{m_0,l-2^{m_0}}& & l=2^{m_0+1}-2^{m_0-1},\mathrm{...},2^{m_0-1}-1\end{array}$$E17

The approximation in Equation (15) cannot be employed for wavelets at the m finest levels $$J-m,\mathrm{...},J-1$$. These wavelets are however, written in terms of both scaling functions and wavelets at the finest level of detail, level J-1, via the block-shifting method as delineated above.

From Equation (11) we have that $\left[f_i|y\right]$ is the convolution of the posteriors of the wavelet coefficients and the scaling coefficient given by,

$\left[f|y\right]=\left[c_{0,0}|{c^{″}}_{0,0}^{}\right]\varphi \left(t_i\right)+{\displaystyle \sum _j{\displaystyle \sum _k\left[d_{j,k}|{d^{″}}_{j,k}^{}\right]}}{\psi }_{j,k}\left(t_i\right).$E18

If X and Y are independent random variables and a and b are real constants, then

$${\kappa }_r\left(aX+b\right)=\{\begin{array}{c}a{\kappa }_1\left(X\right)+b,r=1\\ a^r{\kappa }_r\left(X\right),r=2,3,\cdots \end{array}$$E19

and we have by the additivity property

${\kappa }_r\left(X+Y\right)={\kappa }_r\left(X\right)+{\kappa }_r\left(Y\right),r\in \mathbb{Z}.$E20

for all r. Applying Equations (19) and (20) to Equation (17) shows $\left[f_i|y\right]$ (where $$f_i=f\left(x_i\right)$$) can be estimated from its cumulants as

${\kappa }_r\left(f_i|y\right)={\kappa }_r\left(c_{0,0}|{c^{″}}_{0,0}^{}\right){\varphi }_{j,k}^r\left(t_i\right)+{\displaystyle \sum _j{\displaystyle \sum _k{\kappa }_r\left(d_{j,k}|{d^{″}}_{j,k}^{}\right)}}{\psi }_{j,k}^r\left(t_i\right).$E21

The first cumulant ${\kappa }_1\left(y\right),$ is the mean of $y$, the second cumulant, ${\kappa }_2\left(y\right),$ is the variance of $y$, ${\kappa }_3\left(y\right)/{\kappa }_2^{3/2}\left(y\right)$ is the skewness, and ${\kappa }_4\left(y\right)/{\kappa }_2^2\left(y\right)+3$ is the kurtosis. Note that the third cumulant ${\kappa }_3\left(y\right)$ and the fourth cumulant ${\kappa }_4\left(y\right)$ are zero if $y$ is Gaussian. From Equations (16) and (21), we can now re-write the fourth cumulant as follows,

$$\begin{array}{l}{\kappa }_r\left(f_i|y\right)={\kappa }_r\left(c_{0,0}|c_{0,0}^{*}\right){\varphi }_{j,k}^r\left(t_i\right)+{\displaystyle \sum _j{\displaystyle \sum _k{\kappa }_r\left(d_{j,k}|d_{j,k}^{*}\right)}}{\psi }_{j,k}^r\left(t_i\right)\\ \text{ = }{\kappa }_r\left(c_{0,0}|c_{0,0}^{*}\right){\varphi }_{j,k}^r\left(t_i\right)+{\displaystyle \sum _{j,k}\left\{{\kappa }_r\left(d_{j,k}|d_{j,k}^{*}\right){\displaystyle \sum _l{e}_{j-3,l}{\varphi }_{j-3,l}\left(t_i\right)}\right\}},\\ \text{ = }{\displaystyle \sum _{j,k}{\rho }_{j,k}{\varphi }_{j,k}\left(t_i\right)}\end{array}$$E22

for $${\kappa }_r\left(y\right)$$ the r^th cumulant of $y$, and for suitable coefficients, $${\rho }_{j,k}$$ acquired via the IDWT algorithm which incorporates scaling function coefficients to assess this sum [41, 42, 46]. Bayesian wavelet regression estimates have thus been developed including priors on the wavelet coefficients $d_{j,k},$ which are updated by the observed coefficients ${d^{″}}_{j,k}$ to obtain posterior distributions $$\left[d_{j,k}|{\text{<msup>d <mo>″</mo></msup>}}_{j,k}\right]$$ (refer to Equation (18)). The $${\tilde{d}}_{j,k}$$ (point estimates) can then be computed from such posterior distributions and the Inverse Discrete Wavelet Transform (IDWT) used to estimate $$f\left(x_i\right).$$

The Bayesian wavelet shrinkage rules discussed in this section have used mixture distributions as priors on the coefficients to model a small proportion of the coefficients which contain substantial signal [41]. Indeed the BayesThresh method of [41] assumes independent priors on the coefficients,

$d_{j,k}~{\gamma }_{j}N\left(0,{\tau }_j^2\right)+\left(1-{\gamma }_j\right)\delta \left(0\right),j=0,\mathrm{...},J-1\text{; }k=0,\mathrm{...},2^j-1,$E23

where $0\le {\gamma }_j\le 1.0,\delta \left(0\right)$ is a point mass at zero and $$d_{j,k}$$ are independent. The hyper-parameters are assumed to be of the form ${\tau }_j^2=C_1{2}^{-\alpha j},$${\gamma }_j=\text{min}\left(1,C_2{2}^{-\beta j}\right)$ for non-negative constants C₁ and C₂ chosen empirically from the data and the α and β’s are selected by the user. The choice of α and β corresponds to choosing priors in certain Besov spaces [41] and incorporating prior knowledge about the smoothness of $$f\left(x_i\right)$$ into the prior. See reference [55].

4.2. New Bayesian 90% wavelet probability band metric for tracking (WPB 90%)

Bayesian wavelet shrinkage methods as discussed in section 4.1 can be used to create a wavelet probability band (WPB). A 90% wavelet probability band (WPB 90%) is constructed for each of the 37 patient profiles, and the time and duration of any significant deviations from the wavelet probability band is recorded. A WPB 90% value of 70% implies that for at least 70% of the time (of the time in ICU observation), the estimated mean value of the recorded infusion rate for a given patient lies within its 90% confidence interval of its wavelet probability band. Figure 8 shows the WPB for 4 patients: two good trackers (Patients 8 and 25) and two poor trackers (Patients 9 and 34). The circle symbol represents the hourly recorded infusion rate, the thin line represents the 90% WPB curve and the solid thick line represents the simulated profile (Figure 8). Brief spikes which may occur in the WPB bands are typical of wavelet regression methods. These spikes can be smoothed out by using different values for α and β, but this risks over-smoothing the data due to loss of information. According to [41], setting α=0.5 and β=1 is the best practical approach for Bayesian smoothing. Therefore we set α=0.5 and β=1 and employ Daubeches’ least asymmetric wavelet with eight vanishing moments, namely LaDaub (8), as this is a widely used wavelet and is applicable to a broad variety of data types.

While our WPB approach is graphically very useful (Figure 8), it is however useful to marry this with an objective numerical measure of how close the simulated infusion profile is to the empirical, recorded data. The percentage WPB cannot serve this purpose of objective quantification, because it quantifies visual proximity, by means of artificial hard boundaries, and ignores the fact that the in-band region does not have the same probabilistic significance everywhere. Thus wavelet density numerical metrics, namely ANWD and RANWD, comparing the model outputs to the recorded data were also developed using the posterior densities determined from the smoothed recorded data. The density profile is considered to be informative as unlike the wavelet probability band (WPB) it does discriminate between regions of high or low probability within a band.

5. WTCI based results

5.1. Choice of Wavelet filter and Bootstrap: WTCI

In order to judge the reliability of the wavelet time coverage index (WTCI) for a given patient’s infusion profile, the moving blocks bootstrap was utilized [56]. A total of 1000 bootstrap realizations were generated for each patient’s recorded infusion profiles. A wavelet time coverage index (WTCI), as defined in Equation (12), can then be evaluated for each bootstrap realization, providing a collection of 1000 values of the WTCI. The median WTCI and its standard error, SE, can then be reported for each patient using [57] (see Table 1, where a bold Patient no. indicates a poor tracker by the DWT, WCORR and WCCORR diagnostics in [29]). When the DWT is implemented via the pyramid algorithm [58], an important feature when analysing a given time, is the need to choose the appropriate wavelet filter (basis). The choice of a wavelet basis function is crucial for two reasons. First, the length of a DWT determines how well it approximates an ideal band-pass filter, which in turn dictates how well the filter is able to isolate features to specific frequency intervals. Secondly, as illustrated in the MODWT MRAs shown in [29], the wavelet basis function is being used to represent information contained in the time series of interest and should thus imitate its underlying features. A reasonable overall strategy is to use the shortest width of wavelet filters L= 4, 8 and longer wavelet filters L = 10, 12, as both choices give reasonable results in this ICU A-S application.

Patient no.	Data size; (Min)	Median; WTCI	SE
1	3601	78.55	0.538
2	6421	87.14	0.294
3	6541	87.85	0.106
4	4921	87.94	0.076
5	2941	88.89	0.103
6	5701	88.73	0.104
7	3901	84.78	0.324
8	10561	93.61	0.037
9	8581	93.28	0.046
10	20701	88.46	0.053
11	6721	92.46	0.085
12	8521	91.15	0.323
13	5161	91.37	0.091
14	3001	82.09	0.449
15	4981	92.09	0.072
16	13621	94.57	0.073
17	5941	90.27	0.086
18	4681	93.83	0.036
19	7921	96.34	0.012
20	9661	90.49	0.088
21	3721	83.07	0.685
22	9661	91.85	0.056
23	3481	85.07	0.300
24	8461	92.41	0.058
25	3841	93.44	0.082
26	3901	85.49	0.275
27	13441	93.66	0.039
28	12241	89.47	0.051
29	3241	89.3	0.262
30	3661	85.81	0.092
31	18301	94.34	0.022
32	15181	95.82	0.020
33	25261	95.63	0.036
34	8101	93.77	0.070
35	12721	87.64	0.018
36	3481	92.17	0.059
37	7501	90.79	0.079
Median		90.790	0.079
95%CI		(88.75, 92.39)	(0.058, 0.092)

Table 1.

Wavelet Time Coverage Index (WTCI) per patient. Data size (column 2) indicates the length of the patient’s A-S series. A bolded /shaded Patient indicates a poor tracker by the criteria in [29].

Table 1 presents the data size (time series length) and median bootstrapped WTCI and its standard error (SE) for each of the 37 patients studied. The poor trackers, as classified by the criteria in [29], are bolded in the first column. From Table 1, we note that some poor trackers (Patients 9, 11, 22, 27, 32, 34) have relatively high values of median WTCI and some good trackers (Patients 1, 14) have a low WTCI. The reason for this is padding, used as a reasonable solution to produce data of the size power of two, but this dilutes the signal near the end of the original data set, since the filters are not applied evenly. Hence multiplying by a signal element, constrained to have magnitude zero, is equivalent to omitting the filter coefficient, and then the orthogonality of the transform is not strictly maintained. To overcome this problem we change the current minutes driven length of the A-S data set to an hourly meta-meter, and then apply Bayesian wavelet shrinkage using a universal threshold as developed in section 5.2.

5.2. Alternative WTCI measure via Bayesian Wavelet Thresholding

In Section 5.1 the minimax estimator [59] was used to obtain patient specific WTCI measures (see also section 4.1). In this section we calculate an alternative WTCI using Bayesian wavelet thresholding [53, 44, 50, 41, 42]. Table 2 reports the WTCI measures, which are obtained by employing a Bayesian wavelet thresholding as computed by the WaveThresh software package in R [42] (adopting a LaDaub (8) filter).

The relative total dose, defined as the total drug dose delivered by the simulation (as a percentage of the total recorded drug dose) is also presented for each patient in Table 2. From Table 2 high values of median WTCI clearly indicate the validity of the A-S simulation for a given patient. The median WTCI value (across all patients) is 79.8% with a 95% interpolated confidence interval (CI) of (77.57%, 83.23%) and a range [71.9% to 88.9%]. This indicates some significant merit of the mathematical model of [14] and its physiological validity (Table 2).

Furthermore, the overall patient median relative total dose is 89.1% with a range [77.0% to 95.0%] indicating that the simulated and recorded total drug doses are similar, with the simulator consistently administering slightly less than 100% of the recorded actual sedative dose (Table 2). Slightly decreased levels as such are linked with the sudden-response nature of the recorded infusion profiles, in contrast to the consistent, smooth quality of the simulated infusion. These features are chiefly the result of the consistency of the computer implemented simulation in contrast to the inherent variation between different nurses’ assessment of patient agitation and appropriate feedback of sedation [14, 11, 12, 13].

Overall from Table 2 the values of the WTCI from the bootstrap realizations (per patient) have high median WTCI and a low spread per patient, indicating high reliability of the median WTCI and in turn of the A-S simulation. Larger spread indicates poor reliability, which may be caused by insufficient data, and also by the choice of wavelet filters, the chosen thresholding method, or by the model simulation method itself.

Patient no.	Median WTCI	SE	Relative total Dose (%)
1	74.0	0.378	87.8
2	77.4	0.461	85.9
3	76.6	0.264	86.1
4	79.2	0.650	86.6
5	79.3	0.324	90.3
6	84.6	0.697	87.5
7	75.5	0.083	81.0
8	84.3	0.381	92.1
9	88.1	0.359	89.1
10	74.9	0.175	90.2
11	72.0	0.395	87.8
12	83.8	0.261	86.5
13	82.8	0.166	90.2
14	76.6	0.248	88.5
15	86.3	0.145	90.7
16	84.8	0.431	90.4
17	83.3	0.430	85.2
18	84.2	0.161	95.0
19	86.5	0.341	91.0
20	80.5	0.303	90.4
21	72.1	0.498	87.8
22	77.5	0.171	89.5
23	79.8	0.673	91.7
24	85.1	0.127	89.9
25	88.9	0.075	91.2
26	82.8	0.061	88.5
27	75.7	0.663	87.5
28	75.0	0.202	90.6
29	73.4	0.442	77.0
30	83.3	0.224	94.5
31	82.6	0.146	90.0
32	82.6	0.271	91.2
33	78.5	0.430	90.4
34	85.9	0.193	89.0
35	74.1	0.361	88.1
36	79.1	0.293	81.6
37	78.6	0.283	88.9
75^thpercentile	84.3	0.43	90.4
Median	79.8	0.293	89.1
95% CI: Median	(77.57, 83.23)
25^th percentile	76.14	0.173	87.5

Table 2.

Alternative Wavelet Time Coverage Index (WTCI) summary for the 37 patients. A bolded/shaded patient indicates a poor tracker by the DWT WCORR and WCCORR threshold criteria in [29].

6. Bayesian WPB results

6.1. Wavelet probability band metric for tracking (WPB 90%)

Table 3 presents the time per patient that the simulated infusion rate lies within the 90% wavelet probability band (WPB 90%). Generally high values of WPB 90% are evident across most of the 37 patients (second column of Table 3). With the exception of thirteen (13) patients (Patients 2, 4, 7, 9, 10, 11, 21, 22, 27, 28, 32, 33, 34), all simulated infusion profiles lie within the wavelet probability band at least 75% of the time. These 13 patients were also all deemed to be “poor trackers” according to the WCORR and WCCORR diagnostics developed in [29]. The main reason for the reduced total time within the WPB for these 13 poor trackers seems to be their poor performance throughout the total length of the patient’s simulation. This feature is observed in Figure 8 for patient 9 and for patient 34, and indicates that the simulated infusion profile deviates from the recorded profile over some particular period, and takes some time before tending towards the recorded infusion rate again (see Patients 9 and 34 in Figure 8).

6.2. Comparison of WPB 90% with WTCI measures

Patients 25 and 34 will now be used to illustrate some of the concepts linking and differentiating the different wavelet measures in Table 2 and Table 3. Note that Patient 25 has the maximum WTCI of 88.9% with a high value of relative total dose of 91.2% (Table 2), and patient 25 exhibits low spread in the bootstrapped realization in Figure 7. Table 3 similarly shows that Patient 25’s simulated infusion rate lies within the 90% WPB for 89.06% of the time, indicative of good performance (or tracking), as is evident in the WPB plot (Figure 8). By contrast, Patient 34 has a high WTCI of 85.9%, a relative total dose of 89.0% (Table 2), and exhibits low spread in the bootstrapped realizations, but by contrast Patient 34 has a very low WPB 90% value of 48.44% (Table 3 and Figure 8).

Recall from [29] that Patient 25 is deemed to be a good tracker and Patient 34 a poor tracker by earlier DWT WCORR criteria. Hence whilst the WTCI values of Patient 25 and Patient 34 are both high, 88.9% and 85.9%, respectively; it is only the WPB 90% measure, and not the WTCI measure, that distinguishes between the tracking performance of Patient 25 (WPB 90% = 89.06%) and Patient 34 (WPB 90% = 48.44%). Patient 34’s simulation infusion rate is outside the WPB band for 51.56% of the time indicating that a large maximum departure time between the patient’s recorded and simulated infusion rate occurs in the ICU observation period (see [29]).

6.3. Comparison of WPB 90% with ANWD and RANWD measures

The higher the percentage that the simulated infusion profile lies within the wavelet probability band, for a given patient, the better the simulation model captures the specific dynamics of the agitation-sedation system for that patient. Patient-specific WPB 90% values are reported in Table 3. Columns 3 and 4 in Table 3 present the two novel and alternative performance measures of ANWD and RANWD per patient. Recall that the posterior distribution of the regression curve is used as the density for all patients as described in reference [14]. RANWD thus measures how probabilistically similar the model outputs are to the smoothed observed data, and hence is a measure of the degree of comparability between the simulated and the empirical A-S data. The density profile is used to compute the numerical measures of ANWD and RANWD, so that objective comparisons of model performance can be made across different patients. The ANWD value for the simulated infusion rates is the average of these normalized density values over all time points for a given patient. Similarly, the RANWD value for the smoothed infusion rate is obtained by superimposing the smoothed values by using the first cumulant from a normal posterior distribution onto the same density profile, after which RANWD can be readily computed.

Patient no.	WPB 90%	ANWD	RANWD
1	95.31	0.537	0.553
2	64.06	0.431	0.499
3	96.88	0.632	0.737
4	59.38	0.338	0.475
5	93.75	0.495	0.504
6	95.31	0.659	0.980
7	67.19	0.417	0.455
8	87.50	0.567	0.688
9	57.81	0.343	0.412
10	66.80	0.300	0.388
11	74.34	0.423	0.434
12	84.38	0.622	0.662
13	73.44	0.442	0.504
14	96.88	0.449	0.476
15	89.06	0.702	0.761
16	82.81	0.596	0.770
17	85.94	0.506	0.566
18	93.75	0.548	0.558
19	74.22	0.759	0.780
20	96.88	0.487	0.581
21	65.62	0.407	0.413
22	65.62	0.422	0.455
23	92.19	0.288	0.341
24	71.00	0.655	0.635
25	89.06	0.635	0.670
26	96.88	0.600	0.601
27	47.27	0.368	0.608
28	50.78	0.501	0.540
29	82.81	0.343	0.394
30	96.88	0.554	0.597
31	87.50	0.562	0.669
32	68.36	0.326	0.362
33	58.79	0.373	0.499
34	48.44	0.505	0.551
35	96.10	0.371	0.533
36	75.00	0.573	0.763
37	79.69	0.448	0.607
Min	47.27	0.288	0.341
Median	82.81	0.495	0.552
Max	96.88	0.759	0.981

Table 3.

Wavelet probability band (WPB 90%), ANWD and RANWD measures per patient. A bold patient no. indicates a poor tracker by the RANWD and WPB criteria (developed in section 6.4). (The 13 poor trackers are P2, P4, P7, P9, P10, P11, P21, P22, P27, P28, P32, P33 and P34).

An overall median RANWD of 0.552 (Table 3) with range [0.341 to 0.981] is an objective measure that supports the WPB 90% measures and visual clue of closeness based on the WPB (see Figure 8). It should be noted that as the model is deterministic, its outputs do not belong to the same probabilistic mechanism that generated the data, hence RANWD is an extremely stringent measure.

6.4. WPB 90% and RANWD criteria for poor tracking

Given the conservative nature of the RANWD metric, consistently high RANWD values close to 1.00 are not expected, even for a good simulation model. A reasonable and practical threshold for adequate model performance is RANWD ≥ 0.5, which suggests that the model outputs are more alike than not to the smoothed data. Justification for our 0.5 threshold for RANWD is given in this section, as is a threshold for WPB 90%. Poor trackers according to the metrics developed in this chapter are assumed to satisfy the following:

$$\text{RANWD}\le \text{0}\text{.5 }\text{and/or }\text{a }\text{WPB 90\% }\le \text{70\%}\text{.}$$

Thirteen (13) patients have a low RANWD (deemed below the threshold, RANWD $\le$0.50), namely Patients 2, 4, 7, 9, 10, 11, 14, 21, 22, 23, 29, 32 and 33 (Table 4). Specifically 9 of these patients have RANWD values from 0.412 to 0.499, and 4 have RANWD values between 0.341 and 0.394. Furthermore 11 patients have a WPB 90% value less or equal to 70% (Patients 2, 4, 7, 9, 10, 11, 27, 28, 32, 33, 34). Whilst Patients 14 and 23 have low RANWD values they exhibit very high WPB 90% (> 92%) values (like P29) (Table 3). Note that these 3 patients (P14, P23, P29) were also classified as good trackers according to the criteria developed earlier in [13], [14] and [29]; and as such, given their high WPB 90%, values, will be classified as good trackers in this chapter. Clearly the WBP 90% measure can help find patients (good trackers, say) who have elevated percentage time in the WPB ( (range 83% - 97%) for these 3 patients (P14, P23, P29) even though they exhibit relatively low RANWD values (range 0.34-0.47) (see column 3 Table 6).

Thirteen patients which satisfy our criterion for poor tracking (RANWD $\le$0.5, and/or a WPB 90% $\le$ 70%) are P2, P4, P7, P9, P10, P11, P21, P22, P27, P28, P32, P33 and P34 - all of whom, were also identified as poor trackers by the WCORR and WCCORR DWT diagnostics of [29]. Of these 13 poor trackers 8 have both lower than threshold WPB 90% and low RANWD, 3 exhibit low WPB 90% but above threshold RANWD > 0.5, and 2 exhibit low RANWD but above threshold WPB 90% (see Table 6). This indicates a significant and high agreement between the WPB and the RANWD (dichotomized) criteria (kappa = 0.6679); P (estimated kappa ≤ 0.40) = 0.025).

The resultant, RANWD and WPB 90% thresholds also provide very strong support for the DWT wavelet diagnostics derived in reference [29], in that of the 15 DWT based poor trackers identified in [29] (see Table 6), 13 of these also exhibit a low WPB (WPB 90% < 70%) and/or a low wavelet density based RANWD (RANWD $\le$0.5). Statistically speaking the wavelet probability band and density diagnostics developed in this chapter mirror the DWT based poor versus good classification of [29] (kappa = 0.87, p = 0.0001).

Our 13 poor trackers have WPB and density profiles (not all reported here) which have specific regions where the patient’s DE model did not appear to capture the observed A-S dynamics. In some scenarios, this may occur in the absence of a stimulus or when low drug concentrations coincide with an agitation level that is decreasing (but not close to zero), thus causing the patient’s agitation to remain at a constant non-zero level, despite their recorded infusion rate dropping to near zero.

6.5. Comparison of WPB 90% with Rudge’s physiological model [13] (AND, RAND measures)

The non-wavelets based and earlier performance metrics of AND and of RAND per patient are also given in Table 4. Table 4 thus provides a comparison between the WPB 90%, ANWD and RANWD and TIB, AND and RAND measures from Rudge’s Physiological Model [12, 13]. A highlighted patient is a poor tracker by the WCORR/WCCORR criteria in [29]. Table 4 along with Table 3, allows comparison between Rudge’s [13] values of AND and RAND, with our WPB model diagnostics (WPB) and our wavelet-based estimates of ANWD and RANWD. An underlined patient indicates a poor tracker by our RANWD and WPB criteria (see also Table 6).

Patient no.	WPB model			Rudge’s Physiological Model
	WPB 90%	ANWD	RANWD	TIB 90%	AND	RAND
1	95.31	0.537	0.553	96	0.51	0.62
2	64.06	0.431	0.499	90	0.53	0.66
3	96.88	0.632	0.737	97	0.70	0.83
4	59.38	0.338	0.475	93	0.56	0.62
5	93.75	0.495	0.504	97	0.60	0.80
6	95.31	0.659	0.980	95	0.70	0.84
7	67.19	0.417	0.455	67	0.33	0.43
8	87.50	0.567	0.688	90	0.45	0.59
9	57.81	0.343	0.412	89	0.49	0.62
10	66.80	0.300	0.388	53	0.27	0.34
11	77.34	0.423	0.434	59	0.31	0.38
12	84.38	0.622	0.662	96	0.61	0.77
13	73.44	0.442	0.504	85	0.37	0.45
14	96.88	0.449	0.476	95	0.48	0.56
15	89.06	0.702	0.761	95	0.45	0.60
16	82.81	0.596	0.770	91	0.44	0.57
17	85.94	0.506	0.566	91	0.61	0.72
18	93.75	0.548	0.558	92	0.55	0.68
19	74.22	0.759	0.780	90	0.50	0.66
20	96.88	0.487	0.581	91	0.53	0.65
21	65.62	0.407	0.413	95	0.53	0.72
22	65.62	0.422	0.455	83	0.35	0.45
23	92.19	0.288	0.341	95	0.72	0.85
24	71.10	0.655	0.635	91	0.43	0.54
25	89.06	0.635	0.670	86	0.50	0.66
26	96.88	0.600	0.601	92	0.68	0.88
27	47.27	0.368	0.608	84	0.39	0.49
28	50.78	0.501	0.540	76	0.34	0.44
29	82.81	0.343	0.394	90	0.38	0.45
30	96.88	0.554	0.597	97	0.63	0.82
31	87.50	0.562	0.669	74	0.40	0.51
32	68.36	0.326	0.362	74	0.38	0.50
33	58.79	0.373	0.499	67	0.28	0.36
34	48.44	0.505	0.551	84	0.43	0.55
35	96.10	0.371	0.533	70	0.38	0.46
36	75.00	0.573	0.763	83	0.52	0.64
37	79.69	0.448	0.607	92	0.53	0.59
Min	47.27	0.288	0.341	53	0.27	0.34
Median	82.81	0.495	0.552	90	0.49	0.60

Table 4.

Comparison between the WPB, ANWD and RANWD and TIB, AND and RAND from Rudge’s Physiological Model [12, 13]. A boxed patient is a poor tracker by the WCORR /WCCORR criteria in [29]. A shaded patient indicates a poor tracker by our RANWD and WPB criteria (13 patients: P2, P4, P7, P9, P10, P11, P21, P22, P27, P28, P32, P33 and P34).

6.6. Comparison of WPB, WTCI, ANWD and RANWD across poor versus good tracking groups

Table 5 gives summary statistics of the wavelet density based metrics (WPB, WTCI, ANWD and RANWD) for the poor versus good trackers (classified using the threshold criterion for WPB 90% $\le$ 70% and RANWD$\le$0.50). The poor trackers have significantly lower median values of WPB 90% (64.84% versus 87.50%) (p $\le$0.001); a significantly lower median value of WTCI (76.56% versus 82.79%) (p $\le$0.041); a significantly lower median value of ANWD (0.41 versus 0.55) (p $\le$0.001) and a significantly lower median value for RANWD (0.46 versus 0.59) (p $\le$0.001) compared to the good tracking group.

	WPB 90%	WTCI	ANWD	RANWD
Poor trackers
Min	47.27	71.95	0.03	0.36
Max	77.34	88.05	0.50	0.61
Range	30.07	16.01	0.21	0.25
Mean	61.56	77.97	0.39	0.47
95% CI of Mean	(65.79, 67.32)	(74.72, 81.22)	(0.36, 0.44)	(0.42, 0.51)
Median	64.84	76.56	0.41	0.46
95% CI of Median	(52.63, 67.09)	(74.89, 81.69)	(0.34, 0.43)	(0.41, 0.53)
Good trackers
Min	58.79	73.37	0.29	0.34
Max	96.88	88.94	0.76	0.98
Range	38.09	15.57	0.47	0.64
Mean	85.31	81.36	0.53	0.60
95% CI of Mean	(80.85, 89.77)	(76.62, 83.09)	(0.48, 0.57)	(0.55, 0.66)
Median	87.50	82.79	0.55	0.59
95% CI of Median	(82.72, 93.79)	(79.13, 84.09)	(0.45, 0.60)	(0.53, 0.67)
Kruskal-Wallis; test; P value	0.001	0.041	0.001	0.001

Table 5.

Summary of the wavelet based performance metrics: poor versus good trackers.

6.7. Poor trackers compared across 3 studies (references [29], [14] and [13])

Table 6 summarises the patient numbers of the poor trackers according to the criteria of four studies, including the research described in this chapter (i.e. Kang’s WPB diagnostics, see column 1). The four studies reported across columns 1, 4-6 in Table 6 are Kang’s WPB, WCORR diagnostics [29], Chase diagnostics [14] and the earlier Rudge diagnostics [13].

Kang WPB diagnostics (this chapter)	WPB ≤ 70%	RANWD$\le$0.50	Kang WCORR diagnostics [29]	Chase et al. [14] diagnostics	Rudge el al. [13] diagnostics
-			-	-	-
2	2	2	2	-	-
-	-	-	-	-	-
4	4	4	4	-	-
-	-	-	-	-	-
-	-	-	-	6	-
7	7	7	7	7	7
-	-	-	-	-	-
9	9	9	9	9	-
10	10	10	10	-	10
11	11	11	11	-	11
-	-	-	-	12	-
-	-	-	-	-	13
		14^Φ	-	-	-
-	-	-	-	-	-
-	-	-	-	-	-
-	-	-	-	17	-
-	-	-	-	-	-
-	-	-	-	-	-
-	-	-	-	-	-
21		21	21	21	-
22	-	22	22	-	22
		23^Φ	-	-	-
-	-	-	-	-	-
-	-	-	-	-	-
-	-	-	-	-	-
27	27		27	27	27
28	28		28	-	28
-	-	29^Φ	29	-	29
-	-	-	-	-	-
-	-	-	-	-	-
32	32	32	32	-	-
33	33	33	33	-	33
34	34	-	34	34	-
-	-	-	35	-	35
-	-	-	-	-	-
-	-	-	-	-	-
Total: N1= 13			Total: N2=15	Total: N3=8	Total: N4=10

Table 6.

Patient numbers of the poor trackers according to the criteria of 4 studies. ^Φ P14 and P23 have low RANWD values but high WPB 90% (> 92%) (like P29). P14, P23, P29 were classified as good trackers according to the criteria developed earlier in [13], [14] and [29], and as such, and given their high WPB 90%, are classified as good trackers by our Kang WPB diagnostics.

The resultant ANWD, RANWD, WTCI and WPB 90% thresholds also provide very strong support for the DWT wavelet diagnostics derived in reference [29], in that of the 15 DWT based poor trackers identified in [29], 13 also exhibit a low WPB (WPB 90% < 70%) and/or a low wavelet density based RANWD measure (RANWD $\le$0.5), and are likewise deemed to be poor trackers (Table 6). Statistically speaking the wavelet probability band and density diagnostics developed in this chapter mirror the DWT based criterion of [29] (kappa = 0.87, p = 0.0001). Indeed of the 13 patients assessed by our WPB and RANWD criteria to be poor trackers, all were likewise judged to be poor trackers by the earlier DWT WCORR and WCCORR criteria developed in reference [29] (see Table 6 and also see Tables 4-5 of [29]). This indicates perfect agreement between the RANWD threshold developed in this chapter and the earlier DWT WCORR and WCCORR based criteria for poor tracking in [29] (kappa = 1.00, p = 0.0000).

The performance metrics of AND and RAND and their patient specific values are given in Table 4, which along with Table 3 also allows comparison between Rudge’s [13] (AND and RAND) values with our WPB model diagnostics (WPB, ANWD, RANWD). Rudge’s Physiological Model [12, 13] found 10 of the 37 patients (27%) have values of RAND $\le$ 0.5, with 5 patients with 0.43< RAND<0.49, and 3 with 0.34 < RAND <0.38 (Tables 4 - 5). The model in [63] likewise found that 27 patients (73%) have RAND values greater than 0.57, with 10 poor trackers, with 6 RAND values ranging from 0.43 to 0.49, and with 3 patients exhibiting RAND values between 0.34 and 0.38. The main reason for the reduced total time within the WPB (and the non-significant WCORRs) for this minority group of 10 - 13 poor trackers (of the total 37 patients), is the consistently poor performance of the DE model throughout their total length of the A-S simulation. Of the 13 patients assessed by our RANWD and WPB criteria to be poor trackers, 7 patients (P7, P10, P11, P22, P27, P28, P33) were likewise judged to be poor trackers by the earlier Physiological Model of Rudge [13] (kappa = 0.30, p = 0.03). This shows significant agreement between the physiological Model [13] based criteria for poor tracking and the RANWD and WPB 90% thresholds formulated in this chapter (Table 6).

7. Discussion and conclusions

Agitation management via effective sedation management is an important and fundamental activity in the ICU. However, in clinical practice a lack of understanding of the underlying dynamics, combined with a lack of subjective assessment tools, makes effective and consistent clinical agitation management difficult [14, 12, 13]. The main goal of ICU sedation is to control agitation, while preventing over-sedation and over-use of drugs. Current clinical practice employs subjective agitation and sedation assessment scales, combined with medical staff experience and intuition, to deliver appropriate sedation. This approach usually leads to the administration of largely continuous infusions which lack a bolus-focused approach, and commonly results in either over sedation, or insufficient sedation [12, 13]. Several recent studies have emphasised the cost and health-care advantages of drug delivery protocols based on assessment scales of agitation and sedation. Table 7 gives an overview of recent ICU agitation studies, and provides a brief overview of the equations used for simulations of a patient’s A-S status and also of the methods derived in this chapter (and by other studies) with the aim of establishing the validity of the models in reflecting a patient’s true A-S status.

In this chapter, we successfully developed a density estimation approach via wavelet smoothing to assess the validity of deterministic dynamic A-S models. This wavelet density approach provided graphical assessment and numerical metrics (WTCI and WPB 90%, ANWD and RANWD) to assess the comparability between the modelled and the recorded A-S data per patient. Our new wavelet regression diagnostics identified 13 ICU patients (patients 2, 4, 7, 9, 10, 11, 21, 22, 27, 28, 32, 33 and 34) (out of 37 analysed) [29], whose simulated A-S profiles were poor indicators of their true A-S status, the remaining patients tracked exceptionally well.

All of these 13 poor trackers were also identified as poor trackers by the DWT measures derived in [29]. The WTCI and WPB 90% metrics derived in this chapter thus give strong support for the datawork in [29] and vice versa. Our wavelet regression diagnostics (WTCI, ANWD, RANWD, and WPB 90%) are thus valid for assessing control, as were the wavelet DWT, wavelet correlation (WCORR) and cross-correlation (WCCORR) measures derived in [29]. We have thus successfully assessed the patients A-S by the RANWD cut-point and also distinguished poor trackers, likewise identified by the DWT criteria of [29]. Ten of our 13 so-called poor trackers were also identified as poor trackers by either the kernel smoothing, tracking index and probability band approach of [14] and [13], respectively. Overall the various diagnostics strongly agree and confirm the value of A-S modelling in ICU. Our WPB method is also shown to be an excellent tool for detecting regions where the simulated infusion rate performs poorly, thus providing ways to help improve and distil the deterministic A-S model. The main reason for the reduced total time within the WPB for a minority group, of 13 (of 37) i.e. the poor trackers, is the consistently poor performance of the DE model throughout the total length of the simulation.

Wavelet modelling in this chapter and the earlier work of Kang [29] thereby demonstrate that the models of the recent A-S studies of [14, 11, 12, 13, 63] and of [64], are suitable for developing more advanced optimal infusion controllers. These offer significant clinical potential of improved agitation management and reduced length of stay in critical care. Further details are available in the recent PhD dissertation of Kang (circa 2012) [67]. The A-S time series profiles studied in this chapter are of disparate lengths (with a wide range [3,001 - 25,261] time points in minutes). Our approach is thus generalisable to any study which investigates the similarity or closeness of bivariate time series of, say, a large number of units (patients, households) and of time series of varying lengths and of possibly long length. This chapter demonstrates the value of wavelets for assessing ICU agitation-sedation deterministic models, and suggests new wavelet probability band and coverage diagnostics by which to mathematically assess A-S models. Future work will involve creation of singular spectrum analysis (SSA) based similarity indices following the development of Hudson and Keatley in [68]; and comparing results with the wavelet density-based indices developed in this chapter with the DWT metrics in [29] and in references [69] - [70].

Table 7.

Overview of Studies on ICU

C_c, C_p and C_e are, respectively, the drug concentrations (mg L⁻¹) in the central, peripheral and effect compartments;, U is the intravenous infusion rate; V_d, V_c, V_p and V_e, respectively, the volume of distribution, the distribution volumes (L) of the central, peripheral and effect compartments; A is an agitation index, S is the stimulus invoking agitation; K₁–K₃ are parameters relating to drug elimination and transport and K_ij the transfer rate (L min⁻¹) from compartment i to compartment j; K_CL the drug clearance (L min⁻¹); K_T the effect, and w₁ and w₂ are relative weighting coefficients of the stimulus and drug effect, respectively. Time is represented by t, and τ is the variable of integration in the convolution integral V_c, V_p and V_e, respectively, the distribution volumes (L) of the central, peripheral and effect compartments; U the intravenous infusion rate (mL min⁻¹); A an agitation index; S the stimulus invoking agitation; K_ij the transfer rate (L min⁻¹) from compartment i to compartment j; K_CL the drug clearance (L min⁻¹); K_T the effect time constant (min⁻¹); P^o and P^s are the concentrations of morphine and midazolam, respectively (mgmL⁻¹), where terms with superscript ‘o’ relate to the opioid morphine, and terms with superscript ‘s’ relate to the sedative midazolam. Time is represented by t (min), the variable of integration, and the terms w₁ and w₂ are the relative weights of stimulus and cumulative effect, representing the patient sensitivity. Finally, Ecomb is the combined pharmacodynamic effect of the individual effect site drug concentrations of morphine and midazolam determined using response surface modeling as defined in [66].