The analyses outlined in this study were undertaken using so-called discrete-time Poisson hidden-state Markov models. The basic notion is that the observations are a discrete time-series consisting of seizure counts, observed on a number of consecutive occasions. The central premise is that the observations emanate from an underlying system that switches between states, as illustrated in Figure 1. In some applications this assumption might be based on established theory. An alternative motivation for adopting an HM model for count data is as an empirical approach to dealing with overdispersion, i.e., the condition where the count data has a variance that is greater than the expected variance under the Poisson model. Although the present application belongs to the latter category, the possibility remains that an HM process is a realistic representation of the underlying physiological process that gives rise to the observed seizure behaviour. The key idea is that the observed counts are all that is available in terms of observable data and that the underlying process is hidden.

The HMM assumes that the hidden state switches (transitions) occur according to a probabilistic rule. More specifically it is assumed that the hidden states behave as a discrete-time Markov process. A Poisson HMM assigns differing Poisson rates to each of the underlying hidden states. Thus, as the system switches between states, the rates at which the observed events occur (seizures in the present application) are altered, hence the overdispersion. The challenge is to extract information about the hidden Markov process given the observed counts. Additional information is given in Note S1, which provides a more formal description.

An important feature of epilepsy seizure-count HM modelling is the potential to provide information regarding the state of the subject on a given occasion. In speech recognition and other applications this is referred to as decoding. Two types of decoding algorithms are available, namely local and global. Local decoding provides state probabilities at each occasion, while global decoding determines the sequence of states that is most likely, given the observations. This paper includes decoding results based on the so-called expectation-maximisation (EM) algorithm ⁸ , together with posterior seizure-state probabilities provided by Bayesian HMM analyses. The resulting seizure-state time series data were used to generate seizure-state duration (sojourn time) probability densities and associated cumulative distributions.

The primary objective for this study is the development of statistical models that provide information that is useful from a clinical perspective. Thus the main focus is a characterisation of the seizure states that, according to the model, underlie the observed seizure behaviour. This is motivated by the belief that a model-based description of the process will provide a clue to future behaviour. From a patient management perspective seizure-state duration (sojourn time) prediction is especially relevant, prompting us to focus on sojourn-time cumulative probability distribution results. Thus the main objective is the provision of model parameter estimates. These take the form of a point estimate, usually a median or mode, together with some measure of uncertainty.

Two alternative kinds of HMM analyses have been used in this work, namely traditional likelihood analyses performed using the EM algorithm, and Bayesian analyses implemented using MCMC ^{9, 10} . The main reason for using these two rather different approaches was to assess robustness. Each of the seizure time series reported in this study is relatively short in duration. Parameter estimation is, therefore, compromised because a short time series will include relatively few state transitions, rendering it difficult to obtain the HMM transition probability estimates with good precision. In clinical practice, decisions must sometimes be taken within a short time frame and our approach reflects this clinical reality. By comparing the results provided by two or more differing models/analyses, an indication of robustness is obtained. We stress that, in contrast to traditional model comparison analyses, where the aim is to find the best model, the objective in this study is different. Rather than declaring a given model as being ’best’, the objective is to determine whether related models yield similar parameter estimates. Agreement is taken as an indication of robustness. Sampling methods were used extensively, particularly as a mechanism for obtaining measures of uncertainty. In addition to using MCMC to perform the Bayesian analyses and thus obtain the required Bayesian credible intervals, parametric simulation ¹¹ was used to generate confidence intervals associated with various EM estimates.

In the Bayesian context MCMC is, in effect, a random sampling procedure for evaluating complicated integrals in high dimensional space, as required to calculate various summary statistics. It is widely used in Bayesian analysis, the Gibbs sampler being particularly prominent among the available methods. Several software packages have been developed for performing Gibbs sampling, including WinBUGS ¹² , which was used in the present work (version 1.4.3). Presented with the data, a ‘measurement model’ (likelihood) and prior distributions, WinBUGS constructs the samplers required to generate the MCMC output. In order to avoid confusion we should emphasise that the decision to use MCMC to perform the analyses is entirely unrelated to the fact that the statistical model is based on a Markov process. MCMC is just one of several methods that might be used for HMM parameter estimation, the EM algorithm providing an alternative.

A particular concern in the present study is the adoption of meaningful HMMs. Reasonable care must be taken to ensure the validity/adequacy of these models. Posterior predictive simulation is a Bayesian model checking procedure which is especially useful in modelling applications in which standard model checking methods, based on residuals analysis for example, are problematic. Discrete-data (counts, etc.) regression modelling is an obvious example ¹³ . The method is particularly useful for exploring potential model failure specific to the key research questions. The basic idea underlying posterior predictive simulation is the generation of replicate datasets which, in the MCMC setting, is achieved by using the posterior sample generated under the model in question. One or more key features of the observed data are compared with the replicate datasets. If these features are not seen in a majority of the replicates, this is taken as an indication of model failure. The comparison between the simulated and observed data can be performed using some discrepancy measure or can be undertaken using graphical methods. The latter approach has been suggested as particularly useful in practice ^{9, 14} . In the present epilepsy seizure study the occurrence of extreme counts (either low or high) is an obvious candidate for scrutiny. The model needs to capture these extremes. It is important to note that the replicate datasets incorporate both the uncertainty specified under the model (Poisson variation in the present context) and parameter uncertainty. It might also be noted that the replicate datasets fulfil a dual role in the present study. In addition to their model checking role they provide the clinician with an indication of what might be expected in terms of future seizure pattern in the patient under consideration, assuming no indication of model inadequacy. A key feature of the posterior predictive simulation assessment method is a probabilistic measure of model inadequacy. In this paper we use Bayesian p-values, which are a common choice. Additional information is provided in Note S2.

Hidden Markov models have a history going back to the 1960s, long before Bayesian statistics became accepted by mainstream data analysts. Although HMM parameter estimation can be tackled using standard minimisation (maximisation) methods, the EM algorithm is commonly adopted, based on well-established literature. The monograph by Zucchini and MacDonald (2009) ⁸ , provides a detailed account, together with a set of routines for implementing an HMM analysis based on both minimisation and the EM algorithm. A number of R packages also provide tools for performing HMM analyses, including the R package msm ¹⁵ .

Although maximum likelihood estimation in HM modelling is well established, the traditional approach to confidence interval and standard error estimation is to resort to asymptotic (large sample) approximation. The resulting intervals/standard error estimates can be inaccurate, however, and an alternative strategy is to avoid asymptotic approximation by adopting a sampling method. Zucchini and MacDonald (2009) [ 8, page 53], suggest the parametric bootstrap. In this study we have adopted a sampling approach referred to as parametric simulation [ 11, Chapter 2].

Young Epilepsy is the United Kingdom’s leading provider of education, treatment and residential care for children and young people with severe epilepsy. Located in Surrey just outside Lingfield, it has on-site 24-hour medical services. All young people have very close supervision by carers throughout the days and nights and these carers compile seizure diaries if seizures are observed. (Seizure diaries are used in similar organisations to Young Epilepsy and have been used in previous studies to monitor seizure frequency and severity ¹⁶ .) Evaluation of the diaries plays a major role in carers and medical staff assessing whether there have been substantial changes in seizure frequency and or severity. Thus, much detail and attention is given to their completion and accuracy. Total weekly seizure frequencies over a period of 1 year (starting in September 2004) were ascertained from the diaries for all the participants.

The project was registered at Young Epilepsy as a case note review. All applications for registration are reviewed by a formally constituted internal committee consisting of physicians, nurses and lay people. Questions related to ethics that are raised at that committee are referred to the East Surrey Ethics Committee. The local committee raised no ethical questions about our study. In addition the data collection was deemed to be audit and, under the charity’s official rules, consent was not required. The data were not formally blinded because the study does not take the form of a clinical trial. Instead, the focus is the development of a statistical model, the objective being parameter estimation as an aid to clinical decision making. Thus the purpose is parameter estimation as distinct from hypothesis testing. That said, the analyses were performed by MDK after extracting the seizure time-series data from the original clinical records, thus removing all sources of identification. Prof Rod Scott and Dr Suresh Pujar had access to the clinical information as part of routine clinical care, but took no part in the development of the statistical models or data analysis.

Three children with epilepsy were selected for this study based on an examination of seizure frequency charts, subject to the restriction that no change in medication had occurred during the observation period (inclusion criterion). The three cases were chosen because they differed markedly in their time-dependent seizure behaviour, each presenting a particular challenge from a statistical modelling perspective. The central point to note regarding the choice of patients is the focus of this paper, namely seizure-count time series modelling, aimed at the provision of individual-specific parameter estimates, as opposed to group-average estimation. Our premise is that the analytical procedure can be applied to any individual, recognising that the model details must be tailored to the individual under investigation. This premise is supported by the demonstration that the method yields useful statistical information when applied to patients with disparate clinical presentations.

Two distinct Bayesian seizure-count time-series models were adopted in this study, namely a log-linear HMM and a Poisson-gamma HMM. The details are as follows. It might be noted that Congdon (2006) ¹⁷ provides relevant information on Bayesian HMMs, including their implementation in WinBUGS. It might also be noted that HMMs are a class of mixture model, and that mixture model parameter estimation is not always straight-forward due to a phenomenon known as label switching. The label switching problem is discussed briefly in Note S3.

The following statistical model descriptions are given using an informal notation based on the WinBUGS distributions used in their implementation. This serves the dual purpose of providing a more intelligible outline for a non-mathematical readership and, at the same time, gives an outline of the key components of the WinBUGS code used to implement the models. The observed seizure counts ( y _t ), conditional on state ( Z _t ), are assumed to have a Poisson distribution. In the J-state case and assuming time-independent state-specific Poisson rates (and noting that several model parameters are completely defined within the hierarchical set of equations that follow and do not require further definition)

                                     y _t     ~     Poisson( λ _{Z t}  ), t = 1,...,  N                                                 (1)

                            log( λ ₁ )    ~     N(0, ϵ _p ),  l <  λ ₁ <  u                                                           (2)

                             log( λ _j  )    ~     N(0, ϵ _p ),  λ _j–1 <  λ _j  <  u,  j = 2,...,  J                                 (3)

                                    Z _t     ~     Categorical( P _{Z t–1,1: J} ), t = 2,...,  N                                  (4)

                                    Z ₁    ~     Categorical( P _{0,1: J} )                                                              (5)

                                P _{i,1: J}     ~     Dirichlet( 1 ), i = 0,...,  J                                              (6)

where the symbol ~ indicates ‘distributed as’, N (., .) is the normal distribution with the specified mean and precision, t is an observation occasion label, N is the number of observations, ϵ _p is a small value, l and u are lower and upper constraints, P _{Z t−1,1: J} is a transition matrix row, P _{0,1: J} is a J-length vector of probabilities with sum equal to one, j is the state label, J the number of hidden states and 1 is a J-length vector of ones. Categorical( p) is the distribution of the N mutually exclusive states, as determined by the class probability vector p. (It is implemented in WinBUGS by the dcat distribution (page 352 in 18), which ‘returns’ one of the J states at each MCMC iteration.) The transition matrix row priors are given as Dirichlet distributions mainly for clarity, although the majority of the simulations performed in this study adopted an equivalent Gamma( alpha, 1) construction (pages 235 and 351 in 18). In those cases where the state-specific Poisson rate parameters are a function of time the observation model becomes

           y _t | Z _t , ( λ _{1 t} , . . . λ _{J t} )     ~      Poisson( y _t | λ _{Z t t} ), t = 1,..., N                                          (7)

                             log( λ _jt )     =      λ ^′ _j + f _t , j = 1,..., J                                                            (8)

                             λ ^′ ₁     ~      N(0, ϵ _p ), l < λ ^' ₁ < u                                                             (9)

                              λ ^′ _j      ~      N(0, ϵ _p ), λ ^′ _j–1 < λ ^′ _j < u, j = 2,..., J.                                (10)

Thus, f _t adds a time-dependent offset to each of the Poisson rate parameters in the log domain.

Among the approaches adopted for incorporating time dependence in the state-specific Poisson means is the introduction of quadratic dependence ( f _t = β ₁ t + β ₂ t ²), the introduction of a temporal random effects term and the incorporation of a first order autoregressive term (abbreviated AR(1)). The AR(1) and random effect model details are given in Note S4.

The following Poisson-gamma model was adopted as an alternative to the log-linear model in those cases where the state-specific Poisson rates are time-independent. Using Poisson rate-parameter prior assignments based on George et al. (1993) ¹⁹ and Gelfand and Smith (1990) ²⁰ , the observation model becomes

                                                y _t     ~     Poisson( λ _{z t} ), t = 1,..., N                                           (11)

                                                λ _j     ~     Gamma( α _j , β _j ),   j = 1,..., J                                           (12)

                                                α _j     ~     Exp(1.0), j = 1,..., J                                                   (13)

                                                β _j     ~     Gamma(0.1, 1.0), j = 1,..., J                                        (14)

                                                Z _t     ~     Categorical( P _{Z t–1,1: J} ), t = 2,..., N                           (15)

                                                Z ₁    ~     Categorical( P _{0,1: J} )                                                    (16)

                                            P _{i,1: J}     ~     Dirichlet( 1),   i = 0,..., J.                                        (17)

Hidden-state sojourn times (also referred to as dwell times) are key to the clinical problem. Many parents and patients recognise the fluctuating nature of epilepsy, but find it difficult to judge how long the child will be in any given state. Sojourn-time statistics enable the clinician to provide this important clinical information. The sojourn-time density for the jth state, j = 1, …, J is given by ²¹ d j ( s ) = P ( Z t + s + 1 ≠ j , Z t + s − v = j , v = 0 , … , s − 2 | Z t + 1 = j , Z t ≠ j ) ( 18 ) where J is the number of states, s ϵ {1, …, M _j }, and M _j is the jth state sojourn time upper limit. Thus the sojourn time, s, is the number of observation occasions t + 1, …, t + s, during which the unobserved process remains in the jth state, conditional on being in some other state on occasions t and t + s + 1. Sojourn time probability distributions were derived from the MCMC output obtained for the transition matrix diagonal elements by substitution into the geometric distribution. Specifically, probability densities were generated using g ( P i i ( k ) ) , k = 1 , … , K , where g ( P i i ) = P i i s − 1 ( 1 − P i i ) and P i i ( k ) , k = 1 , … , K is the Gibbs marginal posterior sample obtained for the ith diagonal element of the transition matrix. ( K is the MCMC sample size and ( k) the iteration label. This is based on the standard approach to evaluating the marginal posterior distribution of a variable that is a function of the model parameters; see 22, for example.) To provide the results in a form that will be more familiar and accessible to many clinicians, the duration probability statistics are given as complementary cumulative (probability) distributions, i.e., the probability that the sojourn time is greater than a given value. It might be noted that these complementary cumulative probabilities are equivalent to the survivor function that arises in survival analysis.

Gibbs sampling, convergence analysis and posterior predictive simulation. Gibbs sampling was performed using WinBUGS/GeoBUGS ^{12, 23, 24} , which was downloaded from http://www.mrc-bsu.cam.ac.uk/bugs. The dcat (categorical distribution) and ddirch (Dirichlet distribution) are among the WinBUGS distributions used to implement the statistical models outlined above. Truncated normal samples were generated using the WBDev function djl.dnorm.trunc ²⁵ . Three parallel chains were generated in each analysis, each chain consisting of 5000 samples after thinning. (Thinning is the name given to the common practice whereby a specified proportion of the MCMC output is discarded, the remaining samples being stored for subsequent processing. It has been discussed by a number of analysts, including Carlin and Louis (2009) (Section 3.4.5 in 10). The reason for thinning the MCMC output is to produce a chain with reduced sample autocorrelation. The aim is to reduce both the storage and post-simulation CPU demands of an analysis without suffering much loss in precision. A thinning factor of 10 was used in the majority of the simulations. Thus, 1 in 10 samples was stored and used in subsequent calculations.) The first of each set of three chains was started at an arbitrary position in parameter space, while the other two were started at over-dispersed positions. A burn-in set of samples was acquired prior to storing each chain of 5000 samples. These burn-in samples were also discarded.

Convergence tests are an essential part of any MCMC analysis. Nevertheless, given the exploratory nature of this study in which many models were investigated for each of the three cases, a complete set of convergence tests on all models and all parameters (i.e., both primary and derived parameters) was not undertaken. Instead, the convergence tests were restricted to the log-linear model MCMC output obtained for the primary parameters, i.e., the Poisson rates and transition matrix elements. This provides protection against a complete failure to achieve an adequate sampling of parameter space. With this safeguard in place, a comparison of the statistics provided by the alternative hidden Markov models adopted in each case, e.g. log-linear compared with Poisson-gamma, provides a mechanism for assessing the robustness of the results.

Informal convergence assessment was performed by inspecting overlaid trace plots for visual signs of failure. This was followed by semi-formal analyses (see 26 for a review of methods), which were performed using three convergence test procedures, namely the Gelman-Rubin shrink factor diagnostic and associated shrink factor plots (which is based on an ANOVA-like assessment of the between-chain and within-chain variances), the Geweke Z-score diagnostic and Z-score plots (based on spectral density variance estimation and a Z-score comparison of chain segments), and the Raftery-Lewis diagnostic procedure (which provides a variety of data, including an estimate of the number of iterations required to obtain a given quantile to some specified accuracy, taking into account the correlation between samples). These convergence analyses were performed using the R CODA (Convergence Diagnosis and Output Analysis software) package ²⁷ (R version 2.15.2, CODA version 0.16–1).

As stated in the Introduction, posterior predictive simulation was used extensively to assess the capacity of the various Bayesian models to account for the observed seizure-count time-series data. Bayesian p-values were adopted as a measure of model inadequacy. The procedure was implemented as follows. At each MCMC iteration a simulated dataset (commonly referred to as a replicate dataset) was generated using the parameter values obtained at that iteration. The resulting replicate datasets and the observed data were used to generate a Bayesian p-value estimate, as given by the proportion of replicates satisfying equation 19, based on the chosen discrepancy measure. For example, among the tests that were performed is the capacity of a given model to capture the observed maximum seizure rate for the subject in question. In this case p _B = #( T( y ^rep ) ≥ T( y)) /K, where K is the number of MCMC iterations and T( y) is the number of occasions on which the observed maximum rate occurred. Additional information is given in Note S5, Note S6 and Note S7, while Note S2 provides a brief introduction to posterior predictive simulation.

The EM algorithm and parametric simulation. The majority of the analyses undertaken in this study were performed using Bayesian statistical models, as outlined above. In some instances, however, traditional (non-Bayesian) HMM EM analyses were also undertaken for comparison, the results of which are shown for one of the three cases. These EM calculations, which were restricted to those analyses where the state-specific seizure rates (i.e., Poisson rates) are time invariant, were implemented using the R routines provided by Zucchini and MacDonald (2009) ⁸ . Parametric simulation [ ¹¹ , Chapter 2] was used to construct the associated confidence intervals. (Parametric simulation is a random sampling procedure in which simulated data are generated using the fitted model. Parameters of interest are recalculated from each simulated dataset. Confidence intervals and other statistics can be constructed using the resulting parameter sample distributions.)

All overlaid trace plots were effectively ideal in their visual appearance except for some Poisson rate-parameter plots which showed random and very short-lived excursions to high values. But these excursions were observed in all three chains when they did occur, suggesting satisfactory coverage of the posterior distribution. The Gelman-Rubin shrink factors and associated plots were mostly exemplary; the remainder were all satisfactory. The Geweke Z-scores were also satisfactory, although the Geweke Z-score plots invariably showed one or more scores greater than 2.

The Raftery-Lewis analyses indicated that in two of the three cases, a sample of size 5000 (after thinning) was sufficient to provide the majority of 0.025 quantile estimates with an accuracy of ±0.005, or better, with probability 0.95. Thus the resulting nominal 95% credible intervals have a true coverage of between 94% and 96%, with probability 0.95. In a few instances the accuracy was slightly less, but mostly better than ±0.0055. In Case 1 this level of accuracy was not achieved for every parameter. Nevertheless the accuracy was at least ±0.01, which provides nominal 95% intervals with a true coverage of between 93% and 97%, with probability 0.95.

The seizure count time series observed in the first case is unremarkable except for a single occasion when the rate was close to double that observed on the majority of occasions ( Figure 2). The observed variance suggests overdispersion while a comparison of a 2-state HMM and a standard Poisson model is ambiguous, AIC selecting the 2-state HMM, while BIC suggests that the standard Poisson model is adequate. We note, however, that under the Poisson distribution with the observed mean value, the probability of observing a seizure count of 13 or higher in a series of 53 observations is close to 0.05, providing justification for an examination of the 2-state HMM. Three separate 2-state HMM analyses were undertaken, the first using the EM algorithm combined with parametric simulation confidence interval estimation, the second and third using the MCMC implementation of the Bayesian log-linear and Poisson-gamma models, respectively.

Posterior predictive simulation was performed as part of each MCMC analysis. This provides a powerful method for assessing the potential inadequacy of the model under consideration, as determined by a lack of capacity to capture key features of the data. Visual inspection of graphical displays and Bayesian p-values, both derived from the simulated data, are among the tools used for model evaluation. Similarly, parametric simulation was used to produce equivalent non-Bayesian p-values using the EM parameter estimates. As stated in the Introduction, and taken up in the Discussion, the aim of this exercise is not to compare the alternative models/analyses with a view to selecting one as being preferred. Instead, we have adopted 3 distinct analyses (the two Bayesian analyses and the non-Bayesian analysis using the EM algorithm) as a mechanism for assessing robustness. Broadly, the HMM parameter estimates and associated uncertainty intervals obtained from the three analyses should agree, and each should yield satisfactory predictive simulation results. The parameter estimates are examined in the next section. None of the posterior predictive simulation tests led to the suggestion that the Case 1 2-state HM models under consideration fail in their capacity to capture the observations. Similarly, a graphical inspection of the replicate time series generated from this model assessment exercise ( Figure 3) gives no indication of serious model inadequacy. The conclusion is that it is reasonable to use these models as a basis for characterising the observed seizure behaviour. The details are given in Note S5.

The posterior marginal sample obtained for each parameter of interest provides an estimate of its value (e.g., median and mode) together with a measure of uncertainty. In this paper uncertainty in the primary parameter estimates is summarised using 95% posterior intervals. (Narrower intervals are given for some secondary parameters.) Similarly, parametric simulation was used to generate 95% confidence intervals for the EM parameter estimates. Three sets (EM and two Bayesian model sets) of 2-state HMM primary parameter median estimates (i.e., transition probabilities and Poisson rates) are listed in Table 1, together with their 0.025 and 0.975 quantiles. Unsurprisingly, given only 53 observations, the posterior/confidence intervals obtained for the transition probabilities are large. The Poisson rate-parameter estimates ( λ ₁ and λ ₂) are similarly obtained with a relatively low precision but are, nevertheless, informative from the clinical perspective. The agreement between the three analyses is satisfactory, indicating a degree of robustness in the results. In particular, parametric simulation based on the EM solution is unrelated to MCMC parameter estimation based on a Bayesian model. Thus the level of agreement between the EM and MCMC results supports a claim of robustness.

An important feature of epilepsy seizure-count hidden Markov modelling is the potential to provide information regarding the seizure state of the subject on a given occasion. As stated in the Introduction, considerable insight into the underlying random process can be obtained using the EM algorithm with local and/or global decoding. Local decoding provides information on the state probabilities at each occasion while global decoding provides state-sequence probabilities. Related information is provided by the marginal posterior distribution obtained for the hidden states, Z _t , t = 1, …, N (i.e., the hidden seizure state at each occasion) given by the Bayesian analyses. In addition, mean duration statistics can be derived from the transition probability estimates.

Figure 4 shows the EM local decoding results, with superimposed error bars indicating the 0.025 and 0.975 quantiles given by the Poisson distribution with mean values taken from the EM solution. Hidden-state credibility time-series plots derived from the Poisson-gamma model marginal posterior samples are included in the figure. The agreement between the local decoding results and the Bayesian hidden-state credibility trajectory is striking, given the very different nature of the analyses on which these two results are based. The hidden-state mean-duration median estimates are included in Table 1, together with their 0.25 and 0.75 quantiles. These statistics are derived from the transition matrix diagonal elements, and their imprecision is, in part, a reflection of the lack of precision in the transition probabilities. The uncertainty in the mean duration estimate becomes inflated as the transition probability approaches unity. Nevertheless, the mean-duration 50% intervals are not so large as to render these estimates meaningless in the clinical context.

The main motivation for the HM modelling exercise outlined in this paper is the provision of information on the expected duration of each seizure state. Accordingly, the Case 1 2-state HMM transition-matrix posterior samples were used to calculate the sojourn-time probability densities, together with the associated cumulative distributions ( F( s), where s is sojourn time) and complementary cumulative distributions (1 −F( s)). The cumulative distributions (i.e., cumulative probabilities) provide an estimate of the probability that each seizure state will persist for up to some specified time, while each complementary cumulative distribution gives an estimate of the probability that a given state will persist for longer than the specified time. (As stated in the Methods section, this complementary cumulative probability is equivalent to the survivor function that is well known in survival analysis.) We have chosen to show the complementary cumulative distributions rather than the probability densities because they provide sojourn-time information in a form that is more easily communicated among clinicians. The Case 1 2-HMM log-linear model complementary cumulative probability distributions are shown in Figure 5.

In order to provide a clear clinical insight into the information that can be extracted from these distributions, and taking 2 weeks as an example, the median estimate obtained for the probability that each of the sojourn times exceeds 2 weeks (Prob( s > 2)) is 5.1% (25% and 75% posterior quantiles: 1.2%, 13%), and 68% (25% and 75% posterior quantiles: 56%, 78%) for States 1 and 2, respectively. Two points emerge from an inspection of Figure 5. Firstly, the estimation error is considerable. Secondly, regardless of the estimation error, a considerable range of sojourn times must be expected; this is an inescapable consequence of the stochastic nature of the processes and is unrelated to parameter estimation error. Periods of high seizure frequency lasting many weeks are not unexpected under these HM models. Thus the occurrence of a long period of higher than average seizure activity does not necessarily indicate a change in the underlying pathology or a permanent deterioration. An awareness that clinicians should expect to see a considerable variation in the duration of unusual seizure activity might minimise excessive treatment.

The most striking feature of the time series data acquired from the second patient is an abrupt fall in seizure rate after the 23rd week ( Figure 6), despite no change in treatment. The need to capture this change point, together with the occurrence of two occasions with an unusually high seizure rate, renders this second case somewhat more challenging than the first. Three alternative HMMs were examined in their capacity to capture the observations, namely a 3-state HMM with constant state-specific Poisson parameters (3-HM), a 3-state HMM with AR(1) time dependence in the Poisson parameters (3-HM-AR) and a 3-state HMM based on state-specific rate parameters that include a temporal random effect term with a conditional autoregressive prior distribution (3-HM-CAR). (The distinction between an autoregressive likelihood and an autoregressive prior is important.) Posterior predictive simulation was performed to examine the capacity of these models to capture important features of the data. Visual inspection of the plots given in Figure 6 indicates that the simple 3-state HMM appears to account for the main features of the data although the posterior predictive residuals show some structure. Adding AR(1) dependence offers no substantial improvement (graphical output produced using the AR1 model is not included in this report), while the 3-HM-CAR residuals show less departure from the ideal distribution about zero. Despite these observations, posterior predictive tests suggest that the simple 3-state models are adequate. Details are given in Note S6. Parameter estimates given by three 3-HM model analyses are listed in Table 2. An important finding emerges from the model assessment exercise, namely that the random effect term that is incorporated into the 3-HM-CAR model tends to dominate the fitted time dependence, hidden-state transitions making a relatively small contribution to the seizure rate changes. Thus, although the 3-HM-CAR model provides an important standard against which the simpler models can be compared, in Case 2 we regard the 3-HM-CAR model as uninteresting from a clinical perspective (see Note S6 for additional information).

The Case 2 hidden-state mean-duration statistics are listed in Table 2 together with the associated Poisson-rate parameter and transition probability statistics. The model parameter estimates together with the corresponding hidden-state trajectory plots (this hidden-state graphical output is not included in this report) indicate that transitions between adjacent states have a non-negligible probability both before and after the change-point and that the system is not held in a high-activity state prior to the changepoint and a low activity state afterwards.

Among the other points that emerge from the analysis is the considerable uncertainty in the hidden-state trajectory underlying the observed count data. The identity of the hidden state is particularly uncertain on some occasions. The fact remains that, under the Poisson distribution, a relatively high count can arise by chance even when the system is in a low seizure state, and cannot be attributed with certainty to an underlying high seizure-rate state. Thus, a series of consecutive occasions during which the seizure counts are consistently high does not rule out the possibility of underlying transitions between high and lower activity states. Furthermore, a sequence of consecutive occasions with no change in the most probable state does not preclude the possibility that transitions have occurred. Finally, Table 2 includes the 3-HM-AR Poisson rates and transition-probability estimates, showing that these differ substantially from those obtained with the simple 3-HM models. It is not surprising that introducing time dependence into the Poisson rate parameters substantially alters the transition statistics.

Additional analyses were performed using a number of 2-state and single-state models. Although the 2-HM-CAR model yields very satisfactory posterior predictive simulation results, it happens to exhibit the same random effect dominance described above for the 3-HM-CAR model. For that reason, in Case 2, it is uninteresting from a clinical perspective. Similarly, the single state-CAR model is dismissed as lacking in terms of the clinical insight it provides in Case 2, although it is not ruled out by the posterior predictive simulation analysis. The posterior predictive simulation results obtained with the 2-HM and 2-HM-AR models indicate that these models tend to produce Bayesian p-values that are low, although these do not fall below 0.05 (Bayesian p-value details not included). Based on these observations, the following section focusses on the sojourn time statistics derived from the time-independent log-linear 3-state model.

As in the previous case, we again resort to using the posterior marginal samples obtained for the transition matrix diagonal elements, in conjunction with the geometric distribution, as a mechanism for generating the state sojourn-time probability densities and their associated complementary cumulative distributions. Figure 7 shows the results given by the log-linear 3-HM model. Taking an example duration-of-interest equal to 2 weeks, the median estimate obtained for the probability that the sojourn time exceeds 2 weeks (Prob( s > 2)) is 46% (25% and 75% posterior quantiles: 12%, 71%), 25% (25% and 75% posterior quantiles: 6%, 50%) and 59% (25% and 75% posterior quantiles: 31%, 74%) for States 1, 2 and 3, respectively. Although the median curves might give the impression of a difference between States 2 and 3, the estimation error is too large to allow any definite conclusion other than the inescapable fact that, putting estimation error aside, the duration of each of the states is expected to exhibit considerable variability due to the stochastic nature of the process.

A striking feature of the time series collected from Case 3 is the appearance of a quadratic time dependence in the underlying seizure activity ( Figure 8). It is well known, however, that autoregressive processes can give rise to data that are deceptive in their shape, taking on a profile that appears deterministic. The preceding three 3-state HMMs (i.e., 3-HM, 3-HM-AR, and 3-HM-CAR models) were each investigated for their capacity to capture the observations, together with an HMM supplemented with a second-order polynomial (3-HM-quad). The corresponding 2-state models were also examined, together with the simple single-state quadratic model. Details concerning the various model adequacy tests that were performed are given in Note S7. Here we focus on the 2-state quadratic model (2-HM-quad), which was found to be adequate (Bayesian p-values of 0.21 and 0.34 were obtained for the two curvature tests defined in Note S7). The transition matrix and mean state-duration statistics derived from this model are included in Table 3.

Figure 9 shows the Case 3 sojourn-time complementary cumulative probability distributions given by the 2-HM-quad model. If, as in the previous cases, we take a 2-week duration-of-interest, the median estimates obtained for the probability that the sojourn time exceeds 2 weeks (Prob( s > 2)) are 25% (25% and 75% posterior quantiles: 14%, 39%) and 51.8% (25% and 75% posterior quantiles: 38%, 64%) for States 1 and 2, respectively. The considerable uncertainty in these estimates is a reflection of the limited amount of data available for analysis. This lack of precision must be taken into consideration if sojourn times are to be used in clinical decision making, albeit in conjunction with other information. Putting this imprecision aside, the results do serve to illustrate the important fact that, at any point in time, an additional and considerable uncertainty exists regarding imminent behaviour, due to the fundamentally stochastic nature of the process.

To keep this paper to a reasonable length we have not shown a comprehensive set of model assessment results, as obtained for every model and every case. Instead the content of the paper is designed to illustrate some of the modelling possibilities. The majority of posterior predictive simulation results are based on tests epileptologists will consider relevant to the clinical problem, for example, the number of occasions with an extreme seizure count. A related issue is whether or not the models should be regarded as definitive. We do not regard model selection with a view to finding the ‘best’ model as being among the objectives. In this context there is a need to distinguish between modelling the data and modelling the process that gives rise to the data. In the latter setting model comparison and the search for the best model has a greater relevance. From a clinical perspective interest might be focused on modelling the data while those engaged in basic science research may be more interested in modelling with a view to shedding light on the underlying pathophysiology. It is our view that a clinician whose concern is patient management will be primarily interested in empirical statistical models that provide a useful characterisation of the observed time series and an insight into future seizure activity. Given this objective, we do not seek to find the best model. Instead, a comparison of the results derived from differing models is undertaken as a mechanism for assessing robustness. If different models yield consistent parameter estimates this is taken as an indicator that the estimates are sound. In that sense we do not necessarily regard alternative models as competing. This is the reason for showing the results given by both the Poisson-gamma and log-linear Bayesian models where applicable. As shown, we obtained reasonable agreement between the two models. Given the differing form of the likelihood and priors under these two models, together with the Poisson parameter ordering constraint that is incorporated into the log-linear model, the similarity in the parameter estimates suggests that these are robust, i.e., the results exhibit a desirable insensitivity to the priors and are predominantly determined by the data. Similarly, in those cases in which the EM algorithm was adopted and combined with parametric simulation (an example is given in Table 1), agreement between the resulting confidence intervals and the posterior marginal distributions derived from the corresponding Bayesian analyses is an indication of robustness. Having taken this approach to assessing the reliability of the statistics obtained for a given subject, we have attached less importance to convergence testing and simulation accuracy calculations. For example, agreement between the MCMC results given by a Bayesian model and those produced using parametric simulation based on an EM solution would be unexpected if effective convergence to a stationary distribution had not been achieved in the former case. That said, MCMC convergence test results are reported, as obtained for key model parameters.