{'fraction_non_alphanumeric': 0.05321703956516243, 'fraction_numerical': 0.012893439514599925, 'mean_word_length': 4.528651292802236, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust autoregressive hidden semi-Markov models applied to EEG sleep spindles detection

We propose a generative model for single-channel EEG that incorporates the constraints experts actively enforce while visually scoring recordings. In particular, the framework takes the form of a robust hidden semi-Markov model that explicitly segments sequences into local, reoccurring dynamical regimes. Unlike typical detectors, our approach takes the raw data (up to resampling) without any filtering, windowing, nor thresholding. This not only makes the model appealing to real-time applications, but it also yields interpretable hyperparameters that are analogous to known clinical criteria. We validate the model on stage 2 non-REM sleep recordings that display characteristic sleep spindles. We derive tractable algorithms for exact inference and prove that more complex models are able to surpass state of the art detectors while being completely transparent, auditable, and generalizable.


INTRODUCTION
Sleep spindles are a hallmark of stage 2 non-REM sleep. They are the result of mutual interactions between GABAergic reticular neurons and excitatory thalamic cells [1]. Their proposed functions include memory consolidation [2], cortical development, and potential biomarkers for psychiatric disorders [3]. Hence, proper detection and modeling are crucial.
A sleep spindle is defined as an oscillatory burst in the range 11-15 Hz (sigma band) with duration between 0.5 and 2 s., distinctive waxing-waning envelope, and maximal in amplitude using central Electroencephalogram (EEG) electrodes. EEGers usually adhere to clinical manuals to visually score EEG traces [4,5] (e.g., multi-channel recordings are sequentially inspected in 30-second-long epochs and further categorized according to the estimated temporal dynamics and observed reoccurring patterns). EEGers are also trained to identify non-brain-related artifacts (e.g., eye and muscle This work was supported by NSF CAREER grant 1453756. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. (http://www.tacc.utexas.edu). Code to replicate our results can be found in https://github.com/carlosloza/spindles-HMM activity). Nowadays, with an ever-increasing amount of data, it is advantageous to exploit signal processing techniques and machine learning models to guide principled, human-like automatic EEG scorers.
Typical automatic sleep spindles detectors comprise four stages: pre-processing, decomposition, decision making, and feature extraction [6]. The first stage usually involves bandpass filtering and artifact rejection, while the second stage applies wavelets or other windowing techniques for pseudostationary processes. Decision making takes the form of hard-thresholding plus cross-validation using ground truth (i.e., available labeled EEG). Lastly, feature extraction estimates distinctive traits of sleep spindles for further analysis. A myriad of approaches use slight variations of this pipeline [7,8,9,10] (see [6] and [11] for reviews).
Even though this common methodology provides adequate results, it suffers from two main drawbacks: i) analysis based on filtered EEG, and ii) lack of theoretical foundations for hyperparameter settings. The first point not only differs from human-like scoring, it is also computationally expensive for online applications. The second point is a consequence of a model-less framework; that is, amplitude thresholds are data (and scorer) dependent, which severely limits generalization. In fact, most sleep spindles detectors can be cast as EEG transforms with a large number of hyperparameters (filter, window coefficients, thresholds).
Instead, we propose a generative model with the following components: i) robustness against artifacts via observations modeled as linear generalized t likelihoods, ii) reoccurring modes with distinctive dynamics that characterize non-spindle and spindle regimes, iii) semi-Markovian states that generalize the geometric regime durations of Hidden Markov Models (HMM), and iv) tractable exact inference via message passing routines. The result is a model that encodes the constraints EEGers actively enforce when visually scoring EEG: a robust autoregressive hidden semi-Markov model [12].
We validate the model on the DREAMS sleep spindles database [13] and achieve comparable or better performance than the state of the art. The results open the door to more complex models and sophisticated inference in the future (e.g., real-time detection). The rest of the paper is organized as follows: Section 2 details the model and learning, Section 3 presents the results, and Section 4 concludes the paper.

GENERATIVE MODEL FOR EEG
Let y = {y n } N n=1 be an ordered collection of N samples that represents random variables over time. For our case, y is a single-channel, single-trial EEG trace. The proposed generative model poses each y n as a linear combination of p previous samples plus a noise component parameterized as a zero-mean generalized t distribution (i.e., linear generalized t model). Moreover, the specific parameters of y n (linear weights, scale and degrees of freedom of noise) depend on a bivariate discrete hidden state comprised of labels, z n , and counters, d n , where the former denotes the dynamical mode of the network (out of K possible regimes) and the latter keeps track of the remaining duration in said mode. The generalized t noise accounts for model robustness whereas the bivariate state allows for regime durations beyond the implicit HMM geometric paradigm (i.e., z n is semi-Markovian). Additionally, the autoregressive structure of the observation model is key to capture local temporal dependencies, such as oscillatory sleep spindles. The result is a robust autoregressive hidden semi-Markov model (RARHSMM), which is depicted as a dynamic Bayesian network [14] in Fig. 1.
We adopt a 1-of-K notation for mixture models as in [15]: z n and d n represent K and D-dimensional multinomials, respectively, where only one entry is equal to 1 and the rest are 0 (e.g., z n,k means the k-th entry of z n is 1). The initial hidden label, z 1 , is parameterized by a vector of probabilities, π, while the subsequent z n obey the semi-Markovian dynamics: where A is a K ×K stochastic matrix, c n = d n,1 , and δ(·, ·) is the Kronecker delta function. In words, the hidden label, z n , transitions to a different regime only when the counter variable is equal to 1; otherwise, it remains in the same mode. Similarly, d n samples a new duration only after a regime ends; otherwise, it decreases by one unit. This model is also known as explicit duration hidden semi-Markov [16].
where λ is a collection of K sets of parameters that parameterize the regime durations (e.g., rates of Poisson variables). On the observations side, a convenient parametrization for robust purposes is the Normal-Gamma distribution for the pair (τ n , y n ) (i.e., τ n is marginally Gamma distributed and y n is conditionally Normal given τ n -this results in y n being marginally generalized t distributed). In equations: y n |τ n , z n , σ, a, y n−1 n−p ∼ N a k , y n−1 n−p , y n |z n , σ, ν, a, y n−1 n−p ∼ t a k , y n−1 n−p , σ k , ν k (5) where ν, σ, a are collections of K degrees of freedom, K scale parameters, and K p-dimensional weights, respectively. y n−1 n−p denotes the p observation preceding y n , and ·, · is the inner product. The previous parametrizations and the conditional independences implied by the graphical model in Fig. 1 yield the following complete data log-likelihood over model parameters, θ = {λ, A,π, ν, σ, a}, (ignoring constants): where p(d n |λ k ) is shorthand for p(d n |z n,k , λ), Γ(·) is the gamma function, A j,k denotes the j-th row, k-th column entry of A, and z = {z n } N n=1 (similar for d and τ ). Three tasks are of particular interest: a) given θ, calculate the marginal log-likelihood (LL) of the observations y (i.e., p(y|θ) via marginalization of p(y, τ , z, d|θ)), b) given y and θ, estimate the most likely sequence of hidden states, and c) given y (or a batch of sequences), estimate θ. We refer to the tasks as loglikelihood calculation, inference, and learning, respectively.

Learning model parameters
Maximizing eq. (6) directly is intractable; rather, we exploit a variation of the Expectation-Maximization (EM) algorithm for HMM [17]. We implement coordinate ascent on the loglikelihood: first, we keep θ fixed and obtain expectations of log L(θ) under the posterior distribution of the latent variables given y (underlined variables in eq .(6)). Then, said expectations are fixed and used to optimize the θ components. Both steps constitute one EM iteration; the algorithm continues until convergence to a local optimum (e.g., tolerance between sucessive LLs). For simplicity, we outline the case of a single sequence, yet, the batch case is straightforward [17].

E-step
Let α(z n , d n ) = p(y n 1 , z n , d n ) be the joint probability of the first n observations and the n-th hidden state, and similarly β(z n , d n ) = p(y N n+1 |z n , d n ) . By means of induction and applying d-separation [14] to the graphical model, the following is straightforward (ignoring parameters to avoid clutter): These estimators are the bedrock of all tasks. For example, inference replaces the sums in eq. (7) with "max"" operators plus backtracking (Viterbi algorithm [18]), and LL calculation reduces to compute p(y) = zn dn α(z n , d n ) . Next, by applying Bayes theorem and d-separation: η(z n , d n ) = p(z n , d n |y) = α(z n , d n )β(z n , d n ) p(y) γ(z n ) = p(z n |y) = dn η(z n , d n ) ξ(z n , z n−1 , d n ) = p(z n , z n−1 , d n , d n−1,1 |y) = α(z n−1 , d n−1,1 )p(y n |z n )p(d n |z n )p(z n |z n−1 )β(z n , d n ) p(y) where p(d n |z n ) and p(z n |z n−1 ) refer to regime switching cases. By construction, the entries of z n and d n are binary; hence, their expectations are the probabilities of taking the value 1 (color coded using the same scheme as in eq. (6)). By conjugacy, the posterior conditional distribution of τ n is Gamma with mean ω n,k ; this is used to derive the expectation of log τ n as well [19] (Ψ(·) is the Digamma function): E τ |y,z,θ {log τ n } = log ω n,k + Ψ ν k + 1 2 − log ν k + 1 2 (13)

M-step
Now, we plug in the previous expectations into eq. (6) and optimize each θ component (partial derivatives equal to 0).
where (t) denotes estimates for the t-th EM iteration and γ(z n,k ) = E{z n,k } (analogous for η and ξ). W is a diagonal matrix with entries γ(z n,k )ω n,k . Y p ∈ IR (N −p)×(p) is an embedding matrix of y. Therefore, the optimal a (t) k is calculated via weighted least squares with weights that incorporate label assignment probabilities and robustness. λ k is the solution to the following non-linear problem:

RESULTS
The DREAMS sleep spindles database [13] was used to validate the proposed model. 30-minute-long, single-channel (CZ-A1 or C3-A1) recordings from 8 subjects with corresponding sleep spindles expert scores (2 experts) are available. All traces were resampled to 50 Hz and z-scored prior to feeding them to the model. No bandpass filtering, artifact rejection, nor windowing were implemented. The data was then partitioned into 8 folds of training and test sets (i.e., 7 subjects for training, 1 for testing). The performance measure is the by-sample Matthews Correlation Coefficient (MCC) between the inference output on the test set and the corresponding expert scores. We chose this strict measure to avoid any ambiguity in our results (in contrast to event-based measures). We detail three approaches: a) a supervised scheme in which complete data is used to estimate model parameters (i.e., no EM needed), b) an unsupervised framework in which no labels are provided for learning, and c) an expert-specific supervised approach that allows to compare visual scores.  The supervised case takes the union of both expert scores as ground truth and estimates θ by replacing all M-step expectations with dirac delta probabilities. We set K = 2, (only two regimes are of interest) and p = 5 in order to capture sigma band oscillations. Lastly, we set D = 30 (in seconds) in addition to non-zero self-transitions for the non-spindle mode; this results in a tractable model that also resembles the well-known clinical 30-second-epoch scoring paradigm [4]. Table 1 summarizes the results compared to two previous efforts [8,9]. RARHSMM outperforms both methods also when the by-sample F1 score is the performance measure (mean=0.451 compared to 0.370 and 0.401, respectively). The database authors [7] report an event sensitivity of 0.702 and false positive rate 0.014, while we achieve 0.764 and 0.039, respectively. In [10], the authors use window-based spectral methods to suppress artifacts in a pre-processing stage and achieve an average MCC of 0.447. Therefore, our proposed method is able to perform feature extraction and detection while handling outliers in a principled manner.
For the unsupervised case, we provide sensible initial conditions: i) π = [1, 0] ⊤ , ii) A = 0.5 0.5 1 0 (no self-transitions for spindles), iii) non-spindle autoregressive (AR) coefficients estimated from the first few seconds of EEG, iv) spindle AR weights with a peak in the sigma band, and v) uniform durations for non-spindles and N (1 s., 0.15 s.) for spindles.
The estimated observation parameters (Fig. 2) show clear differences between modes: the non-spindle regime depicts the well-known 1/f spectral distribution [20] while its coun- terpart has clear oscillatory activity in the sigma band. Also, the heavier tails for the non-spindle regime (ν 1 = 4, ν 2 = 9) confirm the expert opinion that most artifacts appeared during non-spindle epochs. We obtain an average MCC of 0.3796, which should be taken with a grain of salt because the supervised case has the advantage of parameter settings and ground truth set according to experts; conversely, the unsupervised case is purely generative and data-driven. As a side note, the approach in [21], autoregressive HMM with Gaussian observations, achieves an average MCC of 0.071. Fig. 2 also details the non-parametric distribution over sleep spindle durations. Its unimodal, skewed nature has also been reported in large-scale studies [11]; however, here λ is an explicit component of the generative model, whereas in classic detectors, it is simply a byproduct of the transform.
The generative framework also allows for more in-depth analysis of the expert visual scores. We fit two different models using complete data according to the labels from each expert, infer the most likely labels on the test set, and evaluate the MCC with respect to the corresponding expert. The expert 1 model achieves an average of 0.351 in comparison to 0.481 for expert 2. This can be interpreted as a consistency measure for visual scorers and opens the door to more sophisticated tools, such as boosting and Bayesian model averaging.
Lastly, different local solutions can be compared also. Table 2 details predictive NLL and shows that the parameters set by expert 2 have stronger predictive power (under the model in question). Also, the results suggest that the unsupervised case outperforms the supervised scheme when it comes to reconstruction error, contrary to what the MCC advises.

CONCLUSION
Probabilistic graphical models over sequences are principled frameworks for robust detection and parametrization of EEG sleep spindles. Future work will go beyond the sigma band to characterize and detect high-frequency rhythms in other structures when ground truth is not available, e.g. hippocampal gamma oscillations and sharp wave-ripples [22].