Modelling Agitation-Sedation (A-S) in ICU: An Empirical Transition and Time to Event Analysis of Poor and Good Tracking between Nurses Scores and Automated A-S Measures

Abstract

Sedation in the intensive care unit (ICU) is challenging, as both over- and under-sedation are detrimental. Optimal sedation and analgesic strategies, are a challenge in ICU and nurses play a major role in assessing a patient’s agitation levels. Assessing the severity of agitation is a difficult clinical problem as variability related to drug metabolism for each patient. Multi-state models provide a framework for modelling complex event histories. Quantities of interest are mainly the transition probabilities e.g. between states, that can be estimated by the empirical transition matrix (ETM). Such multi-state models have had wide applications for modelling complex courses of a disease. In this chapter the ETM of multi-state and counting process (survival analytic) models which use the times for ICU patients to transition to varying states of violations (a violation being a carer’s agitation rating outside so-called wavelet-probability bands (WPB)) confirm the utility of defining so-called trackers and non-trackers according to WPB-based control limits and rules. ETM and multi-state modelling demonstrate that these control-limit scoring approaches are suitable for developing more advanced optimal infusion controllers and coding of nurses A-S scores. These offer significant clinical potential of improved agitation management and reduced length of stay in critical care.

1. Introduction

Pain management is increasingly recognised as a formal medical subspecialty worldwide [1]. Optimal sedation and analgesic strategies, combined with delirium management, are a challenge when caring for critically ill patients. Sedation in ICU aims to provide patient pain reflief, comfort and safety. For sedation monitoring, the most extensively used tools are RASS (Richmond Agitation and Sedation Scale) [2] and SAS (Sedation Agitation Scale) [3]. Despite extensive improvements in analgesia medication there are still barriers to nurses’ assessment, management, documentation, and reassessment of pain [4, 5, 6].

Pain is the most common reason that patients come to the emergency department. Emergency nurses have an indispensable role in the management of this pain [7, 8]. Sedation in the intensive care unit (ICU) is challenging, as both over- and under-sedation are detrimental. Current methods of assessment, such as the Richmond Agitation Sedation Scale (RASS), are measured intermittently and invariable depend on patients’ behavioural response to stimulation, as such may interrupt sleep and rest. A non-stimulating method for continuous sedation monitoring may be beneficial and allow more frequent assessment., noting that appropriate sedation cycling has to accommodate patients’ oscillations between states of agitation and over-sedation, which are detrimental to patient health and increases hospital length of stay [9, 10, 11, 12, 13, 14].

As such there also have been recent studies exploring the impact of augmenting sedation assessment with physiologic monitors [15] and studying the correlation between observational scales of sedation and bispectral index scores [16]. Recently the feasibility of continuous sedation monitoring of ICU patients using the NeuroSENSE was studied and suggested that such a non-stimulating method for continuous sedation monitoring may benefit patient care and allow increased A-S assessment [17]. The authors advocated use of incorporating some degree of automation into sedative drug administration, e.g. closed-loop control based on feedback from a processed EEG monitor, and various studies have suggested the limitations of RASS as a stand-alone measure of sedation levels, and pointed to benefit of adjunct continuous e.g., brain monitoring [17].

Earlier, Rudge, Chase, Shaw, Lee [12] discussed target controlled infusion (TCI) systems to deliver drugs to maintain target plasma concentrations, using a pharmacokinetic model, shown to be feasible when anaesthesia is given over short periods of reduced consciousness and well-known pharmacology is invoked. Infusion systems that regulate the infusion rate to maintain target agitation levels, to regulate the primary metric for longterm sedation, are one approach to improving care in the ICU. The data analysed in this chapter pertains to the scenario and data type studied earlier by [9, 10, 11, 12, 13, 14].

Assessing the severity of agitation is a challenging clinical problem as variability related to drug metabolism for each individual is often subjective. A multitude of previous studies suggest that the assessment accuracy of the sedation quality conducted by nurses tend to suffer from subjectivity and lead to sub-optimal sedation [14, 15, 18]. For example, [19] strongly recommend lighter than deeper levels of sedations. Moreover, [20, 21] argue that sedation should be reviewed and adjusted regularly. Whilst agitation management methods frequently rely on subjective agitation assessment [2, 3] the carers then select an appropriate infusion rate based upon their evaluation of these scales, experience, and intuition [21]. This approach usually leads to largely continuous infusions which lack a bolus-focused approach, commonly resulting in over or under-sedation. The work of [11, 12, 13] aimed to enhance feedback protocols for medical decision support systems and eventually automated sedation administration. A minimal differential equation model to predict or simulate each patient’s agitation-sedation status over time was presented in [12] for 37 ICU patients and was shown to capture patient A-S dynamics. The use of quantitative modelling to enhance understanding of the agitation-sedation (A-S) system and provision of an A-S simulation platform are one of the key tools in this area of patient critical care. A more refined A-S model, which utilised regression with an Epanechnikov kernel was formulated by [12]. A Bayesian approach using densities and wavelet shrinkage methods was later suggested by [9] to assess a previously derived deterministic, parametric A-S model [10, 11, 12, 13, 14], thus successfully challenging the practice of sedating ICU patients using continuous infusions. Wavelets approaches [9, 10] were shown to provide reliable diagnostics and visualisation tools to assess A-S models, giving alternative metrics of A-S control to assess validity of the earlier A-S deterministic models (Table 1 in [10]).

	V1	V2	V3	Total V’s	Time in ICU	WPB%
P18/Good	2	24	26	20	64	93.8%
P28/Poor	1	5	12	114	203	50.8%

Table 1.

Time to the patient-specific 1st violation V1, second violation V2 and third violation V3, total number of violations, total ICU time and WPB% values.

This suite of wavelet metrics based on the discrete wavelet transform (DWT) were able to establish the value of earlier deterministic agitation-sedation (A-S) models against empirical (recorded) dynamic A-S infusion profiles, providing robust performance metrics of A-S control and excellent tools, as based on the classification of patients into poor and good trackers based on Wavelet Probability Bands (WPBs). Importantly, the WPBs were shown as a useful patient-specific method by which to identify and detect regions in the patient’s A-S profile i.e., times whilst in ICU, where the simulated infusion rate performs poorly, thus providing visual and quantified ways to help improve and distil the deterministic A-S model and in practice be a guage to alert carers.

In this chapter Empirical Transition Matrix (ETM) approach of multi-state counting process (survival analytic) models of Allignol and coauthors [22, 23], aligned with the counting process/event history work in [24, 25, 26], which use the times patients transition to varying states of violations (a violation being an A-S measure outside the 90% WPB bands), confirm the utility of defining trackers and non-trackers according to these control limits and wavelet diagnostic rules of Kang et al., [9, 10]. In this chapter ETM and multi-state modelling are found to be valuable for developing advanced optimal infusion controllers and also to assist coding of nurses A-S scores, which potentially offer significant clinical potential of improved agitation management and reduced length of stay, as an augmented approach to also using RASS and SAS. Establishing patient-specific thresholds of poor A-S management and control has significant implications for the effective administration of sedatives, as improved management of A-S states will allow clinicians to improve the efficacy of care and reduce healthcare costs [27, 28, 29].

2. Data and methods

This chapter models the agitation-sedation profiles of Agitation and Sedation (A-S) profiles of 37 patients were collected at the Christchurch Hospital, Christchurch School of Medicine and Health Sciences, NZ. Two measures were recorded for each patient: (1) the nurses’ ratings of a patient’s agitation level and (2) an automated sedation dose (see Figure 1). Infusion data were recorded using an electronic drug infusion device for all admitted ICU patients during a nine-month observation period and required more than 24 hours of sedation. Infusion data containing less than 48 hours of continuous data, or data from patients whose sedation requirements were extreme, such as those with severe head injuries, were excluded [9, 10]. A total of 37 ICU patients met these requirements and were enrolled in the study. Classification of patients into poor and good trackers, as based on the Wavelet Probability Bands (WPB) are given in Table 2. The so-called good tracker delineates the scenario where the nurse’s rating scores remains within the (time based) 90% coverage of wavelet probability band (WPB) based on the simulated dose profiles [9, 10]. Poor tracking delineates the scenario where the nurses rating scores remain outside the (time based) 90% coverage of wavelet probability band (WPB) for a significant portion of time based on the simulated dose profiles [13, 14].

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 Chase et al. [14]).

WPB [9]	WCORR [10]	Chase et al. [14]	Rudge et al. [11]
2	2	—	—
4	4	—	—
—	—	6	—
7	7	7	7
9	9	9
10	10	—	10
11	11	—	11
—	—	12	—
—	—	—	13
—	—	17	—
21	21	21	—
22	22	—	22
27	27	27	27
28	28	—	28
—	29	—	29
32	32	—	—
33	33	—	33
34	34	34	—
—	35	—	35
Total: N1 = 13	Total: N2 = 15	Total: N3 = 8	Total: N4 = 10

Table 2.

Patient numbers of the poor trackers according to the criteria of 4 studies. Developed earlier in [11, 12, 13, 14]. Low WPB 90% indicates a poor tracker by Kang’s WPB diagnostics [9, 10].

By way of illustration we consider four patients from the pool of 37 patients.Tables 1 and 3 summarise each of these four patients’ WPB tracker status, time to first, second and third violation outside the WPB bands [9], their total number of violations over ICU stay and patient’s time in ICU, along with their specific WPB% value. Display of their line profiles of nurses’ rating of A-S in relation to drug infusion dose over time, for each of the 4 patients (P8, P27, P18, P28) are given in Figures 2–4.

	V1	V2	V3	Total V’s	Time in ICU	WPB%
P8/Good	1	2	3	46	128	87.5%
P27/Poor	1	4	5	89	225	43.7%

Table 3.

Time to the patient-specific 1st violation V1, second violation V2 and third violation V3, total number of violations, total ICU time and WPB% values.

Figure 2.

Line plot of nurses’ rating of patient agitation and the automated sedation dose for patient 8 (good tracker).

Figure 3.

Line plot of nurses’ rating of patient agitation and the automated sedation dose for patient 27 (good tracker).

Figure 4.

Line plot and WPB% band for patient 18 (LHS) and 28 (RHS, poor tracker).

The first patient (patient 8) in Table 3 is a good WPB tracker and the second a poor WPB tracker (patient 27), studied in depth in [28], for which upper tail thresholds of the nurses’ scores using copulas were established. We also refer the reader also to Hudson & Tursunalieva’s chapter in this book entitled “Copula thresholds and modelling Agitation-Sedation (A-S) in ICU: analysis of nurses scores of A-S and automated drug infusions by protocol” [27]. The corresponding WPB% values for patient 8 and patient 27 are 87.5% and 43.7%, respectively (Table 3). Overall, the minimum, median and maximum WPB% values for the 24 good trackers is (58.8%, 87.5%, 96.9%) and (47.3%, 64.8%, 77.3%) for the 13 poor trackers (Table 2). Noteworthy also is that the A-S time series of these two patients examined (P8 and P27) were of disparate lengths - patient eight had 10,561 time points and patient 27, 13,441 time points. The full 37 patients studied had a range of [3001–25,261] time points.

Patient 18 (good tracker) with a WPB% of 93.8% and patient 28 (poor tracker) with WPB% of 50.8% (Table 1) were studied in detail in [28], for which both upper and lower tails/thresholds of over or under-estimation of agitation levels by the nurses rating were established using copula dependence analytics [29], refer also to [27].

Patients vary according to their length of stay in ICU and consequently differ in their opportunity for violations to occur. The good trackers generally have shorter ICU time and thus less chance to exibit an increased total number of violations. An indication of how the strata (good versus poor tracker), the patient’s total number of violations and a patient’s time in ICU interact, can be visualied in Figure 5. The total number of WPB based violations is clearly greater for the poor trackers than for the good trackers, and it is the poor trackers that tend to have longer ICU times. Also from the scatterplot in Figure 5 there seems to be three approximate categories of patient ICU time: 50–64, 113–128 and 205–256. The majority of patients (28 (76%)) have ≤40 violations (RHS of Figure 5), 19 (51%) patients have an ICU time of ≤64 (Table 4).

Figure 5.

Total number of violations by WPB tracking status [9] and ICU time.

	ICU category
WPB tracker	0	1	2
Good	15	7	2
Poor	4	4	5

Table 4.

Patient tracker status by ICU time: 0 = 50–64, 1 = 113–128, 2 = 205–256.

Accordingly, for ICU time categorised and coded as: 0 = [50,64], 1 = [113,128], and 2 = [205,256], the total violations profile according to tracking status, displayed in Figure 6, shows that the total number of violations is significantly higher for the poor trackers, particularly when ICU time > 205. Noteworthy, is that the majority of patients 28 (76%) have ≤40 violations (RHS of Figure 6), whereas 19 (51%) patients have an ICU time ≤ 64 (Table 4). Figure 7 displays the histogram of the number of violations where a violation is defined as a nurse’s A-S rating outside the patient’s WPB control band. We note that the majority of the time, in excess of >75% of the 370 violation counts are below a count of five violations (Figure 7).

Figure 6.

Total number of violations by 3 levels of ICU time (LHS) and boxplot of poor good tracker time to 3^rd violation by ICU time: 0 = [50,64], 1 = [113,128], 2 = [205,256].

Figure 7.

Violation counts in bins across all patients.

For the state-space analysis described in Section 3 each patient’s total ICU time is broken into 10 bins, where each bin represents 10% of the patients’ total time in ICU; i.e., Bin 1: 0–10%, Bin 2: 11–20% etc. For the 37 patients, we thus have 370 bins, i.e. 370 counts of violations. The 10% interval approach is used due to the large variation in time in ICU between the WPB-based good versus poor strata [9] - noting that some poor trackers have times up to 256, whereas good trackers are mostly limited to 64–128. Given these bins, patients’ A-S states can then be defined in terms of the total number of violations or jumps outside the WPB bands that occur during each 10% interval of a patient’s total ICU time. The randon, outcome event of A-S status is then the number of violations that over time.

3. Empirical transition matrix state-space & commenges’ test approach

3.1 Mathematical formulation

Multi-state models are known to provide a relevant framework for modelling complex event histories. Quantities of interest are mainly the transition probabilities that can be estimated by the empirical transition matrix, that is also referred to as the Aalen-Johansen estimator [30, 31]. Such multi-state models have had a wide range of applications for modelling complex courses of a disease over the course of time and across applications in medical research (Beyersmann et al. [32], Munoz-Price et al. [33], Andersen & Keiding [34]). We now utilise the Empirical Transition Matrix (etm) approach to model multi-state models of [22] and derive inference tests for such models using the approach of Commenges [25, 35, 36] with a particular focus on Commenges’ test derived in earlier work [25].

Define patient states as follows, any state can transition into any other state (Figure 8).

Figure 8.

State-space system.

• State 1: 0 or 1 violations

• State 2: 2 or 3 violations

• State 3: > 3 violations

A number of different approaches (etm on ICU time, etm on bin time, and log-rank type tests as in Commenges [25], will be used to investigate the difference between good and poor trackers in terms of a devised a 3 state transition formulation as defined below. Differences between transition probabilities between states for the good and poor trackers will be evaluated using Commenges’ [25] chi square test.

The mathematics is well described in the work of Allignol [22], adapted to more complex scenarios in [23]; and in Commenges’ approach [25, 35, 36]. The mathematical formulation of the ETM state-space approach and Commenges’ test are given for general frameworks as follows.

Consider a stochastic process $\left(X_t\right)$ with finite state space $S=\left\{1,\dots ,K\right\}$ where sample paths are right-continuous, and the stochastic process is assumed to be time-inhomogeneous Markov. The transition hazard from state $i$ to state $j$, $i\ne j$ is defined as

$${\alpha }_{\mathit{ij}}\left(t\right)\mathit{dt}=P\left(X_{t+\mathit{dt}}=j|X_t=i\right).$$

The cumulative hazard transitions are defined as,

$$\begin{array}{c}{A}_{\mathit{ij}}\left(t\right)={\int }_0^t{\alpha }_{\mathit{ij}}\left(t\right)\mathit{dt}\\ A_{\mathit{ii}}\left(t\right)=-\sum _{j\ne i}{A}_{\mathit{ij}}\left(t\right).\end{array}$$

Define

$$P_{\mathit{ij}}\left(s,t\right)=P\left(X_t=j|X_s=i\right),\text{for}i,j\in S,s\le t$$

as the probability that an individual who is in state $i$ at time $s$ is in state $j$ at time $t$.

The (K + 1)x(K + 1) probability transition matrix with elements $P_{\mathit{ij}}\left(s,t\right)$ can then be obtained from the transition hazards through product integration. Let $N_{\mathit{ij}}\left(t\right)$ be the number of observed direct transitions from state I to state j up to time t and let $Y_i\left(t\right)$ be the number of individuals under observation in state I just before time t. The Nelson-Aalen estimator is used to estimate the non-diagonal elements of the matrix of cumulative hazards as follows,

$${\widehat{A}}_{\mathit{ij}}\left(t\right)={\int }_0^t\frac{d{N}_{\mathit{ij}}\left(u\right)}{Y_i\left(u\right)},i\ne j$$

and the diagonal elements ${\widehat{A}}_{\mathit{ii}}\left(t\right)$ are obtained as above. The product integration relationship below leads to an estimate of the probability transition matrix as follows,

$$\widehat{\mathit{P}}\left(s,t\right)=\prod _{s<t_k\le t}\left(\mathit{I}+\mathrm{\Delta }\widehat{\mathit{A}}\left(t_k\right)\right),$$

where the product is taken over all possible transition times in time interval (s,t].

An estimator for the covariance of the empirical transition matrix is given by,

$$\widehat{cov}\left(\widehat{\mathit{P}}\left(s,t\right)\right)={\int }_s^t\widehat{\mathit{P}}{\left(u,t\right)}^T⨂\widehat{\mathit{P}}\left(s,u-\right)\widehat{cov}\left(d\widehat{\mathit{A}}\left(u\right)\right)\widehat{\mathit{P}}\left(u,t\right)⨂\widehat{\mathit{P}}\left(s,u-\right){,}^T$$

which is a (K + 1)² x (K + 1)² matrix, with number or rows and columns equal (K + 1)².

3.2 Commenges’ test formulation

We now utilise the framework of the generalised Cochran–Mantel–Haenszel (CMH) test for (I x J x K) tables. The CMH test is based on the hypergeometric distribution. The CMH and the test of Commenges’ for the specific case here, where we have three states (of violations) and two strata (good versus poor trackers) is now described.

In our application then we have 2 × 3 × 3 tables. For each k = 1, 2, 3 we have a 2 × 3 table, where the elements are counts n_ijk (k denotes the departure state j, so the rows, I, are the WPB strata, the columns, J, are the entry states I and the slices, K, are the departure states j), as tabulated below.

	I = 0	I = 1	I = 2	Row total
Good	n11k	n12k	n13k	n1 + k
Poor	n21k	n22k	n23k	n2 + k
Total	n+1k	n+2k	n+3k	n++k

Assume that the row and column marginals are fixed. This implies that there are (I-1) by (J-1) values that are free to vary. For the cell n_ijk the expected value is then given by n_i + k x n_+jk. We are able to express all cell counts by the vector n_k and express all expected cell counts by u_k. The covariance matrix, denoted by V_k then has elements,

$$\mathit{cov}\left(n_{\mathit{ijk}},n_{i^{\prime }{j}^{\prime }k}\right)=\frac{n_{i+k}\left({\delta }_{i{i}^{\prime }}{n}_{++k}-n_{i^{\prime }+k}\right)n_{+\mathit{jk}}\left({\delta }_{j{j}^{\prime }}{n}_{++k}-n_{+j^{\prime }k}\right)}{n_{++k}^2\left(n_{++k}-1\right)},$$

where ${\delta }_{\mathit{ab}}=1$ when a = b, and 0 otherwise. Assuming rows and columns are unordered we sum over the K strata to obtain,

$$\mathit{n}=\sum {\mathit{n}}_{\mathit{k}},\mathit{u}=\sum {\mathit{u}}_{\mathit{k}},V=\sum {\mathit{V}}_{\mathit{k}}.$$

The generalised CMH statistic is then given by,

$$X^2={\left(\mathit{n}\mathit{-}\mathit{u}\right)}^{\prime {\mathit{V}}^{-1}}\left(\mathit{n}\mathit{-}\mathit{u}\right),$$

which follows a chi square distribution with (I-1) by (J-1) degrees of freedom. This test is implemented in R via the mantelhaen.test function. Commenges [25] adapts this concept, with a test which differs to the generalised CMH test in that it does not sum over the K strata, before calculating the relevant chi squared test statistic.

Commenges’ test [25] is as follows,

$$X_k^2={\left({\mathit{n}}_{\mathit{k}}-{\mathit{u}}_{\mathit{k}}\right)}^{\prime {\mathit{V}}^{-1}}\left({\mathit{n}}_{\mathit{k}}-{\mathit{u}}_{\mathit{k}}\right),$$

where $X_k^2$ is chi-square with (I-1) by (J-1) degrees of freedom.

The required total chi squared statistic is then simply obtained by taking $X^2=\sum X_k^2$ which is itself disctributed as a chi square distribution with K(I-1) by (J-1) degrees of freedom.

3.3 Results of the ETM analysis and Commenges’ test on transition states

Conditionally on the number of patients in each state at each step we have 3 x 9 = 27 independent contingency tables (i.e., number of departure states j by the number of time points k–1, recall we have 10, 10% bins for a patient’s time in ICU) and each of these tables has dimension 2 × 3 (good/poor tracker by the number of states i).

The corresponding three specific strata tables are given in Tables 5–7. For example, for departure state j = 0 and time point k = 2 we have a two by three contingency table with an overall total of 7^ϕ violations (labelled ^ϕ in Table 5); the latter informs that, at time k = 2 there are 5^ϕ good WPB based trackers that depart from state j = 0 and enter state 0. Similarly there are two (2^ϕ) WPB based poor trackers that depart state j = 0 and enter state 1 (Table 5).

	Strata	I = 0	I = 1	I = 2	Total
k = 2	Good	5^ϕ	0	0	5
	Poor	0	2^ϕ	0	2
	Total	5	2	0	7^ϕ
k = 3	Good	10	5	1	16
	Poor	0	1	0	1
	Total	10	6	1	17
k = 4	Good	10	1	0	11
	Poor	2	0	0	2
	Total	12	1	0	13
k = 5	Good	9	2	2	13
	Poor	0	3	1	4
	Total	9	5	3	17
k = 6	Good	9	3	2	14
	Poor	0	0	0	0
	Total	9	3	2	14
k = 7	Good	9	2	2	13
	Poor	1	1	0	2
	Total	10	3	2	15
k = 8	Good	5	4	2	11
	Poor	2	0	2	4
	Total	7	4	4	15
k = 9	Good	5	2	0	7
	Poor	3	0	0	3
	Total	8	2	0	10
k = 10	Good	9	1	1	11
	Poor	2	1	2	5
	Total	11	2	3	16

Table 5.

Departure state j = 0: WPB strata [9].

	Strata	I = 0	I = 1	I = 2	Total
k = 2	Good	6	1	1	8
	Poor	1	1	0	2
	Total	7	2	1	10
k = 3	Good	1	2	0	3
	Poor	2	2	2	6
	Total	3	4	2	9
k = 4	Good	2	2	3	7
	Poor	2	2	0	4
	Total	4	4	3	11
k = 5	Good	3	1	1	5
	Poor	0	1	2	3
	Total	3	2	3	8
k = 6	Good	3	1	0	4
	Poor	1	3	1	5
	Total	4	4	1	9
k = 7	Good	2	3	1	6
	Poor	3	2	0	5
	Total	5	5	1	11
k = 8	Good	2	2	3	7
	Poor	1	1	2	4
	Total	3	3	5	11
k = 9	Good	4	3	2	9
	Poor	0	1	0	1
	Total	4	4	2	10
k = 10	Good	4	1	2	7
	Poor	1	1	0	2
	Total	5	2	2	9

Table 6.

Departure state j = 1: WPB strata [9].

	Strata	I = 0	I = 1	I = 2	Total
k = 2	Good	5	2	4	11
	Poor	0	3	6	9
	Total	5	5	10	20
k = 3	Good	0	0	5	5
	Poor	0	1	5	6
	Total	0	1	10	11
k = 4	Good	1	2	3	6
	Poor	0	1	6	7
	Total	1	3	9	13
k = 5	Good	2	1	3	6
	Poor	0	1	5	6
	Total	2	2	8	12
k = 6	Good	1	2	3	6
	Poor	1	2	5	8
	Total	2	4	8	14
k = 7	Good	0	2	3	5
	Poor	0	1	5	6
	Total	0	3	8	11
k = 8	Good	0	3	3	6
	Poor	0	0	5	5
	Total	0	3	8	11
k = 9	Good	2	2	4	8
	Poor	2	1	6	9
	Total	4	3	10	17
k = 10	Good	2	1	3	6
	Poor	1	0	5	6
	Total	3	1	8	12

Table 7.

Departure state j = 2: WPB strata [9].

Three 2 x 3 contingency tables (one for each departure state j) are thus created.

Estimated transition probabilities for the 3 state process are then plotted using the ‘xyplot’ function from the lattice package in R. In the resultant plots (Figures 9–11), the vertical y-axis represents the transition probability value, which is represented by the solid line in each plot region. The numbers in the coloured bar above each plot defines the transition (e.g., 1 2 means transition probability from state 1 to state 2). The dotted lines around the solid line represent the confidence bands based on the covariance as calculated by the etm function. The horizontal x-axis shows the the time i.e., 10% bins (i.e. 2–10, because no transitions occur at time one being the initial state). For each tracker status and possible piairs of state transitions there are three plots, given in the following order, good trackers, poor trackers. Figure 11 displays the probability of being in each of the 3 states (0, 1, 2) given the initial state is state 0.

Figure 9.

Transition probability profiles for WPB good trackers.

Figure 10.

Transition probability profiles for WPB poor trackers.

Figure 11.

Probability of being in each state as time progresses given start state 0. Top is good trackers, bottom is poor trackers: WPB based.

Our procedure results in three 2 x 3 contingency tables (one for each departure state j), see Tables 8–10. The chi squared statistic as derived in [25] can now be calculated in that for the Commenges test the same chi squared calculation is made for each state specific table separately (i.e., without summing over k). The results in this case are χ²(1) = 6.046, χ² (2) = 2.269 and χ²(3) = 9.280. Each of these follows a chi square distribution with 2 degrees of freedom with associated p-values of 0.049, 0.322 and 0.010 (Table 11). Kang WPB (2013) [9].

Tracker strata	I = 0	I = 1	I = 2	Total
Good	71	20	10	101
Poor	10	8	5	23
Total	81	28	15	124

Table 8.

Departure state j = 0 summed over all time points k, k = 2, …,10.

Tracker strata	I = 0	I = 1	I = 2	Total
Good	27	16	13	56
Poor	11	14	7	32
Total	38	30	20	88

Table 9.

Departure state j = 1 summed over all time points k, k = 2, …,10.

Tracker strata	I = 0	I = 1	I = 2	Total
Good	13	15	31	59
Poor	4	10	48	62
Total	17	25	79	121

Table 10.

Departure state j = 2 summed over all time points k, k = 2, …,10.

State	χ^2	p-value
J = 0	6.046	0.049*
J = 1	2.269	0.322
J = 2	9.280	0.010**
Total	17.59	0.007***

Table 11.

Computation of Commenges’ test for the WPB strata [9].

Summing these threeχ ² (j), j = 1,2,3 statistics gives a value of χ^{2 =}17.59 with 6 degrees of freedom and an associated p-value of 0.007 (Table 11). The underlying null hypothesis is that the two nominal variables (strata: good or poor tracker and entry state: 0, 1, or 2) are conditionally independent in each stratum (departure state j; 0, 1 or 2), assuming no three-way interaction. The low p-value of 0.007 suggests that this hypothesis be rejected, i.e., the two variables are not conditionally independent. Thus the Commenges test shows that there is a statistically significant difference between the good versus poor tracker WPB strata, and that this difference is mainly due to transitions out of states 0 and 2, which agrees with the trends based on a graphical inspection of Figures 9–11.

The same procedure and related Commenges’ test is then applied to each of the 3 remaining good/poor tracker definitions of Kang [10, 11, 14] for the three-state context studied in this chapter. These results are reported in Table 12.

State	χ^2	p-value
J = 0	1.724	0.422
J = 1	0.911	0.634
J = 2	7.123	0.028**
Total	9.758	0.135

Table 12.

Computation of Commenges’ test for the remaining A-S studies.

Kang et al., WCORR [10].

Chase et al. [14].

State	χ^2	p-value
J = 0	5.669	0.059
J = 1	6.406	0.041*
J = 2	3.097	0.213
Total	15.172	0.019**

Rudge et al. [11].

State	χ^2	p-value
J = 0	0.367	0.832
J = 1	2.951	0.229
J = 2	5.466	0.065*
Total	8.784	0.186

Transition probability profiles of being in each state as time progresses, given start state 0, for the remaining 3 studies of [10, 14, 11] are given in Figures 12–14. In summary, Figures 9–14 illustrate the trend that good trackers tend to have higher probability of transitioning into state 0 than poor trackers, and the good trackers tend to have lower probability of transitioning into state two than poor trackers, where state two indicates that more violations (>3 violations) are occurring, and state 0 indicates few violations are occurring.

Figure 12.

Probability of being in each state as time progresses given start state 0. Top is good trackers; bottom is poor trackers: WCORR of [10].

Figure 13.

Probability of being in each state as time progresses given start state 0. Top is good trackers; bottom is poor trackers according to Chase et al. [14].

Figure 14.

Probability of being in each state as time progresses given start state 0. Top is good trackers; bottom is poor trackers according to Rudge et al. [11].

Notably also, the probability of transitioning into state 2 overall appears to increase as ICU time increases. This is most likely because poor tracking patients tend to have longer ICU times, and so, as time goes on, it is only poor trackers transitions that are being estimated. By categorising patients according to total ICU time (≤64, >64) as discussed earlier (Figures 5 and 6, Table 4) some of this could be accounted for. The results obtained are still consistent, as shown in the etm profiles using ICU time (≤64, >64) in Figures 15 and 16, respectively. The corresponding ETM probabilities are determined according to etm in R [21] and associated state and strata specific plots given in Figures 15 and 16.

Figure 15.

Transition probability profiles for patients with ICU time ≤ 64. Top panel are the WPB good trackers, and bottom panel the poor trackers.

Figure 16.

Transition probability profiles for patients with ICU time > 64. Top panel are the WPB good trackers, and bottom panel the poor trackers.

3.4 Conclusion regarding the ETM based analysis

The different approaches in Section 3 led to the sasimilar conclusions that there is a difference in the way good trackers and poor trackers transition between states. Most of this difference occurs in states 0 and states 2, as defined. Good trackers tend to have higher probability of transitioning into state 0 than poor trackers, and good trackers tend to have lower probability of transitioning into state 2 than poor trackers, noting that state 2 indicates more violations are occurring, and state 0 indicates fewer violations. The probability of transitioning into state 2 overall appears to increase as ICU time increases. This is most likely due to the fact that poor tracking patients have longer ICU times, and so, as time goes on, it is only the poor trackers’ transitions that are being estimated.

By categorising patients according to their total ICU time (≤64, >64) similar trends were found. The Commenges’ test established a statistically significant difference between the two tracking strata (p = 0.007), and that this difference was mainly due to transitions out of states 0 and 2. For the tracking metric of Chase et al. [14], the Commenges test demonstrated a statistically significant difference between the two good versus poor strata (p = 0.019), with this difference mainly due to transitions out of states 0 and 1. Overall, the WCORR [10] and Rudge [11] classifications of tracking/strata, the transision probability profiles for the 3 state process, good and poor trackers are not significantly different, but exhibited some difference mainly due to transitions out of state 2.

4. Time-to-response bias as a counting process mult state model

4.1 Mathematical formulation

In this section we analyse the WPB violation data and investigate the times to the third violation given times to second violation for both poor trackers and good trackers where tracking status is defined by WPB diagnostics. The process can be thought to have three states. State 1 corresponds to less than two violations, state 2 means two violations and if a patient is in state 3 then three violations have occurred. This is a sequential three state process shown schematically in Figure 17. Events of interest are a transition from state 2 to state 3, i.e. the occurrence of a third violation.

Figure 17.

States of a patient’s agitation (violations) defined by certain levels of violations or jumps outside of the patient’s WPB bands - a 3 state process.

Time one is the patient’s entry time into state 2 and time two is the patients exit time from state 2 (i.e., entry time to state 3, so-called ‘death’ state). Time one is the patient’s entry time into state 2 and time two is the patient’s exit time from state 2 (i.e., entry time to state 3, so-called ‘death’ state).

In the case of multiple events of interest, the process can be treated as a Markov chain. Let N_ij (t) be the process counting the number of observed transitions from state i to state j in the interval [0, t]. The transition intensity from state i to state j at time t is then λ_ij (t) and gives the instantaneous risk of transition from state i to j.

N_ij (t) has intensity process of the form λ_ij(t)N_i(t) where Y_i(t) is the number of individuals in state i just before time t. This is the setup of Simon and Makuch [37] who considered 4 states and two transitions of interest.

The concept of time in our ICU application represents time on-study (i.e. time at ICU) rather than calendar time. In the case of our 3 state process (Figure 8) the hazard functions of the two transitions of interest are λ₁₃(t) and λ₂₃(t) and the number of individuals in the states just before t are N₁(t) and N₂ (t), respectively. A chi squared test is conducted to test for independence between response and non-response as in the development formulated in [37].

This same test can be conducted to assess the association between strata and hazard rate. If λ₂₃(0)(t) is the hazard rate (from state 2 to 3) for the good trackers and λ₂₃(1) is the hazard rate (from state 2 to 3) for the poor tackers. Since the focus here is to test the effect of response (prior 2nd violation) on the hazard function, the null hypothesis of interest is H₀: λ₂₃(0) (t) = λ₂₃(1). The hypothesis is tested via a log-rank type test following [37] which tests for the time-to-response bias. Now the equivalent of Table 2 in Simon and Makuch [37] can be constructed for both the good trackers and the poor trackers. Let N₁(t) be the number of patients in state 1 at time t, and let N₂ (t) denote the number of patients in state 2 at time t in Table 13.

Table 13 presents the WPB data as state-specific patient counts for each event time t. Events are a transition from state 2 to state 3, i.e. the occurrence of the event of interest i.e. a third violation. Time represents time to third violation (so-called end-state/death in terms of a counting process). Note that d_ij (t) are the so-called end-state “deaths” i.e., third violations. The hazard function for transfers between states i and j at time t is denoted by λ_{i j} (t) and here time represents time on study at ICU. Also T_ij denotes the set of times at which a transition from state i to state j occurs and N_i(t) is the number of patients in state i just before time t, or in other words, N_i(t) is the number of patients at risk of a transfer out of state i at time t. (Note our Ni (t) is equivalent to Y_i(t) in the Aalen’s notation). The symbol λ_ij (t) denotes the intensity, or hazard function, for a transfer from state i to state j at time t.

Time	N1	N2	Events
3	12	12	7
4	7	10	5
6	6	6	2
9	5	5	1
12	4	5	1
18	2	6	1
19	2	5	1
26	0	6	1
27	0	5	1
28	0	4	1
33	0	3	1
36	0	2	1
43	0	1	1
(a) Good trackers.

Time	N1	N2	Events
3	9	4	3
4	8	2	1
5	5	4	2
6	3	4	2
8	1	4	2
9	1	2	1
12	1	1	1
23	0	1	1
(b) Poor trackers

Table 13.

Simon and Makuch’s [37] represeantation and formulation of the WPB data.

Mathematically the cumulative hazard function is conventionally estimated instead of the hazard function λ(t), as the latter is difficult to estimate. The cumulative hazard function and survival function is then given as,

$$\widehat{A_{\mathit{ij}}}\left(x,t\right)=\sum -log\left(1-d_{\mathit{ij}}\left(u\right)/N_i\left(u\right)\right)$$

$$\widehat{S_{\mathit{ij}}}\left(x,t\right)=\prod \left(1-d_{\mathit{ij}}\left(u\right)/N_i\left(u\right)\right)$$

$$=exp\left(-\widehat{A_{\mathit{ij}}}\left(x,t\right)\right)$$

Note that d_ij(u) above are the so-called end-state “deaths” i.e., third violations, the number of transitions from state i to state j in time interval [x, t]. The estimated survival and cumulative hazard curves are shown as in Figure 18.

Figure 18.

Survival function of time to 3rd violation given the 2nd violation for good tracker and poor (non-)trackers (LHS) and cumulative hazard functions (RHS).

Survival curves and cumulative hazard functions were calculated according to Simon and Makuch’s method [37]. In essence, this counting process formulation keeps track of the number of patients in state 1 and 2 and event times (i.e., transitions into state 3). The survival package is used for estimation, where two times are used. Time one is the patients entry time into state 2 and time two is the patients exit time from state 2 (i.e. entry time to state 3, ‘death’).

The log rank test for H₀: λ₂₃(0) (t) = λ₂₃(1), based on the counting process which utilises the number of individuals in the states just before t, these are, N₁(t) and N₂ (t), was performed. Accordingly, it is shown that the good tracker and poor tracker hazard rates/ (survival curves) time to the 3rd violation, given a 2nd violation has occurred, are statistically significantly different (p-value = 0.044), see left hand side of Figure 18. Notably, the hazard rate for the poor trackers is 2.1 times that of good trackers, 95% confidence interval (CI) [1.01, 4.38].

Further interpretation of the hazard function can be made by assessing the slope of the cumulative hazard function. Figure 18 (RHS), shows that the cumulative hazard increases faster for the poor trackers than the good trackers indicated by a much steeper slope. This suggests it takes less time for the poor trackers to reach their third violation than for the good trackers, this is also confirmed by the 95% confidence bands for the survival curves shown in Figure 19 for the good tracker and poor trackers. Note that the interpretation of Kaplan–Meier curves here is not as straight-forward as for conventional survival analysis. In our ICU A-S process formulation the curves do not correspond to fixed cohorts, as patients can contribute to different states/curves at different times (Table 13). Thereby the curves may be considered to represent hypothetical cohorts whose values remain constant after follow-up [38, 39].

Figure 19.

Survival curves (95% CIs) for good (left) and poor (right) (non-)trackers.

A Cox proportional hazards model (CPHM) was then fitted with tracking status and a patient’s number of violations as covariates. The general CPHM hazard function is,

$$\lambda \left(\mathrm{t}|X\right)={\lambda }_0\exp \left({\beta }_1,X_1,+,\dots ,+,{\beta }_p,X_p\right)$$E13

In our application we model two covariates: X₁ (0 for a good tracker, 1 for a poor tracker), and X₂ being the patient’s total number of violations in the CPHM. The log rank test associated with the CPHM confirmed that the good tracker and poor tracker hazard rates, and the survival curves were significantly different (p-value = 0.0496), with the hazard rate for the poor trackers being 2.1 times that of good trackers, with a 95% confidence interval of [1.01, 4.38]. The associated hazard rate for poor trackers is shown to be 1.87 times that of good trackers, with a 95% confidence interval [0.75, 4.70]. By inclusion of the total number of second time violations the effect of tracking status has only reduced slightly, and it remains significant (1.87 versus 2.1).

4.2 Asseement of times to different violation counts and patient’s last jump

Log-rank tests were likewise conducted to assess times to different violation counts. Let VX denote the violation times for the Xth violation and the DX’s the associated event indicators (0 censored, 1 event of interest). Log rank tests for the two WPB tracking strata for selected violation times (X = 5, 10, 15, 20, 25, 30) showed significant differences between good and poor WPB trackers regarding the time to the patient’s time to 10th violation (p = 0.027), their 15th (p = 0.025) and their 25th violation (p = 0.011). Likewise, significance at the 10% level was demonstrated for times to the patient’s 5th, 20th and 30th violation (non-violatory lifetimes). All survival curves (not shown here) are significantly different or are close to being significant at the 5% level of significance. This confirms that the difference in time to violations between the good and poor trackers are consistently different, for these varying number of violations (VX, for X = 5, 10, 15, 20, 25, 30).

The time to a patient’s last violation event was also investigated using log-rank tests and Kaplan–Meier curves. We examined nine levels of the effect of the following covariate, which categorises the counts the patient levels of violations as follows: 0–5 violations, 5–10, 10–15, 15–20, 20–25, 25–30, 30–40, 40–50 and >50 violations. A histogram of the time to a patient’s last violation with boxplots of the times for each of these nine levels of categorisation shown is given in Figure 20 (the number above the boxplots gives the number of patients in each of the nine categories).

Figure 20.

Histogram and boxplots of time to last violation. The numbers above the boxplot (RHS) specify the number of patients in each violation level.

Using the patient’s time to their final violation/jump, as the event of interest, and implementing log-rank based tests using this covariate adjustment, the log rank test demonstrated a statistically significant difference between the survival curves of time to last violation (p-value <0.000001) across the above nine different total number of violation levels, {0–5, 5–10, 10–15, 15–20, 20–25, 25–30, 30–40, 40–50, >50 violations}.

Notably, the levels that most contribute to the difference between trackers and non-trackers are those patients who have a total number of violations between 10 and 15, between 20 and 25 and >50, in that order. A log rank test on time to last violation (time of last jump outside the WPB bands) as the outcome of interest by tracking status, also establishes that there is a difference between the two survival curves (p-value = 0.045) (Figure 21). Clearly the WPB-based poor trackers tend to take longer to reach their last violation than good trackers. Corresponding Kaplan–Meier estimated curves are given in Figure 21. We note that up to ICU time 64 (≤64), 40% of the good versus 70% of the poor trackers are still violating, whereas after time point, 130, the corresponding percentages violating are 15% versus 40%, of the good versus poor trackers (Figure 21).

Figure 21.

Estimated survival curves for time to last violation by tracking status.

5. Conclusion

A log-rank test from the counting process formulation [37, 38, 39] established a significant difference between the hazard curves of the WPB-based good and poor trackers (p-value = 0.044). Similarly log-rank tests performed for a variety of violation numbers to test for differences between good and poor trackers times to their 5th, 10th, 15th, 20th, 25th and 30th violation, showed evidence of a significant difference between good and poor trackers for a selection of these violation times (namely patient’s time to their 10th, 15th and 25th violation). In regard to analysing the patients time to their last recorded violation, log-rank tests and Kaplan–Meier curves showed that poor trackers tend to have a higher probability of still violating as time progresses in ICU compared to good trackers (p-value = 0.045).