Open access

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

Written By

In Kang, Irene Hudson, Andrew Rudge and J. Geoffrey Chase

Submitted: 26 February 2012 Published: 06 February 2013

DOI: 10.5772/52434

From the Edited Volume

## Discrete Wavelet Transforms

Edited by Awad Kh. Al - Asmari

Chapter metrics overview

View Full Metrics

## 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]).

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,

ϕ(t)={1, if 0t<10, otherwise} ϕj,k(t)={2j/2, 2-jkt<2-j(k+1)0, otherwise}.           E1

We can then count the number of data points that lie within a particular interval, say, [2-jk, 2-j(k+1)) using the quantity 2-j/2i=1nϕj,k(ti).

Now for any tR and jZ,

2-j[2jt]t<2-j(2jt+1),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-ji=1nϕj,[2jt](ti)=i=1nϕ(2jti2jt).E3

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

f˜J(t)=1n2J/2i=1nϕJ,[2Jt](ti).E4

This estimator can be regarded as being the best estimator of the density f on the approximation space VJ, where VJ is defined as length N/2J vector scaling coefficients associated with averages on a scale of length 2J = 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:

f˜J(x)=kc˜j0,kΦj0,k(x)+j=j0+1J-1kd˜j',kΨj',k(x).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 {Vj}j., then f(x), which is a square integrable density function, is

f(x)=kcj0,kϕj0,k(x)+j>j0kdj,kψj,k(x),E6

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

c^j,k=f~,ϕj,k=1ni=1nϕj,k(Xi)E7
d^j,k=f~,ψj,k=1ni=1nψj,k(Xi).E8

From Equations (7) and (8) above, the wavelet estimator for f at level Jj0 is given by

f^J(x)=kc^j0,kϕj0,k(x)+j>j0kd^j,kψj,k(x)           = kc^J,kϕJ,k(x).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 (=210) 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.

To quantify the relationship between the two variables (xi, yi), we can employ the standard regression model as follows,

yi=f(xi)+εi, i=1,...,n,E10

where the εi’s are independent and identically distributed N(0, σ2) random variables. It will be assumed that the design points (x1,...,xn) are equally spaced, and, without further loss of generality, that they lie on the unit interval: xi=1/n, i=1,...,n. Our approach constitutes projecting the raw estimator f onto the approximating space VJ, 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 yi, 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 yi to the wavelet domain by applying a DWT. If d is the DWT of f, and d=[c0,0,d0,0,...,dJ1,2J1]T the vector of empirical coefficients, we then have a sequence of wavelet coefficients d=d+ε, where ε is a vector of n independent N(0,σ2). In the second step, the true coefficients d are estimated by applying the thresholding rule to the empirical coefficients d to obtain estimates, d˜. Finally, the sampled function values f are estimated by applying the inverse DWT (IDWT) to obtain f˜=WTd˜, where WT is the transpose of an orthonormal n×n matrix. We can then represent the DWT as the sum

f˜(x)=c˜0,0ϕ(x)+j=0Jk=02j1d˜j,kψj,k(x)E11

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 2J 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 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:

WTCI={1j,k|d˜j,kdj,k|j,kd˜j,k}×100E12

where d˜j,k is given in Equation (11) and dj,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 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.

Let yt={y1, y2,...,yn} be the output data produced by a proposed model. These are often called the simulated data. Define the average normalized wavelet density (ANWD) of yt as

ANWD(yt)=1nt=1nf˜t(yt)max(f˜t)E13

where max(f˜t)denotes the maximum value of the wavelet density function 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 f˜t at yt normalized by max(f˜t). At time t the normalized wavelet density equals 1 when Xt coincides with the point where 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 y˜ giving the relative average normalized wavelet density (RANWD):

RANWD=ANWD(y)ANWD(y˜)E14

RANWD indicates the value of ANWD (y) relative to a typical realisation in the form of 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 y1,y2,...,yn satisfying

yi=f(xi)+εi,       i=1,...,n,

where the εi’s are independent and identically distributed N(0,σ2) random variables, assuming that (x1,...,xn) are fixed points.

We now consider a method to approximate the posterior distribution of each f (xi), 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 ϕ(t)=I(0t<1) and ψ(t)=I(0t<1/2)I(1/2t<1), respectively, where I() is the indicator function. Clearly the square of the Haar wavelet is just the Haar scaling function, ψj,k2(t)=2j/2ϕj,k(t), ψj,k3(t)=2jψj,k(t); ψ3(t)=ψ(t) and ψ2(t)=ψ4(t)=ϕ(t) and ψj,k4(t)=23j/2ϕj,k(t). 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 ψj0,0r,(0j0Jm), by

ψj0,0r~tej0m,lϕj0m,l(t)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 ϕ(t) and the wavelets at levels 0, 1,…, j-1. Moreover, if scaling functions ϕj,k(t)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 C2νj, where C is some positive constant, (see, for example reference [44]).

To approximate ψj0,kr(t) for some fixed j0, we simply compute yj0,0r(t) using the pyramid algorithm [45], then take the DWT and set the coefficients em0,l to be equal to the scaling function coefficients em0,k at level m0, where m0=j0+m. Recall that the wavelets at level j are simply shifts of each other yj,k(t)=yj,0(t2jk), hence

yj0,kr(t)yj0,0rtem0,l2mkϕm0,l(t)E16

As we are assuming periodic boundary conditions, the em0,l can be cycled periodically. Given the localised nature of wavelets, the coefficients em0+1,l employed to approximate ψj01,0r(t) can be found by inserting 2m0 zeros into the vector of em0,l

em+1,l={2em0,l l=0,...2m011 0l=2m01,...,2m0+12m0112em0,l2m0l=2m0+12m01,...,2m011E17

The approximation in Equation (15) cannot be employed for wavelets at the m finest levels Jm,...,J1. 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 [fi|y] is the convolution of the posteriors of the wavelet coefficients and the scaling coefficient given by,

[f|y]=[c0,0|c0,0]ϕ(ti)+jk[dj,k|dj,k] ψj,k(ti).                E18

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

κr(aX+b)={aκ1(X)+b, r=1 arκr(X), r=2,3,E19

and we have by the additivity property

κr(X+Y)=κr(X)+κr(Y), r.E20

for all r. Applying Equations (19) and (20) to Equation (17) shows [fi|y] (where fi=f(xi)) can be estimated from its cumulants as

κr(fi|y)=κr(c0,0|c0,0)ϕj,kr(ti)+jkκr(dj,k|dj,k) ψj,kr(ti).    E21

The first cumulant κ1(y), is the mean of y, the second cumulant, κ2(y), is the variance of y, κ3(y)/κ23/2(y) is the skewness, and κ4(y)/κ22(y)+3 is the kurtosis. Note that the third cumulant κ3(y) and the fourth cumulant κ4(y) are zero if y is Gaussian. From Equations (16) and (21), we can now re-write the fourth cumulant as follows,

κr(fi|y)=κr(c0,0|c0,0*)ϕj,kr(ti)+jkκr(dj,k|dj,k*) ψj,kr(ti)               = κr(c0,0|c0,0*)ϕj,kr(ti)+j,k{κr(dj,k|dj,k*)lej3,lϕj3,l(ti)},               = j,kρj,kϕj,k(ti)E22

for κr(y) the rth cumulant of y, and for suitable coefficients, ρ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 dj,k, which are updated by the observed coefficients dj,k to obtain posterior distributions [dj,k|d j,k] (refer to Equation (18)). The d˜j,k (point estimates) can then be computed from such posterior distributions and the Inverse Discrete Wavelet Transform (IDWT) used to estimate f(xi).

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,

dj,k~γjN(0,τj2)+(1γj)δ(0), j=0,...,J1k=0,...,2j1,   E23

where 0γj1.0, δ(0) is a point mass at zero and dj,k are independent. The hyper-parameters are assumed to be of the form τj2=C12αj,γj=min(1,C22βj) for non-negative constants C1 and C2 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(xi) 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 75th percentile 84.3 0.43 90.4 Median 79.8 0.293 89.1 95% CI: Median (77.57, 83.23) 25th 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:

RANWD0.5 and/or WPB 90% 70%.

Thirteen (13) patients have a low RANWD (deemed below the threshold, RANWD 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 0.5, and/or a WPB 90% 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 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% 70% and RANWD0.50). The poor trackers have significantly lower median values of WPB 90% (64.84% versus 87.50%) (p 0.001); a significantly lower median value of WTCI (76.56% versus 82.79%) (p 0.041); a significantly lower median value of ANWD (0.41 versus 0.55) (p 0.001) and a significantly lower median value for RANWD (0.46 versus 0.59) (p 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-WallistestP 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≤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 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 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

Cc, Cp and Ce are, respectively, the drug concentrations (mg L−1) in the central, peripheral and effect compartments;, U is the intravenous infusion rate; Vd, Vc, Vp and Ve, 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; K1K3 are parameters relating to drug elimination and transport and Kij the transfer rate (L min−1) from compartment i to compartment j; KCL the drug clearance (L min−1); KT the effect, and w1 and w2 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 Vc, Vp and Ve, respectively, the distribution volumes (L) of the central, peripheral and effect compartments; U the intravenous infusion rate (mL min−1); A an agitation index; S the stimulus invoking agitation; Kij the transfer rate (L min−1) from compartment i to compartment j; KCL the drug clearance (L min−1); KT the effect time constant (min−1); Po and Ps are the concentrations of morphine and midazolam, respectively (mgmL−1), 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 w1 and w2 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].

## References

1. 1. GoupillaudPGrossmannAMorletJ1984Cycle-octave and related transforms in seismic signal analysis. Geoexploration, 23(1), 85-102.
2. 2. MorletJ1983Sampling theory and wave propagation. Issues in Acoustic Signal/Image Processing and Recognition, 1233261
3. 3. BarberSNasonG. PSilvermanB. W2002Posterior probability intervals for wavelet thresholding. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 64(2), 189-205.
4. 4. MallatS1989A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern analysis and Machine Intelligence, 11674693
5. 5. MeyerF. G2003Wavelet-based estimation of a semiparametric generalized linear model of fMRI time- series. IEEE Transactions on Medical Imaging, 22(3), 315-322.
6. 6. KumarPFoufoula-georgiouE1993A multicomponent decomposition of spatial rainfall fields. 1. Segregation of large- and small-scale features using wavelet transforms. Water Resources Research, 29(8), 2515-2532.
7. 7. KumarP1994Role of coherent structures in the stochastic-dynamic variability of precipitation. Journal of Geophysical Research-Atmospheres, 101(D21)(26), 393-404.
8. 8. DonohoD. L1995De-noising by soft-thresholding. IEEE Transactions on Information Theory, 41(3), 613-627.
9. 9. ChangS. GYuBVetterliM2000aSpatially adaptive wavelet thresholding based on context modeling for image denoising. IEEE Transactions on Image Processing, 9(9), 1522- 1531.
10. 10. ChangS. GYuBVetterliM2000bAdaptive wavelet thresholding for image denoising and compression. IEEE Transactions on Image Processing, 9(9), 1532-1546.
11. 11. RudgeA. DChaseJ. GShawG. MLeeDWakeG. CHudsonI. Let al2005Impact of control on agitation-sedation dynamics. Control Engineering Practice, 13(9), 1139-1149.
12. 12. RudgeA. DChaseJ. GShawG. MLeeD2006bPhysiological modelling of agitation-sedation dynamics. Medical Engineering and Physics, 28(1), 49-59.
13. 13. RudgeA. DChaseJ. GShawG. MLeeD2006aPhysiological modelling of agitation-sedation dynamics including endogenous agitation reduction. Medical Engineering Physics, 28(7), 629-638.
14. 14. ChaseJ. GRudgeA. DShawG. MWakeG. CLeeDHudsonI. Let al2004Modeling and control of the agitation-sedation cycle for critical care patients. Medical Engineering and Physics, 26(6), 459-471.
15. 15. FraserG. LRikerR. R2001bMonitoring sedation, agitation, analgesia, and delirium in critically ill adult patients. Crit Care Clin, 17(4), 967-987.
16. 16. JaarsmaA. SKnoesterHVan RooyenFBosA. P2001Biphasic positive airway pressure ventilation (pev+) in children. Crit Care Clin, 5(3), 174-177.
17. 17. RamsayM. ASavegeT. MSimpsonB. RGoodwinR1974Controlled sedation with alphaxalone-alphadolone. Br Med J, 2(920), 656-659.
18. 18. RikerR. RPicardJ. TFraserG. L1999Prospective evaluation of the Sedation-Agitation Scale for adult critically ill patients. Critical Care Medicine, 27(7), 1325-1329.
19. 19. SesslerC. NGosnellM. SGrapM. JBrophyG. MONealP. VKeaneK. A., et al2002The Richmond agitation-sedation scale: validity and reliability in adult intensive care unit patients. Am J Respir Crit Care Med, 166(10), 1338-1344.
20. 20. KressJ. PPohlmanA. SHallJ. B2002Sedation and analgesia in the intensive care unit. Am J Respir Crit Care Med, 166(8), 1024-1028.
21. 21. KangIHudsonI. LRudgeA. DChaseJ. G2005Wavelet signature of Agitation-Sedation profiles of ICU patients. In: Francis AR, Matawie KM, Oshlack A, Smyth GK (eds) Statistical Solutions to Modern Problems. 20th International Workshop on Statistical Modelling, Sydney, July 10-15, University of Western Sydney (Penrith), 2932961-74108-101-7
22. 22. KangIHudsonI. LKeatleyM. R2004Wavelet analysis in phenological research-the study of four flowering eucalypt species in relation to climate. International Biometrics Conference (IBC 04), July, Cairns, Australia.
23. 23. HudsonI. L2010Interdisciplinary approaches: towards new statistical methods for phenological studies, Climatic Change, 100(1), 143-171.
24. 24. HudsonI. LKeatleyM. R2010Phenological Research: Methods for Environmental and Climate Change Analysis: Springer Dordrecht. 523p.
25. 25. HudsonI. LKangIRudgeA. DChaseJ. GShawG. M2004Wavelet signatures of agitation and sedation profiles: a preliminary study. New Zealand Physics Engineering in Medicine (NZPEM) Christchurch New Zealand, Nov 2223
26. 26. HudsonI. LKeatleyM. RRobertsA. M. I2005Statistical Methods in Phenological Research. Proceedings of the Statistical Solutions to Modern Problems. Proceedings of the 20th International Workshop on Statistical Modelling, eds. AR Francis, KM Matawie, A Oshlack GK Smyth, Sydney, Australia, 10-15 July, 2592701-74108-101-7
27. 27. HudsonI. LKangIKeatleyM. R2010bWavelet Analysis of Flowering and Climatic Niche Identification. In I. L. Hudson IL, Keatley MR, editors. Phenological Research: Methods for Environmental and Climate Change Analysis. Springer Dordrecht. 361391
28. 28. HudsonI. LKeatleyM. RKangI2011Wavelet characterisation of eucalypt flowering and the influence of climate. Environmental & Ecological Statistics. 183513533
29. 29. KangIHudsonI. LRudgeAChaseJ. G2011Wavelet signatures and diagnostics for the assessment of ICU Agitation-Sedation protocols, In: Olkkonen H (editor) Discrete Wavelet Transforms, InTech. 321348978-9-53307-654-6
30. 30. OgdenR. T1997Essential Wavelets for Statistical Applications and Data Analysis. Boston: Birkhäuser.
31. 31. SilvermanB. W1986Density Estimation for Statistics and Data Analysis, New York: Chapman Hall, 4547
32. 32. WalnutD. F2004An Introduction to Wavelet Analysis: Birkhauser.
33. 33. ClineD. B. HEubankR. LSpeckmanP. L1995Nonparametric estimation of regression curves with discontinuous derivatives. Journal of Statistical Research, 291730
34. 34. DelouilleVFrankeJVon SachsR2001Nonparametric stochastic regression with design-adapted wavelets. Sankhy : The Indian Journal of Statistics, Series A, 63(3), 328-366.
35. 35. VidakovicB1999Statistical modelling by wavelets. New York: John Wiley Sons.
36. 36. CencovN. N1962Evaluation of an unknown distribution density from observations. Soviet Mathematics, 315991562
37. 37. EngelJ1990Density estimation with Haar series. Statistics Probability Letters, 9(2), 111-117.
38. 38. DonohoD. LJohnstoneI. M1994Ideal spatial adaptation via wavelet shrinkage. Biometrika, 81425455
39. 39. WeaverJ. BYansunXHealyD. M. JCromwellL. D1991Filterring noise from images with wavelet transforms. Magnetic Resonance in Medicine, 24288295
40. 40. DonohoD. L1995De-noising by soft-thresholding. IEEE Transactions on Information Theory, 41(3), 613-627.
41. 41. AbramovichFSapatinasTSilvermanB. W1998Wavelet thresholding via a Bayesian approach. Journal of the Royal Statistical Society. Series B, Statistical Methodology, 725749
42. 42. NasonG. P2008Wavelet Methods in Statistics with R: Springer Verlag.
43. 43. HerrickD. R. M2000Wavelet Methods for Curve and Surface Estimation. Thesis. University of Bristol, Bristol, U.K.
44. 44. VidakovicB1998Nonlinear Wavelet Shrinkage with Bayes Rules and Bayes Factors. Journal of the American Statistical Association, 93(441), 173-179.
45. 45. DaubechiesI1992Ten lectures on wavelets. Philadelphia: Society for Industrial and Applied Mathematics.
46. 46. BarberSNasonG. PSilvermanB. W2002Posterior probability intervals for wavelet thresholding. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 64(2), 189-205.
47. 47. HillI. DHillRHolderR. L1976AlgorithmA. SFitting Johnson curves by moments. Appl. Statist, 25180189
48. 48. HillI. D1976AlgorithmA. SNormal-Johnson and Johnson-Normal transformationns. Appl. Statist, 25190192
49. 49. JohnsonN. L1949Systems of frequency curves generated by methods of translation. Biometrika, 36149176
50. 50. ChipmanH. AKolaczykE. DMccullochR. E1997Adaptive Bayesian wavelet shrinkage. Journal of the American Statistical Association, 92(440).
51. 51. ChipmanH. AWolfsonL. J1999Prior elicitation in the wavelet domain. Lecture Notes in Statistics, 1418394
52. 52. ClydeMParmigianiGVidakovicB1998Multiple shrinkage and subset selection in wavelets. Biometrika, 85(2), 391-402.
53. 53. ClydeMGeorgeE. I2000Flexible empirical Bayes estimation for wavelets. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 62(4), 681-698.
54. 54. CrouseM. SNowakR. DBaraniukR. G1998Wavelet-based statistical signal processing using Hidden Markov models. IEEE Transactions on Signal Processing, 46(4), 886-902.
55. 55. ClydeMGeorgeE. I2000Flexible empirical Bayes estimation for wavelets. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 62(4), 681-698.
56. 56. ChipmanH. AWolfsonL. J1999Prior elicitation in the wavelet domain. Lecture Notes in Statistics, 1418394
57. 57. EfronBTibshiraniR. J1993An introduction to the Bootstrap, (57of Monographs on Statistics and Applied Probability).
58. 58. HettmanspergerT. PMckeanJ. W1998Robust Nonparametric Statistical Methods. London: Arnold.
59. 59. MallatS1989A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern analysis and Machine Intelligence, 11674693
60. 60. DonohoD. LJohnsonN. L1998Minimax estimation via wavelet shrinkage. Annals of Statistics, 26879921
61. 61. GencayRSelcukFWhitcherB2002An introduction to wavelets and other filtering methods in finance and economics: Academic Press.
62. 62. IKangI. LHudsonA. DRudgeJ. GChaseWavelet signatures and diagnostics for the assessment of ICU Agitation-Sedation protocols, for book entitled Discrete Wavelet Transforms-Biomedical Applications, InTech Publishers, http://www.intechweb.org,Vienna, Austria, 978953307654
63. 63. PercivalD. BWaldenA. T2000Wavelet Methods for Time Series Analysis. Cambridge, England: Cambridge University Press.
64. 64. RudgeA. DChaseJ. GShawG. MWakeG. C2003Improved agitation management in critically ill patients via feedback control of sedation administration. In: Proc World Congress on Medical Physics and Biomed Eng, August, Sydney, Australia.
65. 65. LeeD. SRudgeA. DChaseJ. GShawG. M2005A new model validation tool using kernel regression and density estimation. Computer Methods and Programs in Biomedicine, 807587
66. 66. WandM. PJonesM. C1995Kernel smoothing: Chapman Hall/CRC.
67. 67. MintoC. FSchniderT. WShortT. GGreggK. MGentiliniAShaferS. L2000Response surface model for anesthetic drug interactions. Anesthesiology, 92(6), 1603-1616.
68. 68. KangI2012Wavelets, ICA and Statistical Parametric Mapping: with application to agitation-sedation modelling, detecting change points and to neuroinformatics. PhD thesis University of Canterbury, Christchurch, New Zealand. 338 p.
69. 69. HudsonI. LKeatleyM. R2010Singular Spectrum Analysis: Climatic Niche Identification. In: Hudson IL, Keatley MR (eds) Phenological Research: Methods for Environmental and Climate Change Analysis, Springer, Dordrecht, 393424
70. 70. HudsonI. LKeatleyM. RKangI2011Wavelets and clustering: methods to assess synchronization. In: del Valle M, Muñoz R, Gutiérrez JM (eds) Wavelets: Classification, Theory and Applications. Nova Science Publishers. Chapter 5. 97124978-1-62100-252-9
71. 71. HudsonI. LKeatleyM. RKangI2011Wavelet Signatures of Climate and Flowering: Identification of Species Groupings In: Olkkonen H (ed) Discrete Wavelet Transforms Biomedical Applications, InTech, publishers. Vienna. 267296978-9-53307-654-6

Written By

In Kang, Irene Hudson, Andrew Rudge and J. Geoffrey Chase

Submitted: 26 February 2012 Published: 06 February 2013