Robust heart rate variability analysis by generalized entropy minimization

Computational Statistics and Data Analysis 82 (2015) 137–151

Contents lists available at ScienceDirect

Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda

Robust heart rate variability analysis by generalized entropy minimization Davide La Vecchia a,b,∗ , Lorenzo Camponovo b , Davide Ferrari c a

Department of Econometrics and Business Statistics, Monash University, Australia

b

School of Economics and Political Science, University of St. Gallen, Switzerland

c

Department of Mathematics and Statistics, University of Melbourne, Australia

article

info

Article history: Received 24 February 2014 Received in revised form 18 August 2014 Accepted 3 September 2014 Available online 16 September 2014 Keywords: Heart rate variability Time series analysis Influence function M-estimators Robust estimation

abstract Typical heart rate variability (HRV) times series are cluttered with outliers generated by measurement errors, artifacts and ectopic beats. Robust estimation is an important tool in HRV analysis, since it allows clinicians to detect arrhythmia and other anomalous patterns by reducing the impact of outliers. A robust estimator for a flexible class of time series models is proposed and its empirical performance in the context of HRV data analysis is studied. The methodology entails the minimization of a pseudo-likelihood criterion function based on a generalized measure of information. The resulting estimating functions are typically re-descending, which enable reliable detection of anomalous HRV patterns and stable estimates in the presence of outliers. The infinitesimal robustness and the stability properties of the new method are illustrated through numerical simulations and two case studies from the Massachusetts Institute of Technology and Boston’s Beth Israel Hospital data, an important benchmark data set in HRV analysis. Crown Copyright © 2014 Published by Elsevier B.V. All rights reserved.

1. Introduction An electrocardiogram (ECG) is a representation of the heart electrical activity over a period of time, as detected by electrodes attached to the outer surface of the skin. The ECG is an important diagnostic tool for noninvasive investigation of heart rate and beat regularity, heart damages, effect of drugs or artificial devices on the cardiovascular system and many other medical applications. ECG signals are periodic waves characterized by the presence of sharp upward spikes, called R-peaks; the time between two consecutive R-peaks is called an RR interval (see Fig. 1). This paper is concerned with the analysis of HRV time series where an observation is given by HRV = 60/(RR interval length), measured in beats per minute. Over the years, HRV analysis has gained a major role in modern medicine for both diagnosis and modeling of the neurocardiovascular system; e.g., see Morillo (2012) for an introduction on this topic. Despite the importance of HRV analysis, however, relatively little has been written on robust statistical analysis of HRV times series. Relevant contributions in this area include Guo et al. (2001), Li (2008), Spangl and Dutter (2012), and Lu (2012). Typically, HRV times series are cluttered with anomalous observations, some of which have medical significance. Outliers can be divided into two main types: representative and non-representative outliers (Chambers, 1984). We regard representative outliers as atypical observations pertaining to the actual HRV process; they are correctly recorded values with medical meaning and there is no reason to exclude the chance of observing again similar values. Non-representative outliers are incorrectly recorded values or observations corresponding to unique events. In HRV analysis, measurement errors

∗

Corresponding author at: Department of Econometrics and Business Statistics, Monash University, Australia. E-mail address: [email protected] (D. La Vecchia).

http://dx.doi.org/10.1016/j.csda.2014.09.001 0167-9473/Crown Copyright © 2014 Published by Elsevier B.V. All rights reserved.

138

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

Fig. 1. RR intervals excerpt for Patient #100 from the MIT-BIH data set. Top panel: raw data (minus their mean) corresponding to 10 min from the central part of the ECG recording. Bottom panel: the same observations filtered by the nguess routine available on the MIT-BIH website.

(e.g., due to sudden movements of the patient) and artifacts (anomalous observations generated by the recording device) are non-representative outliers; ectopic beats (extra beats between two regular beats) can be either non-representative or representative outliers and can potentially convey information about meaningful but rare HRV patterns. In the context of time series models, such as the AR, ARMA and ARIMA models, an important example of representative outliers are the innovation outliers (IO), which affect the current observation as well as subsequent observations. Additive outliers (AO), defined as errors affecting a single observation, are examples of non-representative outliers. See Fox (1972) and Maronna et al. (2006, Chapter 8) for a more detailed discussion on IOs and AOs in the context of time series. Clinicians are well-aware of these types of outliers and usually attempt to deal with them by data preprocessing, often involving visual and ad hoc outlier detection techniques. For example, currently available R and MATLAB statistical software packages for HRV analysis mostly deal with outliers by allowing the user to manipulate directly the HRV records (Perakakis et al., 2010; Rodríguez-Liñares et al., 2011). This kind of pre-processing, however, comes with a major risk: the filtered trajectory might be clear of evident outliers, but other atypical observations incompatible with the assumed model remain hidden within the processed sample. For example, in Fig. 1 we show a sample for Patient #100 from the Massachusetts Institute of Technology and Boston’s Beth Israel Hospital (MIT-BIH) data set filtered by a standard routine.1 In the bottom panel, the filtered records between 400 and 700 (shaded area) have shifted location and large variability compared to the rest of the observations. In stationary time series, this anomalous behavior is typically produced by outliers (Maronna et al., 2006, Chapter 8). Autoregressive (AR), autoregressive moving average (ARMA), or bilinear models are standard model choices for HRV data (Christini et al., 1995; Spangl and Dutter, 2005). One of the most popular estimators for these models is the Gaussian pseudo maximum likelihood estimator (PMLE), where the Gaussian likelihood is used for the true but possibly unknown likelihood (Gourieroux et al., 1984). However, it is well-known that a handful of observations deviating from the parametric assumptions (e.g., observations violating the assumptions about either the first or the second conditional moment) can affect the PMLE, potentially leading to severely biased estimates (Mancini et al., 2005). Li (2008), Kaufmann et al. (2011), Clifford et al. (2012), and Behar et al. (2013) analyzed ECG records, showing that bias arises due to different types of anomalous records, such as measurement errors, artifacts, and ectopic beats. In this paper, we advocate the use of robust inference methods for HRV data. Robust methods have a clear advantage over non-robust alternatives in this context: they are capable of producing parameter estimates resistant to outliers; hence, they enable clinicians to detect anomalous physiological patterns with increased confidence. From this perspective, we study a new estimator in the context of time series models, which we refer to as tilted and conditional maximum Lq-likelihood estimator (TC-MLqE). Then, we use the TC-MLqE to analyze the MIT-BIH data set. The MIT-BIH data set has been studied extensively for assessing arrhythmia detectors by manufacturers, as well as for comparing basic heart-related research, and it represents one of the most important benchmarks in HRV analysis (Lin and Tian, 2012). 1 The MIT-BIH data set and the filtering routine nguess are freely available on http://www.physionet.org.

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

139

The TC-MLqE is obtained by minimizing a generalized pseudo-likelihood function. Such a criterion function can be viewed as a sample-based approximation of a generalized information measure often referred to as q-entropy (see Tsallis, 1988 and subsequent papers). The behavior of the resulting estimator is controlled by the constant q, which tunes the degree of robustness of the estimator. When the reference likelihood function is the Gaussian likelihood, the TC-MLqE defines a robust version of the PMLE: for q = 1 the TC-MLqE coincides with the PMLE, while, for q < 1, the estimator gains considerable robustness and is robust to evident or hidden outliers. Our estimation strategy is analogous to that first proposed by Windham (1995) in the i.i.d. setting to render robust estimating equations. The TC-MLqE boils down to correct the traditional likelihood scores by robust weights proportional to the assumed conditional density. A simple re-scaling transformation is applied to the parameter estimates to recover Fisher consistency. Ferrari and La Vecchia (2012) studied the robustness properties of this kind of estimator in the i.i.d. setting and pointed out its connection to a flexible family of generalized information measures. Note that in this paper we define a robust estimation procedure on raw data and avoid any additional data pre-processing or filtering methods. The proposed robust M-estimators have the same spirit as the Hampel-type infinitesimally robust M-estimators of Künsch (1984), which are resistant to local deviations from the reference model. If global robustness properties are desired (high BP), then additional filtering techniques may be applied before parameter estimation, such as the robust filtering procedures for ARMA processes proposed by Bianco et al. (1996) and Muler et al. (2009). See Spangl and Dutter (2012) for an application of similar methods in the context of HRV analysis. As it is customary in the literature of clinical signals, we study the empirical performance of our approach both in time and frequency domains (e.g., see Kannathal et al. (2007)). When applied to both simulated and real HRV time series, the TC-MLqE shows remarkable stability compared to other classic estimators. Model fitting by the TC-MLqE allows for sharp detection of representative and non-representative outliers in both time and frequency domains. We believe that this property makes our approach a valid diagnostic tool in HRV analysis, which helps clinicians identify anomalous HRV patterns. Differently from other traditional robust estimation methods, such as those involving Huber-type weights (see, e.g., La Vecchia and Trojani, 2010), the TC-MLqE is very easy to implement, since it avoids the usual numerical integration required to preserve Fisher consistency. The remainder of the paper has the following structure. In Section 2, we introduce the model setup, some basic concepts in robust inference for time series, and describe our robust estimation approach. In Section 3, we address robust model selection and propose a data-driven method to select the tuning constant q. In Section 4, we illustrate the properties of our method by several Monte Carlo simulation studies. In Section 5, we apply our method to two case studies within the MIT-BIH data set. Section 6 concludes the paper with final remarks. 2. Methodology Let Y := (Yt )t ∈Z be a real-valued strictly stationary and ergodic time series, defined on the probability space (R∞ , F , P∗ ). Let P := {Pθ , θ ∈ 2 ⊂ Rp } be a parametric family of models for P∗ . The process dynamics is represented by the following location-scale model: Yt = m(Yt −1 , Xt , εt −1 , ϑ) + σ εt ,

0 ≤ t ≤ n,

(2.1)

where εt is the error (shock) at time t, with zero mean and unit variance, θ = (ϑ , σ ) is the location/scale parameter vector. The function m(·) is differentiable in θ and, as in Drost et al. (1997), it parameterizes the conditional mean of the process. The mean function m depends on the past values of the process, Yt −1 , on the past shocks, εt −1 , and on a vector of exogenous covariates, Xt , encompassing age, gender, and other health indicators. For clarity of presentation, we drop the dependence on Xt and εt −1 , but our results can be easily extended to the general case. s v Expression (2.1) includes important time series models for different choices of m. For instance, i=1 ϑi Yt −i + l=1 βl Xt −l T

+

k

j =1

T

ϕj εt −j gives an ARMAX(s, v, k) process, where the AR and ARMA models are special cases. Nonlinear processes such s

as the exponential AR(s) of Haggan and Ozaki (1981) are obtained by setting i=1 (υi + ϑi exp(−cYt2−d ))Yt −i , where c ≥ 0, d ∈ N and υi , ϑi ∈ R. Threshold AR and self-exciting threshold AR models are also special cases of (2.1). 2.1. Robust inference: M-functionals and M-estimators From a robust viewpoint, model (2.1) is interpreted as an approximate description of the actual data generating process, P∗ . Our aim is to define computationally feasible estimators for the parameters in (2.1), when P∗ is in a neighborhood of a reference model Pθ 0 . To formalize this idea, let Yts := (Yt , . . . , Yt −s ), s ∈ N, denote the finite random sequence of Y and Pθs 0 (or P∗s ) the s-dimensional marginal distribution of Yts induced by Pθ 0 (or P∗ ). We allow the reference model Pθs to be misspec0 ified for the true model P∗s , by considering local deviations from Pθs . In particular, we assume that P∗s is in the nonparametric 0 s s neighborhood Uη (Pθs ) of Pθs , defined as Uη (Pθs ) = {Pηs = (1 − η)Pθs + ηQ s , η ≤ b, b ∈ [0, 1], Q s ∈ Mstat }, where Mstat de0 0 0 0 η s notes the class of s-dimensional marginals of strictly stationary and ergodic processes. The neighborhood U (Pθ ) contains 0 all the distributions Pηs with Kolmogorov distance from the reference model Pθs bounded by η. A number of standard con0 tamination models, such as isolated, replacement, or additive outliers, are well-approximated by distributions in Uη (Pθs ) 0 for small η; see Künsch (1984) for more details.

140

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

To specify a robust estimator for θ 0 when local departures from Pθs are allowed, e.g., due to data contamination, artifacts 0

and/or recording device failures, we work within the framework of M-estimation theory. Let ψ : Rs+1 × Rp → Rp be a martingale difference such that Eθ 0 [ψ(Yt ; Yts−1 , θ 0 )|Ft −1 ] = 0, for all t = 1, . . . , n, where Eθ 0 [·|Ft −1 ] denotes the expectation s → 2 is defined by the unique root under Pθ 0 conditional to the σ -algebra up to time t − 1. The M-functional T (·) : Mstat of the system of moment conditions Eθ 0 [ψ(Yt ; Yts−1 , T (Pθ 0 ))] = 0 and is such that T (Pθ 0 ) = θ 0 . Given a sample, Y0 , . . . , Yn ,

the M-estimator θˆ n = T (Pn ), with Pn representing the empirical measure, is defined by the zero of the following equations: n 

ψ(Yt ; Yts−1 , θˆ n ) = 0.

(2.2)

t =1

√

s +1 This is equivalent to maximizing Hn (θ) = , θ), where ρ satisfies ψ = ∇θ ρ . Under standard conditions, n(θˆ n − t =1 ρ(Yt −1 θ 0 ) converges in distribution to N (0, J (θ 0 ) K (θ 0 )J (θ 0 )−1 ), where J and K are p × p matrices defined by J (θ) = −Eθ [∇θ ψ (Yt ; Yts−1 , θ)] and K (θ) = Eθ [ψ(Yt ; Yts−1 , θ)ψT (Yt ; Yts−1 , θ)]. Assume that the conditional transition density function for Y is specified by the reference parametric model y → pθ (y; u), y ∈ R, and u ∈ Rs . For estimating   model (2.1), typically the pseudo maximum likelihood estimator (PMLE) is computed assuming Yt |Ft −1 ∼ N (m Yts−1 , ϑ , σ ), so that pθ is the Gaussian density function. Thus, Gaussian PMLE is obtained when using the martingale difference

n

ψPML (Yt ;

Yts−1

, θ) = ∇θ log

1

σ

 φ

Yt − m Yts−1 , ϑ





σ

in (2.2), where φ(·) denotes the standard Gaussian probability density function. A desirable property for an estimator is that minor deviations (due to the presence of outliers) of the data distribution from the ideal Gaussian model imply only minor changes of θˆ n . However, this requirement is typically not satisfied by the Gaussian PMLE, which is known for being extremely sensitive to outliers and, hence, severely biased (Mancini et al., 2005). To see this, it is useful to introduce the notion of conditional influence function. Definition. The conditional influence function (IFc ) for the M-functional T (·) is defined by the function IF c : Rs+1 ×2 → Rp such that:

∂ T (Pηs )|η=0 = ∂η  c

 Rs+1

IF c (y, u; θ 0 )Q s (dy, du)

(2.3)

and satisfying Eθ 0 IF (Yt , Yts−1 ; θ 0 )|Ft −1 = 0.



For an M-functional implied by a martingale difference ψ , La Vecchia and Trojani (2010) show that IFc exists and it is unique for ergodic and stationary Markov processes. From Corollary 5 in La Vecchia and Trojani (2010) it follows that for the class of models in (2.1) IF c (y, u; θ 0 ) = J (θ 0 )−1 ψ (y; u, θ 0 ) .

(2.4)

In first-order infinitesimal robustness, the main interest is to obtain estimators with uniformly bounded asymptotic bias, where the bias is characterized by the linear approximation B(η) = η∂η T (Pηs )|η=0 . From (2.3) we argue that a bounded IFc ensures a bounded linearized asymptotic bias (Hampel et al., 1986). However, since ψ PML is unbounded over R, the Gaussian PMLE is not infinitesimally robust. This aspect explains the sensitivity of the Gaussian PMLE even to mild deviations from the model assumptions. 2.2. Tilted conditional maximum Lq-Likelihood estimator ∗

For a reference density pθ (Yt ; Yts−1 ), we consider the estimator θˆ n , which minimizes the surrogate objective Hn∗ (θ) = −

n 1

n t =1

logq [pθ (Yt ; Yts−1 )],

q ≤ 1,

(2.5)

over 2, where logq (x) = (x1−q − 1)/(1 − q), if q ̸= 1, and logq (x) = log(x), if q = 1. This is equivalent to solving n 1

n i =1

ψq (Yt ; Yts−1 , θ) = 0 with ψq (Yt ; Yts−1 , θ) = ∇θ log pθ (Yt ; Yts−1 )pθ (Yt ; Yts−1 )1−q .

(2.6)

The expression in (2.5) can be interpreted as a sample-based version of the conditional q-entropy function (Tsallis, 1988; Ferrari and La Vecchia, 2012). Under correct model specification, Hn∗ (θ 0 ) → H ∗ (θ 0 ) = Eθ 0 {logq [pθ 0 (Yt , Yts−1 )]|Ft −1 }, as

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

141

Table 1 Re-centered parameters by the transformation τ q for common examples of parametric densities. pθ

θ

τ q (θ)

Normal (mean, st.dev.) Non-standardized t-student (mean, scale, df) Exponential (rate) Laplace (mean, scale)

µ, σ µ, σ , ν β µ, σ

√ µ, σ / q µ, σ [(ν + 1)q − 1]ν −1 , (ν + 1)q − 1 β/q µ, σ /q

n → ∞, where H ∗ represents the q-entropy measure, conditional to Ft −1 . The above estimator is generally biased and corresponds to the maximum Lq-likelihood method of Ferrari and Yang (2010) within the i.i.d. framework. In model (2.1) with ∗

correct specification of the transition density, θˆ n is consistent for ϑ , but not for the scale parameter σ . Hence, the definition of a Fisher consistent M-estimator requires an appropriate scale transformation τ q , which we examine next. 2.2.1. Fisher consistency by tilting When q is a fixed constant, asymptotic normality of the M-estimator defined by ψ q follows from standard results (Bollerslev and Wooldridge, 1992). However, the definition of a Fisher consistent M-estimator requires additional care. The approach followed here is motivated by the need to avoid the usual re-centering strategies in robust M-estimation, since they require involved (and often computationally intractable) multi-dimensional numerical integration. q Let τ q : Θ → Θ be a bijective function such that pτ q (θ) (y; u) ∝ pθ (y; u), for all pairs (y, u). A consistent estimator is

∗ ∗ θˆ n = τ q (θˆ n ), where θˆ n is the solution of (2.6). Note that this is equivalent to minimize the re-centered empirical q-entropy 1 ∗ Hn (θ) = Hn∗ (τ − q (θ)), where H is the surrogate objective function (2.5). The above re-scaled estimator is analogous to the

robust approach proposed by Windham (1995) in the i.i.d. setting, where weights proportional to the assumed model are introduced to render the log-likelihood scores robust and model parameters are rescaled to obtain consistency. Here, we adapt both Windham’s weighting scheme and parameter re-scaling to the time series setting of model (2.1), where the aim is to make the conditional pseudo-likelihood function robust. The final estimator θˆ n is Fisher consistent. This can be seen by taking the expectation of (2.7) under model pθ 0 :



pθ 0 (y; yts−1 )ψq (y; yts−1 , θ)dy =



pθ 0 (y; yts−1 )

∇θ pθ (y; yts−1 ) dy. q pθ (y; yts−1 )

1 Then, evaluating the above expression in θ ∗ = τ − q (θ 0 ) we obtain



∇θ∗ pθ∗ (y; yts−1 )dy = ∇θ∗



pθ ∗ (y; yts−1 )dy = ∇θ ∗ 1 = 0.

The above derivation shows that when the true model is a member of P , the TC-MLqE is consistent for the true parameter θ 0 . For common families of error distributions – such as Gaussian, Laplace, exponential, or t-distributions – the transformation τ q has a closed form, regardless of whether or not m(·) is linear; Table 1 shows the function τ q for different error densities. Hence, a major advantage of the TC-MLqE compared to Huber-type M-estimators, such as the approaches in Künsch (1984), Mancini et al. (2005), or La Vecchia and Trojani (2010), is that the function preserving the Fisher consistency admits an easy-to-implement analytical expression, which avoids the difficulties related to numerical integration. For example, the M-estimator in Mancini et al. (2005) uses Laplace method to preserve Fisher consistency when m(·) is nonlinear. We believe that closed-form re-centering is an important aspect for the practitioner, since it simplifies considerably the implementation of the estimator and reduces its computational burden. We illustrate this point in the following example. Example 1. When εt ∼ N (0, 1) in (2.1), ψPML defines a root-n consistent M-estimator for (ϑ T , σ ) and achieves full parametric efficiency. A consistent Gaussian TC-MLqE for both location and scale is obtained by solving n 1

n i=1

 ψPML (Yt ;

Yts−1

, (ϑ, σ ))φ

Yt − m Yts−1 , ϑ



 1−q

σ

= 0,

(2.7)

and then applying the transformation τ q (·) to the solution. For the mean and scale parameters we have τ q (ϑ) = ϑ , and √ τ q (σ ) = σ / q, respectively. 2.2.2. Re-descending conditional influence function 1 −q

The IFc of the TC-MLqE is proportional to ψ q = pθ ∇θ log(pθ ), where ∇θ log(pθ ) is the likelihood score. The usual maximum likelihood estimator (MLE) is obtained when q → 1, while q < 1 defines robust M-estimators since, typically, the 1 −q term pθ goes quickly to zero as |y| → ∞ and dominates ∇θ log(pθ ). Therefore, the IFc is re-descending and the extreme observations cannot influence much the estimates. This implies that our estimation approach rejects completely gross outliers and mitigates the negative impact of anomalous records more effectively than the Huber-type M-estimators, which do

142

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

not completely down-weight large observations. Although the influence function is re-descending for many useful models, such as for the models in the exponential families (e.g., for Gaussian reference density), this property should be verified on a case-by-case basis; see Ferrari and La Vecchia (2012) for a related discussion. The breakdown point for the TC-MLqE is expected to decrease as q → 1, achieving its maximum when q is far away from 1. Thus, one should avoid choosing q too close to 1 in the presence of a relatively large number of outliers. This understanding is confirmed by exploratory analyses on simulated data. For example, simulations for the ARMA(1, 1) model within the setup given in Section 4.2 (sample size n = 1000) showed that for q between 0.75 and 0.95 the mean squared error grows exponentially fast when the fraction of contaminated data is larger than 20%. When q = 0.99 the mean squared error grows rapidly when the fraction of contaminated data is larger than 15%. A data-driven strategy for choosing q which reflects the above considerations is discussed in the next section. A theoretical breakdown analysis of our estimator is not pursued in this paper; however formal results may be derived by a methodology analogous to that of Basu et al. (1998). 3. Computational aspects 3.1. Selection of the tuning parameter q The selection of the degree of robustness for the TC-MLqE is an important aspect in applications. In robust statistics, there are two popular approaches for selection of analogous tuning constants. A first option is to select the constant that achieves a target asymptotic relative efficiency of the robust estimator compared to the maximum likelihood estimator, at the reference model (Mancini et al., 2005). A second approach is to control the maximal bias of the test based on the robust estimating function, assuming a theoretical level of contamination (Ronchetti and Trojani, 2001). Although both approaches result in robust estimators, by construction they are sub-optimal when the model is correctly specified. An additional drawback of asymptotic approaches is that they rely on the prior specification of the (theoretical) contamination level, which is hard to determine. We propose a data-driven method for the selection of q, which ensures full efficiency in absence of contamination. Our method relies on the empirical stability of the estimates rather than asymptotic criteria. The main idea is to select the tuning constant q closest to 1 (i.e., closest to MLE) such that the point estimates of the parameters of interest are sufficiently stable. In a sample without anomalous observations, q = 1 would suffice to obtain stable estimates, and such a choice would also give the same efficiency as MLE. In the presence of outliers, values of q close to 1 give unstable estimates, so we move q away from 1 until sufficient stability is recovered. This rationale suggests the following steps. Algorithm 1 (Selection of q). (1) Define an equally spaced grid 1 ≥ q1 > q2 > · · · > qn ≥ 0 and compute the correspondent point estimates, θˆ qi , i = 1, . . . , n. (2) Compute the absolute variations AVqi = ∥θˆ qi − θˆ qi+1 ∥. (3) Select the optimal value q∗ = {max qi : AVqj < ϱ, for all qj ≤ qi }, where ϱ > 0 is some threshold value. By definition, q is the tuning constant in the grid 1 ≥ q1 > q2 > · · · > qn ≥ 0 closest to 1 such that the variation of the point estimates is smaller than an acceptable value ϱ. From (2.4) and (2.7) we note that for extreme observations the influence function of the TC-MLqE decreases as q decreases. Therefore, by selecting qn close to 0 we can tune the gross error sensitivity (defined as the supremum of the conditional influence function) of our estimator. Based on our simulations, we recommend a grid between q1 = 1 and qn = 0.9, with step 0.01, and ϱ = 5% · ∥θˆqn ∥. Our numerical results suggest that choices around 0.9 typically yield considerable robustness. If only a small percentage of the data are contaminated, it may be useful to consider progressive refinements for the step size towards 1 to improve the efficiency of the estimator (for example, qi = 1 − (i/n)2 , i = 1, . . . , n). In the presence of a large amount of contaminated data, our selection of q might yield a biased M-estimator so other procedures such as the most bias-robust estimator (MBRE) could be investigated (see, e.g., Hampel et al., 1986). Example 2. We generate 5000 samples (Y1 , . . . , Y600 ) from a AR(3) model with N (0, 1) errors and parameters: ϑ1 = 0.5, ϑ2 = −0.6, ϑ3 = 0.75, and σ = 1. Outliers are introduced by replacing Ymax by 5 · Ymax in each sample, where Ymax = max(Y1 , . . . , Y600 ). In Fig. 2, we plot values of q ranging from 0.9 to 1 against parameter estimates. Without outliers, the point estimates are nearly constant for all values of q. In this case, q = 1 gives stable estimates compared to other choices of q, so our procedure selects q∗ = 1. In the presence of contamination, the point estimates are very unstable as q approaches 1. As q moves away from 1, the point estimates become increasingly stable. When q < 0.96, the estimates converge to their true values and our procedure selects q∗ = 0.95. 3.2. Model selection A popular method for model selection is Akaike’s Information Criterion (AIC), which relies on the likelihood function. However, in the presence of contamination, AIC is typically not trustworthy due to the non-robust nature of the likelihood score. Many authors, including Ronchetti (1997), Shi and Tsai (1998), have highlighted these negative features of AIC and

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

143

Fig. 2. Selection of q for the AR(3) model with and without outliers. In each plot, the central dotted line is the median of the Monte Carlo estimates, the dashed lines are the maximum and minimum Monte Carlo estimates, the horizontal solid line is the true parameter value. The plots correspond to non-contaminated and contaminated samples. Results are based on 5000 Monte Carlo samples of size n = 600.

stressed the importance of robust model selection procedures. Following Ronchetti (1997), we select the model which minimizes





AICR = Hn (θˆ n ) + tr J (θˆ n )−1 K (θˆ n ) .

(3.1)

1 ∗ For the TC-MLqE, we have Hn = Hn∗ (τ − q (θ)), where Hn is the empirical q-entropy function defined in (2.5), and J and K are the matrices defined in Section 2.1.

4. Monte Carlo simulations 4.1. AR process We illustrate the stability properties of our procedure through a Monte Carlo study where we compare results obtained from clean and contaminated processes. The clean process is the AR(3) of Example 2. As in Künsch (1984), the contaminated process is obtained by the replacement model Y (t ) = Htϵ X (t ) + (1 − Htϵ )ξ (t ), where: X (t ) is the clean AR(3) process; Htϵ is a binary ergodic, stationary 0–1 process; and η = P (Htϵ = 1) is the contamination level. In our Monte Carlo setting, the contamination scheme generates different types of outliers: (i) Innovation outliers (IO)—observations generated by the reference AR(3), but with larger error variance than that for the clean process; (ii) Replacement outliers (RO)—white noise observations with up-shifted mean and a larger variance; and (iii) Additive outliers (AO)—observations having the AR(3) dynamics summed to another stationary process, such that the resulting observations have down-shifted means and slightly larger variances. More details on these outlier types are described in Maronna et al. (2006). Fig. 4 (top panel) shows a sample generated by our simulation setting. Note that the resulting simulated dynamics is similar to that of Patient #100 in the MIT-BIH data set (Fig. 7, top). In our simulations, we compare the following M-estimators: Gaussian PMLE (PMLE), Gaussian TC-MLqE (TC-MLqE), and Huberized Gaussian PMLE (H-PMLE). Prompted by an anonymous referee, we also compared our approach with the high breakdown-point and high efficiency robust MM-estimator by Salibian-Barrera and Van Aelst (2008) (with Tukey’s bi-weight objective function and scale estimated by re-scaled MAD of the residuals) and the fast robust S-estimator by Salibian-Barrera and Yohai (2006). The TC-MLqE is implemented as in Example 1, with Gaussian transition density and linear mean function m. The H-PMLE obtains robustness by imposing Huber-type weights on the PMLE errors. The robust H-PMLE weights take the form wb (x) = min{1, b/∥x/σˆ R ∥}, where ∥ · ∥ is the Euclidean norm and σˆ R is the median absolute deviation (MAD) from the median. 2 ˆ Time domain. For each Monte Carlo run, we compute the squared error ∥θˆ −θ∥2 = j=1 (ϑj −ϑj ) for the three M-estimators under clean and contaminated samples. We consider 5000 Monte Carlo samples of size n = 800 and n = 2400 and include 52 outliers (this corresponds to contamination levels η = 6.5% and 2.2%, respectively). To obtain a fair comparison, all the tuning parameters for all the robust methods are calibrated so that they all achieve nearly identical Monte Carlo mean squared error at the model (η = 0) with 95% efficiency compared to PMLE. The boxplots in Fig. 3, show the distribution of the Monte Carlo mean squared error under contaminated models (η > 0). Not surprisingly, the PMLE does poorly compared to the other estimators, while the performance of the TC-MLqE is comparable to that of the high-breakdown estimators.

3

144

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

0

0

-2 -2 -4 -4 -6 -6 -8

-8

-10

Fig. 3. Distribution of the Monte Carlo errors ∥θˆ − θ∥2 (on a log scale) for contaminated and clean AR(3) processes. The model parameters are estimated by PMLE (1), H-PMLE (2), TC-MLqE (3), MM-estimator with Tukey’s biweight objective function and scale estimated by re-scaled MAD of the residuals (4), and Fast S-estimator (5). The clean process is AR(3), εt ∼ N (0, 1), ϑ1 = 0.5, ϑ2 = −0.6, ϑ3 = 0.75, and σ = 1. The contaminated process contains 52 observations representing different patterns of patchy outliers (IO, AO, RO). The tuning constants for the robust methods are calibrated so that the robust methods give a 5% efficiency loss with respect to PMLE under the clean process. Plot based on 5000 Monte Carlo samples.

Frequency domain. We compare the spectral density estimates for clean and contaminated AR(3) models. The parameters of the underlying process are estimated using PMLE, H-PMLE, TC-MLqE, MM- and Fast S approaches described above. Fig. 4 shows the estimated spectral densities and their 95% confidence bands, for PMLE, H-PMLE and TC-MLqE. The plots for the MM and Fast S estimators are similar to those obtained by TC-MLqE, so they are not to be presented here. Under the clear model, the robust estimators are calibrated in time domain to achieve a 5% small efficiency loss with respect to the PMLE, so their estimated spectral densities are similar. In the presence of contamination, the PMLE gives an almost flat – and therefore strongly biased – estimated spectral density, while H-PMLE and TC-MLqE clearly detect the spectral peak around 0.2 Hz. The median estimated spectral density for the TC-MLqE is closer to the target reference model than the standard H-PMLE robust method. 4.2. ARMA process We generate data from the ARMA(1, 1) model Yt = ϑ Yt −1 + ϕϵt −1 + ϵt , θ = (ϑ, ϕ). We consider sample sizes n = 1000 and n = 2000 and simulate both clean and contaminated samples. In the contaminated samples we include 32 and 132 outliers, corresponding to 3% and 6% contamination levels, respectively. The clean samples are generated using ϵt ∼ N (0, 1), while the contaminated samples are generated using the same kind of replacement model as described in Section 4.1. We compare the performance, in time and frequency domains of Gaussian PMLE, H-PMLE, and TC-MLqE. As in Section 4.1, we compared our approach with the high breakdown-point and high efficiency robust MM-estimator by Salibian-Barrera and Van Aelst (2008) (with Tukey’s biweight objective function and scale estimated by re-scaled MAD of the residuals) and the fast robust S-estimator by Salibian-Barrera and Yohai (2006). Time domain. We consider 5000 Monte Carlo runs and compute the squared error ∥θˆ − θ∥2 = (ϑˆ − θ )2 + (ϕ − ϕ) ˆ 2 for each run. To obtain a fair comparison, the tuning parameters for all the robust methods are calibrated so that the corresponding estimators achieve nearly identical Monte Carlo mean squared error at the model (η = 0) with 90% efficiency compared to PMLE. In Fig. 5, we show the resulting box plots. In the presence of outliers, the Gaussian PMLE gives large errors for both the AR and MA parameters and is outperformed by all the robust methods. TC-MLqE appears to be more stable than H-PMLE in the presence of contamination and the distribution of the TC-MLqE errors is close to zero with a relatively small dispersion. In the considered setting, our approach performs comparably to the MM- and fast S-approaches, with all the robust methods outperforming H-PMLE. Frequency domain. Let f (ω, θ) be the ARMA(1, 1) spectral density (see, e.g., Maronna et al., 2006 for the analytical expression). We evaluate the impact of outliers in the frequency domain by evaluating the integral of f (ω, θ) at high frequencies;  0.5 particularly, as an illustration, we consider HF(θ) = 0.45 f (ω, θ)dω. Similar spectral measures related to the integral of the spectrum over a given set of frequencies are routinely applied in HRV analysis (e.g., see Malik et al. (1996)). Fig. 6 shows violin plots (boxplots with density estimates) representing Monte Carlo estimates for HF( θ n ), where the estimates,  θ n , are obtained from samples of size n = 2000 with and without outliers. The value of HF evaluated at the true parameter is represented by the horizontal line. All the methods show similar HF values under correct model specification, but the two robust Mestimators have a slightly larger variability, due to some efficiency loss compared to PMLE. In the presence of contamination, the PMLE is clearly unreliable, since it shows bias and large variability. In contrast, the two robust M-estimators yield stable HF estimates, with the TC-MLqE being closest to the true value of HF. Analogous plots (not reported here) and conclusions are obtained when the sample size n = 1000 is considered. The plots for the MM- and fast S-estimators resemble those obtained for the TC-MLqE, so they are not presented here.

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

145

Fig. 4. Monte Carlo simulations for n = 800. Top panel: Simulated clear (dotted) and contaminated (solid line) trajectories. The contaminated trajectory contains 52 outliers. Bottom panel: Monte Carlo spectral density estimates for clear and contaminated processes fitted by Gaussian: PMLE, H-PMLE and TC-MLqE. In each plot, the solid line represents the spectral density under the uncontaminated model, the dotted lines represent the Monte Carlo median, the dashed lines are 95% Monte Carlo confidence bands. Plots based on 5000 Monte Carlo runs.

5. Case studies 5.1. Arrhythmia data In our first example, we analyze the MIT-BIH arrhythmia data set consisting of a collection of ECG signals from patients with different heart conditions. The data are freely available on http://physionet.org/physiobank/database/mitdb/. We focus on Patient #100 (male, 69 years old), a benchmark subject for which a number of analyses in the frequency domain are available in the literature; see Goldberger et al. (2000) and subsequent related papers.

146

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

Fig. 5. Distribution of the Monte Carlo errors ∥θˆ − θ∥2 for contaminated and clean ARMA(1,1) process (on a log scale). The model parameters are estimated by PMLE (1), H-PMLE (2), TC-MLqE (3), MM-estimator (4), and S-estimator (5). The clean process is a ARMA(1, 1) with parameters ϑ = 0.2 and ϕ = 0.75 and error N (0, 1). The tuning constants for the robust methods are calibrated so that the robust methods give a 10% efficiency loss with respect to PMLE under the clean process. Plot based on 5000 Monte Carlo samples.

Fig. 6. Monte Carlo estimates of HF( θ n ) for the ARMA(1, 1) process with parameters ϑ = 0.2 and ϕ = 0.75 for clear and contaminated samples, having size n = 2000. The violin plots (boxplots with density estimates) correspond to Gaussian PMLE, H-PMLE, and Gaussian TC-MLqE. The horizontal line represents the true HF. The tuning constants for the robust methods are calibrated so that they give a 10% efficiency loss with respect to PMLE under the clean process.

Our sample contains 3600 records, corresponding to approximately 15 min of ECG signal. Previous work agrees that 0.33 Hz is the dominant frequency for Patient #100. This frequency is important because of its relationship with respiratory sinus arrhythmia (RSA), a type of arrhythmia related to the HRV dynamics during a breathing cycle. A list of artifacts and the corresponding frequencies for Patient #100 is available at the MIT-BIH website, as well as other information related to atrial/ventricular ectopic beats.

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

147

Fig. 7. Patient #100 data. The top panel shows the demeaned HRV series, while the bottom panel shows normalized weights w ˆ i∗ = w ˆ i / maxi w ˆ i , based on TC-MLqE. In both plots, the horizontal axis represents the observation number.

a

b

Fig. 8. First 10 s of ECG signal for Patient #100. Panel (a): first 10 s of raw signal with an ectopic event emphasized by the shaded area (short RR interval followed by long RR interval). Panel (b): 10 s for the HRV containing the atrial ectopic event (top) and corresponding weights w ˆ i∗ (bottom).

Time domain: HRV anomalies and ectopic beats. In this section, we illustrate our procedure by analyzing the HRV for Patient #100. Following Kannathal et al. (2007), we consider a simple AR(s) process with Gaussian errors. Model parameters are estimated using PMLE, H-PMLE and TC-MLqE. For the TC-MLqE, we selected q by the procedure described in Section 3, obtaining q = 0.95. The AR order, s, for PMLE is selected using the standard AIC criterion, resulting in s = 50. For H-PMLE and TCMLqE, we applied the robust AICR selection method obtaining a much more parsimonious model with s = 25. A choice of the autoregressive order near s = 20 is deemed as satisfactory from a prediction viewpoint and give good performances in applications (Kannathal et al., 2007). A more aggressive selection of the autoregressive order may be obtained using the robust BIC criterion for M-estimators proposed by Machado (1993). The simple AR(s) model is a standard choice in the HRV literature and allows us to compare our method. However, more sophisticated applications of our method to other models, including ARMA or ARMAX models, may be explored to improve the trade-off between empirical performance and model complexity. Fig. 7 shows the (demeaned) HRV signal and the estimated weights, respectively. The weights w ˆ i = p(y; u, θˆ n )1−q , i = 1, . . . , n, detect the presence of anomalous observations incompatible with the AR model. For clarity of presentation, we show the rescaled weights w ˆ i∗ = w ˆ i / maxi w ˆ i : weight close to one indicates that the observation is in agreement with the model assumptions, while weight close to zero implies incompatibility with the Gaussian AR model. Note that the observations between 1500 and 2000 receive small weights. These records coincide with values between 400 and 700 in Fig. 1, which we tagged as worrisome in Section 1 due to their anomalous location and scale compared to the rest of the data. Interestingly, our robust procedure gives almost zero weight to initial observations, which indicates an atypical heart activity. A closer inspection of the first ten seconds of HRV reveals the presence of an anomalous pattern, which can be clearly seen in Fig. 8. In Fig. 8, Plot (a) and Plot (b) (top) show the occurrence of an atrial ectopic event in terms of the raw and HRV series. An atrial ectopic event is characterized by one short RR interval between two regular RR intervals in the ECG signal (shaded area in left plot) and results in a sudden change of the HRV (shaded area in the right plots). The bottom plot of Fig. 8(b) shows

148

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

Fig. 9. Spectral density estimates for Patient #100. The plots correspond to Gaussian PMLE, H-PMLE with b = 1.345 and Gaussian TC-MLqE with q = 0.05. The location of artifacts related to recording devices is given by the dashed vertical lines. The frequency for the respiratory sinus arrhythmia (RSA) is highlighted by the vertical solid line.

that TC-MLqE detects this anomalous pattern by assigning very near-zero weights to the observations related to the cardiac ectopic event. After the ectopic beats, the weights move away from zero, meaning that a typical cardiac activity is restored. Frequency domain: arrhythmia and artifacts. Fig. 9 shows the AR spectral density estimated using Gaussian PMLE, H-PMLE, and TC-MLqE. In the same figure, we also show the frequencies (vertical dotted lines) of the artifacts described by Goldberger et al. (2000). The spectral density of PMLE detects RSA, but it also shows three spurious peaks related to different artifacts. The spectral density from H-PMLE detects the RSA. Although the H-PMLE is not influenced much by low-frequency artifacts, it is still sensitive to other high-frequency artifacts. Differently from the other two M-estimators, the TC-MLqE yields a spectral density which clearly detects RSA and, at the same time, it is not much influenced by the artifacts at any other frequency. 5.2. Congestive heart failure data In our second example, we study Patient #chf01 (male, 71 year old) from the MIT-BIH congestive heart failure data set, which includes long-term ECG recordings from subjects with severe congestive heart failure. The individual recordings are each about 20 h in duration. We present here results based on a central excerpt of 1950 HRV observations, corresponding to approximately nine minutes of ECG records. The results presented here do not change sensibly when considering different time intervals, for example at the beginning or at the end of the sample. Time and frequency domain analysis. We compare three estimators: PMLE, H-PMLE, and TC-MLqE. Following previous work in the literature, we fix the number of lags as s = 14 for all the considered estimators (Kannathal et al., 2007) in order to facilitate comparisons. We estimate the parameters in the time domain, and then use them to construct spectral densities. Fig. 10 (left panels) shows the demeaned HRV and the weights for the TC-MLqE. For TC-MLqE, we selected q = 0.95 based on algorithm presented in Section 3.1. Fig. 11 shows spectral densities for PMLE, H-PMLE, and TC-MLqE. Although the three spectra detect the RSA around 0.3 Hz, the density implied by the PMLE partially ignores the low frequencies from 0.01 Hz to 0.2 Hz, while both the robust M-estimators detect relevant heart activity at low frequencies. The results from TC-MLqE confirm the findings by Goldberger et al. (1988), who claim that two of the main aspects of the HRV dynamics for this patient are: (i) sustained low frequency oscillations in heart rate; (ii) abrupt HRV changes. The spectrum implied by PMLE simply fails to capture aspect (i), while TC-MLqE detects the low frequency components. As far as aspect (ii) is concerned, the robust weights implied by TC-MLqE help us detect observations in the series where the HRV dynamics has abrupt changes. For example, the records between observations 500 and 1000 receive weights smaller than 0.6.

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

149

Fig. 10. Robust TC-MLqE weights for Patient #chf01. HRV series with and without AR(3) contamination (top) and corresponding normalized weights w ˆ i∗ = w ˆ i / maxi w ˆ i (bottom). The shaded areas correspond to AR(3) contaminations.

Fig. 11. Spectral density estimates for Patient #chf01. Estimated densities based on original sample (solid lines) and sample with additional AR(3) noise (dashed lines) for Gaussian PMLE, HMLE with b = 1.345, and Gaussian TC-MLqE with q = 0.95.

A closer inspection of the HRV trajectory shows that these observations are related to abrupt changes in the HRV pattern. We conjecture that the TC-MLqE can detect low frequency oscillations by setting relatively small weights for such observations. Sensitivity analysis via noise stress test. A crucial aspect in HRV analysis is the stability of spectral estimates in the presence of noise, since lack of stability is the principal source of error in arrhythmia detection and heart diseases classification. To test

150

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

the stability of our findings, we perform a standard sensitivity analysis, following the protocol applied to test the arrhythmia detectors illustrated in Clifford et al. (2012). We follow 3 steps: (i) we contaminate the original HRV series by adding realistic noise, similar to that typically encountered in ECG recordings; (ii) we analyze original and noisy data by a given method; (iii) we compare the stability of the outcomes from noisy and clean samples. In our test, we use noise generated by an AR(3) with parameter θ = (0.8, −0.6, 0.1), and σ 2 = 0.01, which mimics oscillations due to electrode motion artifacts (Goldberger et al., 1988). In Fig. 10, we show the HRV dynamics and TC-MLqE weights based on TC-MLqE with and without noise. The shaded areas emphasize observations with additional noise. The TC-MLqE detects most of the additional artificial records, by assigning very small weights to such observations. The reduced weights imply more stable fits for both HRV dynamics and HRV spectral density. This can be clearly seen in Fig. 11, where we compare the spectral densities corresponding to estimates based on original (solid lines) and noisy data (dashed lines). The spectral density implied by TC-MLqE is fairly stable and it preserves the main low frequency behavior and the RSA. In comparison the PMLE is unstable in the presence of noise. 6. Final remarks This paper advocates the use of robust statistical procedures for the analysis of HRV times series and proposes the TCMLqE for time series models. In HRV applications, robust model fitting should be preferred to other pre-processing methods routinely used in combination with PMLE. The latter approach typically produces biased estimates, as well as loss of medical evidence related to representative outliers. We show that the TC-MLqE represents a reliable pattern recognition tool in the presence of ECG signal anomalies. The results presented in Sections 4 and 5 on real and simulated data suggest that our approach exploits correctly the information in HRV data by detecting influential outliers and reducing their negative impact on parameter estimates. Differently from other classic robust estimators for time series, the TC-MLqE does not require intensive multivariate Monte Carlo integration needed to recenter the estimating Eqs. (2.2). This simplicity is preserved also in the multivariate setting for standard Gaussian models, thus making the TC-MLqE a suitable method to be implemented in a fast software package for real-time HRV analysis. This would be a valuable tool in cardiology, but other biomedical disciplines would also benefit from such a software. For instance, robust inference for model (2.1) can be applied to electroencephalogram, or electrogastrogram signals (Verhagen et al., 1999; Nagarajan et al., 2007; Kim et al., 2012; Park et al., 2013). Our framework is compatible with a number of extensions potentially useful for HRV analysis. For example, Model (2.1) can be improved by considering the heteroscedasticity in Yt ; this would help our understanding of the heart behavior during failure events, where it is reasonable to assume that both the conditional HRV mean and variance depend on the most recent observations. Finally, here the focus is on univariate ECG signals, but our method suits well multivariate HRV time processes. For instance, our approach can be applied to vector autoregressive moving average (VARMA) models for the analysis of multivariate ECG signals for discriminant analysis, as in Maharaj and Alonso (forthcoming). Acknowledgments We thank Martin Gael, Fashid Vahid-Araghi, Mervyn Silvapulle, David Harris, Don Poskitt, Peter Phillips, and all the other participants to the 2013 Time Series workshop at Monash University, for useful comments. Davide La Vecchia is particularly grateful to Damian Hutter, medical doctor at Pediatric Critical Care Medicine and Pediatric Cardiology, Berne University Children’s Hospital, for helpful discussions about the application of statistical procedures in cardiology. References Basu, A., Harris, I.R., Hjort, N.L., Jones, M.C., 1998. Robust and efficient estimation by minimising a density power divergence. Biometrika 85, 549–559. Behar, J., Oster, J., Li, Q., Clifford, G.D., 2013. ECG signal quality during arrhythmia and its application to false alarm reduction. IEEE Trans. Biomed. Eng. 60 (6), 1660–1666. Bianco, A.M., Martinez, E.J., Garcia Ben, M., Yohai, V.J., 1996. Robust procedures for regression models with arima errors. In: COMPSTAT 96, Proceedings in Computational Statistics, pp. 27–38. Bollerslev, T., Wooldridge, J.M., 1992. Quasi-maximum likelihood estimation and inference in dynamic models with time-varying covariances. Econometric Rev. 11 (2), 143–172. Chambers, R.L., 1984. Outlier robust finite population estimation. J. Amer. Statist. Assoc. 81, 1063–1069. Christini, D.J., Bennett, F.H., Lutchen, K.R., Ahmed, H.M., Hausdorff, J.M., Oriol, N., 1995. Application of linear and nonlinear time series modeling to heart rate dynamics analysis. IEEE Trans. Biomed. Eng. 42 (4), 411–415. Clifford, G.D., Behar, J., Li, Q., Rezek, I., 2012. Signal quality indices and data fusion for determining clinical acceptability of electrocardiograms. Physiol. Meas. 33, 14–19. Drost, F.C., Klaassen, C., Werker, B., 1997. Adaptive estimation in time-series models. Ann. Statist. 25, 786–817. Ferrari, D., La Vecchia, D., 2012. On robust estimation via pseudo-additive information. Biometrika 99 (1), 238–244. Ferrari, D., Yang, Y., 2010. Maximum Lq -likelihood estimation. Ann. Statist. 38 (2), 753–783. Fox, A.J., 1972. Outliers in time series. J. R. Stat. Soc. Ser. B Stat. Methodol. 350–363. Goldberger, A.L., Amaral, L.A.N., Glass, L., Hausdorff, J.M., Ivanov, P.Ch., Mark, R.G., Mietus, J.E., Moody, G.B., Peng, C.-K., Stanley, H.E., 2000. PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. Circulation 101 (23), e215–e220. Goldberger, A.L., Rigney, D.R., Mietus, J., Antman, E.M., Greenwald, S., 1988. Nonlinear dynamics in sudden cardiac death syndrome: heart rate oscillations and bifurcations. Cell. Mol. Life Sci. 44 (11), 983–987. Gourieroux, C., Monfort, A., Trognon, A., 1984. Pseudo maximum likelihood methods: theory. Econometrica 52, 681–700.

D. La Vecchia et al. / Computational Statistics and Data Analysis 82 (2015) 137–151

151

Guo, M., Huang, M.N.L., Bai, Z., Hsieh, K.S., 2001. Important ECG diagnosis-aiding indices of ventricular septal defect children with or without congestive heart failure. Stat. Med. 20 (7), 1125–1141. Haggan, V., Ozaki, T., 1981. Modelling nonlinear random vibrations using an amplitude-dependent autoregressive time series model. Biometrika 68 (1), 189–196. Hampel, F.R., Ronchetti, E.M., Russeeuw, P.J., Sthael, W.A., 1986. Robust Statistics: The Approach Based on Influence Functions. Wiley, New York. Kannathal, N., Acharya, U.R., Min, L.C., Suri, J.S., 2007. Prediction of cardiac signals using linear and nonlinear techniques. In: Advances in Cardiac Signal Processing. Springer, pp. 83–107. Kaufmann, T., Süterlin, S., Schultz, S., Vogle, C., 2011. Artiifact: a tool for heart rate artifact processing and heart rate variability analysis. Behav. Res. Methods http://dx.doi.org/10.3758/s13428-011-0107-7. Kim, J.H.K., Pullan, A.J., Bradshaw, L.A., Cheng, L.K., 2012. Influence of body parameters on gastric bioelectric and biomagnetic fields in a realistic volume conductor. Physiol. Meas. 33 (4), 545. Künsch, H., 1984. Infinitesimal robustness for autoregressive processes. Ann. Statist. 12, 843–863. La Vecchia, D., Trojani, F., 2010. Infinitesimal robustness for diffusions. J. Amer. Statist. Assoc. 105, 703–712. Li, T.H., 2008. Laplace periodogram for time series analysis. J. Amer. Statist. Assoc. 103 (482), 757–768. Lin, T., Tian, S., 2012. An ECG signal processing system based on matlab and MIT-BIH. In: Recent Advances in Computer Science and Information Engineering. Springer, pp. 787–792. Lu, K., 2012. An efficient and robust analysis of covariance model for baseline adjustment in parallel-group thorough QT/QTc studies. Stat. Med.. Machado, J.A.F., 1993. Robust model selection and M-estimation. Econometric Theory 9 (03), 478–493. Maharaj, E.A., Alonso, A.M., 2014. Discriminant analysis of multivariate time series using wavelets. Comput. Statist. Data Anal. (forthcoming). Malik, M., Bigger, J.T., Camm, A.J., Kleiger, A., Malliani, R.E., Moss, A., Schwartz, P., 1996. Heart rate variability standards of measurement, physiological interpretation, and clinical use. Eur. Heart J. 17 (3), 354–381. Mancini, L., Ronchetti, E., Trojani, F., 2005. Optimal conditionally unbiased bounded-influence inference in dynamic location and scale models. J. Amer. Statist. Assoc. 100 (470), 628–641. Maronna, R.A., Martin, D.R., Yohai, V.J., 2006. Robust Statistics: Theory and Methods. Wiley, New York. Morillo, C.A., 2012. Heart Rate Variability (HRV) Signal Analysis: Clinical Applications. CRC Press. Muler, N., Peña, D., Yohai, Víctor J., 2009. Robust estimation for ARMA models. Ann. Statist. 37 (2), 816–840. Nagarajan, S.S., Attias, Hagai T., Hild, I.I., Kenneth, E., Sekihara, K., 2007. A probabilistic algorithm for robust interference suppression in bioelectromagnetic sensor data. Stat. Med. 26 (21), 3886–3910. Park, S.A., Hwang, H.J., Lim, J.H., Choi, J.H., Jung, H.K., Im, C.H., 2013. Evaluation of feature extraction methods for EEG-based brain–computer interfaces in terms of robustness to slight changes in electrode locations. Med. Biol. Eng. Comput.. Perakakis, P., Joffily, M., Taylor, M., Guerra, P., Vila, J., 2010. Kardia: a matlab software for the analysis of cardiac interbeat intervals. Comput. Methods Programs Biomed. 98 (1), 83–89. Rodríguez-Liñares, L., Méndez, A.J., Lado, M.J., Olivieri, D.N., Vila, X., Gómez-Conde, I., 2011. An open source tool for heart rate variability spectral analysis. Comput. Methods Programs Biomed. 103 (1), 39–50. Ronchetti, E., 1997. Robustness aspects of model choice. Statist. Sinica 7, 327–338. Ronchetti, E., Trojani, F., 2001. Robust inference with GMM estimators. J. Econometrics 101 (1), 37–69. Salibian-Barrera, M., Van Aelst, S., 2008. Robust model selection using fast and robust bootstrap. Comput. Statist. Data Anal. 52 (12), 5121–5135. Salibian-Barrera, M., Yohai, V., 2006. A fast algorithm for S-regression estimates. J. Comput. Graph. Statist. 15 (2). Shi, P., Tsai, C.L., 1998. A note on the unification of the akaike information criterion. J. R. Stat. Soc. Ser. B Stat. Methodol. 60 (3), 551–558. Spangl, B., Dutter, R., 2005. On robust estimation of power spectra. Austral. J. Statist. 43, 199–210. Spangl, B., Dutter, R., 2012. Analyzing short-term measurements of heart rate variability in the frequency domain using robustly estimated spectral density functions. Comput. Statist. Data Anal. 56 (5), 1188–1199. Tsallis, C., 1988. Possible generalization of Boltzmann–Gibbs statistics. J. Stat. Phys. 52 (1–2), 479–487. Verhagen, M., Van Schelven, L.J., Samsom, M., Smout, André J.P.M., 1999. Pitfalls in the analysis of electrogastrographic recordings. Gastroenterology 117 (2), 453–460. Windham, M.P., 1995. Robustifying model fitting. J. R. Stat. Soc. Ser. B Stat. Methodol. 57, 599–609.