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

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. 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. This wavelet paradigm has been applied in a wide range of fields, such as signal processing, data compression and image analysis.


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][12][13][14].The goal of the research specifically in reference [14] was to develop a physiologically representative model that captures the fundamental dynamics of the agitationsedation system.The resulting model can serve as a platform to develop and test semiautomated 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][26][27][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][33][34][35][36].The local character of wavelet functions is the basis for their inherent advantage over projection estimators -specifically that wavelets Discrete Wavelet Transforms -A Compendium of New Approaches and Recent Applications 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.

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, We can then count the number of data points that lie within a particular interval, say, 2 -j k, 2 -j (k + 1) ) using the quantity 2 -j/2 ∑ i=1 n ϕ j,k (t i ).
Now for any t∈ R and j∈Z, ( ) 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 The histogram density estimator with origin 0 and bins of width 2 -J is given by 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: ( ) ( ) ( ) 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).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 {V j } j∈ℤ ., then f (x), which is a square integrable density function, is , where j 0 represents a coarse level of approximation.Haar coefficients are estimated using ( ) .
From Equations ( 7) and ( 8) above, the wavelet estimator for f at level J ≥ j 0 is given by 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 10 ) 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 (x i , y i ), we can employ the standard regression model as follows, ( ) , 1,..., , where the ε i 's are independent and identically distributed N ( 0, σ 2 ) random variables.It will be assumed that the design points (x 1 , ..., x n ) are equally spaced, and, without further loss of generality, that they lie on the unit interval: x i = 1 / n, i = 1, ..., 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 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 ˜= W T d ˜, where W T is the transpose of an orthonormal n×n matrix.We can then represent the DWT as the sum This procedure, as formulated in [38], is schematised below: Wavelet Estimate: 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 undersmoothed 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.

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

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 1 100 where 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 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.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 y t = {y 1, y 2 , ..., y n } 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 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 ( ) ( ) 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.

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).
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 , respectively, where I ( ⋅ ) is the indicator function.Clearly the square of the Haar wavelet is just the Haar scaling function, . 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 ψ j 0 ,0 r , (0 ≤ j 0 ≤ J − m), by ( ) 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 ψ j 0 ,k r (t) for some fixed j 0 , we simply compute y j 0 ,0 r (t) 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 (t) = y j,0 ( t − 2 − j k ) , hence  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 ψ j 0 −1,0 r (t) can be found by inserting 2 m 0 zeros into the vector of e m 0 ,l 0 0 0 0 0 The approximation in Equation ( 15) cannot be employed for wavelets at the m finest levels J − m, ..., 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 f i | y is the convolution of the posteriors of the wavelet coefficients and the scaling coefficient given by, ( ) ( ) Density Estimation and Wavelet Thresholding via Bayesian Methods: A Wavelet Probability Band and Related… http://dx.doi.org/10.5772/52434 If X and Y are independent random variables and a and b are real constants, then and we have by the additivity property for all r.Applying Equations ( 19) and (20) to Equation (17) shows f i | y (where be estimated from its cumulants as The first cumulant κ 1 (y), is the mean of y, the second cumulant, κ 2 (y), is the variance of y, κ 3 (y) / κ 2 3/2 (y) is the skewness, and κ 4 (y) / κ 2 2 (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, for κ r (y) the r th 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 d j,k , which are updated by the observed coefficients d ″ j ,k to obtain posterior distributions d j,k | 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 (x i ).
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, where 0 ≤ γ j ≤ 1.0, δ(0) is a point mass at zero and d j,k are independent.The hyper-parameters are assumed to be of the form τ j 2 = C 1 2 −αj , γ j = min(1, C 2 2 −βj ) for non-negative constants C 1 and C 2 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 (x i ) into the prior.See reference [55].

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.

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] [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 metameter, and then apply Bayesian wavelet shrinkage using a universal threshold as developed in section 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

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

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

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 Discrete Wavelet Transforms -A Compendium of New Approaches and Recent Applications 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.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.

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: RANWD ≤ 0.5 and/or a WPB 90% ≤ 70%.
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.

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).[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).

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 RANWD≤ 0.50

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].diagnostics 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.
Density Estimation and Wavelet Thresholding via Bayesian Methods: A Wavelet Probability Band and Related… http://dx.doi.org/10.5772/52434 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).

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 Discrete Wavelet Transforms -A Compendium of New Approaches and Recent Applications 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 socalled 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].
onal scaling function, : mother wavelet : wavelet estimator for at level : Portion of synthesis due to scale , th 'detail' : ' mooth' of th order Kernel smoothing: Chebychev's inequality for the probability band [14] Relative average normalised density (RAND) Average normalised density (AND) Develop a physiologically representative model that incorporates endogenous agitation reduction (EAR).Use performance measures as follows: 1. RTD: relative total dose (RTD) expresses the total dose administered in the simulation as a percentage of the actual total recorded dose.2. Relative average normalised density (RAND) measures how probabilistically similar the model outputs are to the smoothed data, and hence the degree of comparability between the model and the empirical data.

Percentage time in band (TIB).
Rudge et al. [11] The agitation-sedation system model: Phamarcokinetic model adding patient agitation as a third state variable   Develop a nonparametric approach for assessing the validity of deterministic dynamics models against empirical data.Use performance measures as follows: 1. Kernel regression and density estimation to yield visual graphical display of data.2. Construct a probability band for the nonparametric regression curve and check whether the proposed model lies within the band.3. Average normalised density (AND) to measure how well the simulated values coincide with the maximum density at every time point and relative average normalised density (RAND).
Chase et al [14] The agitation-sedation Phamarcokinetic model adding patient agitation as a third state variable     Infinite Impluse Response (IIR) filter Proportional-Derivative (PD) control with agitation for infusion rate (U) Tracking Index (TI) Chebychev's inequality for probability band [14] Develop a mathematical model to capture the essential dynamics of the agitation-sedation system and test for statistical validity using the recorded infusion data for the 37 ICU patients.Use performance measures as follows: 1. Kernel smoothing using the uniform kernel.2. Tracking Index (TI).3. Moving blocks bootstrap to gain an understanding of the reliability of the TI for a given patient's infusion profile 4. 90% Probability Band -by definition the range within at least 90% of the time, the estimated mean value of the recorded infusion rate lies within the band.Chebychev's inequality for the probability band [14] Develop a nonparametric approach for assessing the validity of deterministic dynamics models against empirical data.Use performance measures as follows: 1. Kernel regression and density estimation to yield visual graphical display of data.2. Construct a probability band for the nonparametric regression curve and check whether the proposed model lies within the band.3. Average normalised density (AND) to measure how well the simulated values coincide with the maximum density at every time point and relative average normalised density (RAND).
Chase et al [14] The agitation-sedation Phamarcokinetic model adding patient agitation as a third state variable     Infinite Impluse Response (IIR) filter Proportional-Derivative (PD) control with agitation for infusion rate (U) Tracking Index (TI) Chebychev's inequality for probability band [14] Develop a mathematical model to capture the essential dynamics of the agitation-sedation system and test for statistical validity using the recorded infusion data for the 37 ICU patients.Use performance measures as follows: 1. Kernel smoothing using the uniform kernel.2. Tracking Index (TI).3. Moving blocks bootstrap to gain an understanding of the reliability of the TI for a given patient's infusion profile 4. 90% Probability Band -by definition the range within at least 90% of the time, the estimated mean value of the recorded infusion rate lies within the band.
C c , C p and C e are, respectively, the drug concentrations (mg L −1 ) 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 1 -K 3 are parameters relating to drug elimination and transport and K ij the transfer rate (L min −1 ) from compartment i to compartment j; K CL the drug clearance (L min −1 ); K T the effect, and w 1 and w 2 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 −1 ); A an agitation index; S the stimulus invoking agitation; K ij the transfer rate (L min −1 ) from compartment i to compartment j; K CL the drug clearance (L min −1 ); K T the effect time constant (min −1 ); P o and P s 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 w 1 and w 2 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].

Figure 1 .
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]).
Figures 2-3 are clearly more informative than the histogram

Figure 2 .
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 .
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 .Figure 5 .
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.

Figure 6 .
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 .
Figure 7. Box and whisker plot of the WTCI index for each of the 37 patients.
Uniform kernel with bandwidth h[65] Uniform kernel with bandwidth h[65]

Table 1 .
, 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.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

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.

Table 2 .
Alternative [29]let 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].

Table 4 .
Density Estimation and Wavelet Thresholding via Bayesian Methods: A Wavelet Probability Band and Related… http://dx.doi.org/10.5772/52434Comparison between the WPB, ANWD and RANWD and TIB, AND and RAND from Rudge's Physiological Model

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