RNN-BOF: A Multivariate Global Recurrent Neural Network for Binary Outcome Forecasting of Inpatient Aggression
^††thanks: This research was supported by the Australian Research Council under grant DE190100045.

We have identified two studies that directly employed psychometric instruments and clinical covariates alongside machine learning techniques. Galatzer-Levy et al. [3] tested the efficacy of various models including Support Vector Machines (SVM), Random Forests (RF) and kernel ridge regressions trained to predict the presence of a traumatic stress disorder diagnosis within 15 months of a traumatic event. The models were trained on patient data collected at the time of the event, general demographic data, physiological responses, and an assessment on the Acute Stress Disorder scale psychometric instrument. Their approach was able to outperform the benchmark predictive power of the scale itself. A similar work [4] assessed high-risk psychiatric inpatients using the Suicide Crisis Inventory instrument which is composed of 49 five-point ordinal scales. They trained models including RFs and gradient boosted decision trees (GBDT) to predict the occurrence of suicide attempts over a one month period although they made no comparison to the baseline psychometric evaluation.

Statistical local forecasting methods have consistently out-performed neural networks and other machine learning techniques in forecasting competitions for decades [11]. This changed only relatively recently in the M4 forecasting competition, which was won by a hybrid model that employed both statistical and neural network forecasting, training models globally across all time series [12]. Following this success and a subsequent rise in popularity, Montero-Manso and Hyndman [6] established theoretical bounds of effectiveness for global cross-series approaches as compared to local (per series) methods, suggested that the multivariable equivalent would inherit the beneficial properties of global cross-series methods, and offered empirical evidence for the lack of necessity of homogeneity across series. Hewamalage et al. [13, 14] supported these conclusions and identified RNNs with LSTM cells (RNN-LSTM) and Light Gradient Boosted Machines (LGBM) as strong performers on both real world and simulated datasets.

However, to the best of our knowledge there is limited literature on multivariable time series forecasting where models are trained globally across many homogeneous multivariable datasets. The two multivariable global forecasters we have identified are DeepGLO [7] and TEDGE [8]. Both of these models learn from a set of time series, where each series has covariate time series and static categorical features. Both use initial matrix factorisation to extract global trends; DeepGLO then uses a temporal convolution network to learn local patterns, while TEDGE makes use of a multivariable RNN-LSTM architecture. Either of these forecasters could suit our needs if not for the probabilistic or scored binary forecast requirement.

The majority of models for binary time series forecasting are from the auto-regressive family, with neural network approaches being scarce [15]. gbARMA [16], GLARMA [17] and PC-ARIMA [18] are capable forecasters given a multivariate forecasting problem with a binary target series. However, they are all locally trained models [7], meaning for a patient without any instances of the adverse event, the model would have no historical reference and therefore could not predict future instances. Bautu et al. [15] used a probabilistic hyper-network to predict the direction of a stock price. Although theoretically globally trainable, the model’s implementation is univariate, making it again unsuitable for our dataset.

We first define the learning tasks as used in supervised learning classification works with similar feature sets [3] and similar learning tasks [4] to the TEH dataset. For a population of patients, $p\in P$, we have the dataset with features $f\in\mathcal{F}$, where each $f$ can be binary, $f_{b}\in\{0,1\}$, the value of some bound ordinal scale $s$, $f_{s}\in\{1,2,...,s_{max}\}$, or a value of continuous measurement, $f_{c}\in\mathbb{R}$. Examples of $f_{b}$ are general clinical data such as ‘history of drug use’ or an outcome on a specific psychometric measure, such as ‘irritability’ on the DASA assessment [1]. Each $f_{s}$ could represent the summation of the binary measure components of an assessment such as DASA, or a score on a Likert scale [19], while continuous measures, $f_{c}$, could refer to age or blood pressure at time of assessment. Each patient would have associated label, $y^{p}\in\{0,1\}$, representing the presence of a future event within a fixed future time frame. From this data we would learn a model that predicts a score or probability for the positive outcome (presence of the event) within the fixed time frame for any given patient, $\hat{y}_{p}\in\mathbb{R}|0\leq\hat{y}_{p}\leq 1$.

The TEH dataset consists of daily longitudinal psychometric and clinical data with a feature set, $\mathcal{G}$, that exhibits the same feature types as those described in $\mathcal{F}$, alongside a daily recording representing the presence of an instance of aggression over the same daily periods. We propose modelling the future presence of aggression using a global multivariable binary time series forecaster, where the target series is a patient’s binary time series indicating presence of aggression and associated covariate time series are features that vary over time. The task is, for any given patient, to forecast a score or probability associated with the presence of aggression in their next time interval. We note that the TEH dataset does not contain equal history length for each of the patients, contains static features that do not vary over time and is relatively small with only 83 patients total and a mean series length of 140 with a standard deviation of 106.

We formally define the new learning task as follows: For a patient $p\in P$, we represent the time series of the history of a dynamic feature, $d^{p,j}$ as the vector $\mathbf{d}^{p,j}_{1:Tp}=(d^{p,j}_{1},d^{p,j}_{2},...,d^{p,j}_{Tp})$ where $Tp$ is the total series length of a patient, and $j$ is a dynamic feature $j\in\mathcal{F}_{d}$. We then define the set of all of a given patient’s dynamic feature time series as $D^{p}_{1:Tp}=\{\mathbf{d}^{p,1}_{1:Tp},\mathbf{d}^{p,2}_{1:Tp},...\}$. The static features for a patient, $s^{p,l}$, $l\in\mathcal{F}_{s}$, are represented as the set $S^{p}=\{s^{p,1},s^{p,2},...\}$. The specific adverse event for a patient is represented by a binary vector $\mathbf{y}^{p}_{1:Tp}=(y^{p}_{1},y^{p}_{2},...,y^{p}_{Tp})$, where $y_{t}^{p}\in\{0,1\}$ indicates the presence of the event. The complete dataset then contains the sets: $D=\{D^{p}_{1:Tp}\}^{p\in P}$, $S=\{S^{p}\}^{p\in I}$ and $Y=\{\mathbf{y}^{p}_{1:Tp}\}^{p\in P}$. The objective is to learn a model $M$ that takes a fixed length window of length $n\leq Tp$ from a given patient to produce a probability associated with the future positive patient outcome, namely $\hat{y}^{p}_{Tp+1}=M(D^{p}_{Tp-n:Tp},S^{p},\mathbf{y}^{p}_{Tp-n:Tp})$, where $\hat{y}_{Tp+1}\in\mathbb{R}|0\leq\hat{y}_{Tp+1}\leq 1$. We note that when only a length of $n=1$ is considered by the model, this learning task matches the traditional future classification style approach with a bound of one day placed on the future time period of when the event can occur.

RNN-BOF uses a stacked Recurrent Neural Network architecture similar to that described by Bandara et al. [20], but with only the last LSTM cell’s hidden state vector from the last stacked layer being fully connected to a single sigmoid output neuron for the single probabilistic binary forecast point (Fig.  1). For input data, we concatenate the dynamic, static, and adverse event series elements per patient per time step, so an input from time $t$ is defined as $x_{t}^{p}=y^{p}_{t}\mathbin{+\mkern-10.0mu+}D^{p}_{t}\mathbin{+\mkern-10.0mu+}S^{p}$, where $\mathbin{+\mkern-10.0mu+}$ denotes vector concatenation and $x_{t}\in\mathbb{R}^{m}$ where $m=|\mathcal{F}_{s}|+|\mathcal{F}_{d}|+1$. For a fixed input window size of $n$, a stacked RNN architecture passes hidden state outputs from the layer $l$, $h_{t,l}$, as inputs to the next layer $x_{t,l+1}$, where $1\leq t\leq n$. For our RNN-BOF with $Z$ stacked layers, the output corresponds with the hidden state vector represented by $h_{n,Z}$, which is fully connected to a single neuron with sigmoid activation for learning the binary probability. This differs from other multivariate time series forecasters [20] where the last layer is instead fully connected to the output window for the learning of many continuous values. This also differs from stacked RNNs as implemented for binary classification tasks, where all hidden states from the final layer, $\{h_{t,l}\}^{l\leq Z}$, would be flattened and connected to a dense layer before again being connected to a single dense neuron with sigmoid activation. The output of the final neuron is defined by the equation: $P(y^{p}_{n+1}=1)=\sigma(\mathbf{w}_{d}\cdot h_{n,Z}+\beta)$, where $\mathbf{w}_{d}\in\mathbb{R}^{d}$ is the learned weight vector, and $\beta\in\mathbb{R}$ is the learned bias and $P(y^{p}_{n+1}=1)$ is the learned probability that the next day in the adverse event series of a given patient belongs to the positive class.

We make use of two techniques to deal with potential overfitting on small datasets. Dropout layers [21] are added after each LSTM layer (Fig.  1) which set values in the hidden layer output vectors, $h_{t,}$, to 0 at a fixed rate, and L2 kernel regularisation which has been used in RNNs as a tool to reduce overfitting [22] calculated by $L2=\frac{\lambda}{2}\sum w_{b}^{2}$, where $w_{b}$ is a specific hidden weight during back propagation and $\lambda$ is a learned hyper-parameter.

We train our model using a supervised learning sliding window training scheme, again similar to that described by Bandara et al. [20], but where the forecasting point is only the next value in the $\mathbf{y}^{p}_{t}$ series, which we will refer to as the label associated with any given input window. We employ this training scheme as it completely maximises use of the relatively small dataset. For each patient, $p\in P$, we define their ordered set of $n$-length dynamic feature windows as $\mathbb{D}^{p}_{1:Tp-n}=(D^{p}_{1:1+(n-1)},D^{p}_{2:2+(n-1)},...,D^{p}_{Tp-n:Tp-1})$, similarly, their adverse event windows are defined as $\mathbb{Y}^{p}_{1:Tp-n}={(\mathbf{y}^{p}_{1:1+(n-1)},\mathbf{y}^{p}_{2:2+(n-1)},...,\mathbf{y}^{p}_{Tp-n:Tp-1})}$, and static covariates are still the set $S^{p}$. Given windows, w, for patient p, $D^{p}_{w:w+n-1}$ and $\mathbf{y}^{p}_{w:w+n-1}$ have associated label $y^{p}_{w+n}$ where $w\in\mathbb{N}|w\leq(Tp-n)$. We can then define the label vector associated with $\mathbb{D}^{p}_{1:Tp-n}$, $\mathbb{Y}^{p}_{1:Tp-n}$ and $S^{p}$ as $\mathbb{L}^{p}_{1:Tp-n}=(y^{p}_{n},y^{p}_{n+1},...,y^{p}_{Tp})=\mathbf{y}^{p}_{n:Tp}$ (Fig.  2).

For training and evaluation purposes we split each individual patient’s windows sequentially into the train and evaluation sets instead of using a held-out patient group. This represents the real world case where a forecaster is implemented on a group of series to predict future values of those series and is common practice in global time series forecasting evaluation [9, 20]. We define this for a 80-20 train test split by taking the last $\rho^{p}$ windows per patient, where $\rho^{p}=\left\lfloor{0.2\cdot|Tp|}\right\rfloor$ so that across different window lengths $n$, the set of test windows will match exactly. The train and test window sets are then defined as $\mathbb{D}_{train}=\{\mathbb{D}^{p}_{1:\rho^{p}-1}\}^{p\in P}$, $\mathbb{Y}_{train}=\{\mathbb{Y}^{p}_{1:\rho^{p}-1}\}^{p\in P}$, with associated labels, $\mathbb{L}_{train}=\{\mathbb{L}^{p}_{1:\rho^{p}-1}\}^{p\in P}$ and $\mathbb{D}_{test}=\{\mathbb{D}^{p}_{\rho^{p}:Tp}\}^{p\in P}$, $\mathbb{Y}_{test}=\{\mathbb{Y}^{p}_{\rho^{p}:Tp}\}^{p\in P}$, with associated labels, $\mathbb{L}_{test}=\{\mathbb{L}^{p}_{\rho^{p}:Tp}\}^{p\in P}$. Thus, we maintain temporal integrity as no window labels are shared across the sets.

Unlike time series forecasting problems with continuous outputs where loss functions such as L2 or L1 are used to learn model weights, we use weighted binary cross entropy batch loss $wbl_{B}$:

We evaluate RNN-BOF on the Thomas Embling Hospital risk assessment dataset (TEH dataset) provided by Swinburne University of Technology³³3Ethics approval granted by Swinburne University of Technology (Project #1258). The dataset contains 40 static demographic and general clinical features, 73 features that vary daily or over longer periods, and daily recordings indicating the presence of an aggressive incident. The features include binary representations of general clinical data such as ‘history of drug use’, all components of START and DASA assessments, ordinal summations of both measures as well as general demographic data such as age.

We categorised the following three types of missing values found in the dataset: all patients had entries where the collection of dynamic features, psychometric assessments and presence of aggressive incidents spanned more than 24 hours; some patients had components of or full psychometric measures missing for some days; and a small subset were missing the label for presence of aggressive events. The frequency and length of the data entries lasting more than 24 hours were greater at the start and end of a given patient’s admission, with one or two long spans of daily assessments. We defined these long frequent assessment periods as a sequence of at least 30 days, starting and ending with 5 daily entries, and with an average of 1.1 days or less between entries over the entire period. As the contents of assessments lasting more than 24 hours were still accurate, no imputation or duplication was performed but a flag was added. For the dynamic missing values, nearest value linear interpolation was used and another flag was added, while patients that were missing values in the aggression series were dropped entirely. We performed z-score standardisation calculated per feature across the entire population’s training set. The cleaned frequent assessment periods were then windowed and split 80-20 into train and test window sets as described in Section III-C. Descriptive statistics for the dataset after pre-processing can be seen in Table I.

Previous works in the clinical psychology domain have commonly evaluated the predictive performance of their psychometric instruments using the ROC curve and its respective area [4, 5]. The ROC curve is not attenuated by imbalanced datasets such as ours and as such optimistically evaluates performance on the minority class, masking poorly performing models [23]. We instead use the Precision Recall Gain curve (PRG) [10], a variation on the Precision Recall (PR) curve which has been shown to be a stronger predictive evaluator than the ROC on imbalanced datasets [23]. While representing the same relationship between precision and recall, the Area Under the PRG Curve (AUC-PRG) is valid under cross-validation averaging and the curve does not require non-linear interpolation between points. The AUC-PRG is defined by the following equation:

We chose benchmark models that have seen use in clinical psychology machine learning papers as probabilistic or scored predictive classifiers. This list includes Random Forests (RF) [4, 5], Support Vector Machines (SVM) [3], Feed Forward Neural Networks (FFNN) [5] and logistic regression [4, 5]. We also included Light Gradient Boosted Machines (LGBM) due to their probabilistic output capability and for being strong global time series forecasters [13]. All benchmark models were trained using the same global sliding window training scheme across all patients as described in our modelling approach (Section III-C). The benchmarks were trained and tested using the same input window lengths of 10 and 20 days and therefore are also acting as binary times series forecasters. An additional one day length input window was also tested, representing the traditional classification approach described in Section III-A.

For psychometric instrument benchmarks we used the two measures collected in the TEH dataset; the DASA assessment [1] and the vulnerability component of START assessment [2]. The DASA total score is calculated by summing seven individual binary scored components, each representing the presence of a risk indicating behavior. The START vulnerability assessment contains 20 components, each scored 0, 1 or 2 on a scale representing vulnerability that is summed to give a score between 0 and 40.

The complete set of windows were split sequentially 80-20 per patient, then pooled as the train and test set respectively (Section III-C). For each model’s hyper-parameter tuning we used 5-fold expanding sequential cross-validation on the pool of patients’ training windows [24], with folds instead of single instances and pooling across a population. The cross-validation procedure splits the training windows per patient into six sequential and equal size window pools. A validation model is then trained on the first pool and validated on the second, and a second validation model is then trained on the first two pools and validated on the third and so on. This still maintains the temporal integrity described by Hyndman and Athanasopoulos [24] as no labels are shared between train and validation folds. The AUC-PRG is averaged across validation folds to inform a Bayesian hyper-parameter tuning algorithm (namely the Hyperopt algorithm, for 100 iterations). After parameter selection, each model is trained on the full training set [9] with forecasts then performed on the held-out test windows for final performance evaluation.

We chose the hyper-parameters and their ranges through existing literature [21, 22, 9] and manual exploration on the training set cross-validation folds. The tuned parameters and their respective ranges identified for RNN-BOF were: number of hidden nodes per layer, 10 - 100; number of stacked layers, 1-4; dropout percent, 10% - 50%; L2 regularisation, 0.000001 - 0.001; and number of epochs, 10 - 100. We used a fixed batch size of 512 and found that too high of a learning rate led to especially poor performance so we fixed it at 0.0001. We note that all benchmarks except the FFNN showed worse performance with the weighted loss despite the imbalanced nature of the output windows, so standard binary cross entropy loss was used.