Disentangling synchrony from serial dependency in paired event time series
DDisentangling synchrony from serial dependency in paired event time series
Adrian Odenweller
1, 2, 3, ∗ and Reik V. Donner
1, 4, † Potsdam Institute for Climate Impact Research (PIK), Germany Center for Earth System Research and Sustainability (CEN), University of Hamburg, Germany The Land in the Earth System, Max Planck Institute for Meteorology, Hamburg, Germany Department of Water, Environment, Construction and Safety,Magdeburg–Stendal University of Applied Sciences, Magdeburg, Germany (Dated: March 20, 2020)Quantifying synchronization phenomena based on the timing of events has recently attracted agreat deal of interest in various disciplines such as neuroscience or climatology. A multitude ofsimilarity measures has been proposed for this purpose, including Event Synchronization (ES) andEvent Coincidence Analysis (ECA) as two widely applicable examples. While ES defines synchronyin a data adaptive local way that does not distinguish between different time scales, ECA requiresselecting a specific scale for analysis. In this paper, we use slightly modified versions of both ES andECA that address previous issues with respect to proper normalization and boundary treatment,which are particularly relevant for short time series with low temporal resolution. By numeri-cally studying threshold crossing events in coupled autoregressive processes, we identify a practicallimitation of ES when attempting to study synchrony between serially dependent event sequencesexhibiting event clustering in time. Practical implications of this observation are demonstrated forthe case of functional network representations of climate extremes based on both ES and ECA,while no marked differences between both measures are observed for the case of epileptic electroen-cephalogram (EEG) data. Our findings suggest that careful event detection along with diligentpreprocessing is recommended when applying ES while less crucial for ECA. Despite the lack of ageneral modus operandi for both event definition and detection of synchronization, we suggest ECAas a widely robust method, especially for time resolved synchronization analyses of event time seriesfrom various disciplines.
PACS numbers: 05.45.Tp, 05.45.Xt, 92.60.Ry, 87.19.La
I. INTRODUCTION
Along with the rising availability of empirical dataacross many scientific fields, in the past decades a varietyof statistical methods have been newly developed to dealwith ever larger datasets. As specific events, particularlyextremes, in both nature and society attract a great dealof attention from the academic world as well as the gen-eral public [1], a methodologically sound analysis of eventtime series in general, and synchrony between event inparticular, is crucial not only for research progress, butalso for informed decision making relying on confident re-sults. Accordingly, in different fields of science, the quan-tification of event synchronization has recently become afocal point of a plethora of different studies and meth-ods. Among the methodological developments aimed atserving this purpose, event synchronization (ES) [2] andevent coincidence analysis (ECA) [3, 4] stand out as twoconceptually simple nonlinear measures that are poten-tially applicable to a broad variety of problems of suchdiverse fields like neuroscience and climatology.Being originally motivated by the emerging nonlin-ear dynamical analysis of electroencephalogram (EEG)recordings in terms of spike train synchrony between dif-ferent brain areas [2, 5–12], ES has recently been applied ∗ [email protected] † [email protected], [email protected] to problems outside the neurosciences as well, includinggroup dynamics in both humans and animals [13, 14],econophysics [15], and climate extremes [16–28]. No-tably, the thorough application to climate problems hasbeen mainly governed by a methodological combinationwith the paradigm of functional network analysis [29, 30],as will be further detailed in the course of this work. EShas the important advantage of automatically classifyingpairs of events at two distinct spatial locations as (not)synchronized without the need to manually select any al-gorithmic parameters, particularly a maximum tolerablemutual delay to consider two events synchronized.On the other hand, ECA has been recently introducedbased on the general idea of capturing event synchronybetween point processes, which do not necessarily sharethe common properties of neuronal spike trains like a rel-atively well expressed pacemaker. Successful applicationsof the method can be found across various disciplines, in-cluding paleoclimatology [3], the climate-security nexus[4, 31], plant sciences [32–36], modern day climatology[37, 38] and even seismology [39]. As opposed to ES,ECA commonly requires at least one input parameter(the maximum possible delay) to be selected, which how-ever also entails the potential advantage of a more refinedanalysis by isolating certain time scales based on a pri-ori knowledge or specific research questions. Unlike ES,ECA has so far not been used to analyze EEG data orgenerate functional network representations of large scalespatiotemporal climate data, making this a novelty worth a r X i v : . [ phy s i c s . d a t a - a n ] M a r exploring in the forthcoming sections.Both ES and ECA share the fundamental property ofbasically counting synchronous events based on pairwisecomparison and subsequent aggregation. Yet, they differin the specification of the tolerance window for identi-fying synchronous events, with ES relying on a dynamic(data adaptive) and hence local approach, while ECA re-quests a static (global) parameter to be selected. Eventhough both methods have been demonstrated to be ap-plicable to a wide range of research problems, they haveexhibited a tendency to remain used by rather disjointscientific communities, as no in-depth comparison, whichelaborates on the (dis)advantages of each method, hasbeen conducted so far.Accordingly, in the present paper, we attempt to pro-vide a thorough comparison between ES and ECA bymeans of numerical simulation of simple coupled stochas-tic processes along with real world applications to two ex-emplary climate and EEG datasets. Firstly, in Section II,we introduce formally correct variants of both associationmeasures that address previous ambiguities in the count-ing procedures, thereby re-establishing proper normal-ization, which is especially relevant for short time serieswith low temporal resolution. In Section III, we considercoupled autoregressive AR(1) processes to demonstratethat ES has structural difficulties to capture synchronyin the case of events that are temporally clustered, i.e.,serially dependent.This clustering is no typical feature of EEG recordingsbut emerges rather commonly in climate datasets [40],as we illustrate with a functional network analysis of theIndian monsoon system that replicates previous studiesand highlights their methodological deficiencies (SectionIV B). Therefore, we argue that previous research resultsneed to be interpreted with caution. On the contrary,we provide evidence that ECA is not markedly affectedby event clustering, but provides the additional benefitof allowing for testing physical hypotheses via systemat-ically varying the associated parameter settings (whichmay be guided by a priori knowledge of the system).Building on these conceptual concerns, we propose ECAas a promising robust alternative to ES if temporal eventclustering cannot be ruled out. By analyzing epilepticrat EEG signals (Section IV C), we find that for time se-ries with relatively regularly spaced peaks, ES and ECAare practically interchangeable. However, if the integra-tion of different time scales into one measure is contex-tually justified and events can be clearly marked-off, ESmay still be the favorable method, since it does not re-quire any parameter selection by the user. In Section V,we discuss the repercussions of event definitions on thechoice of association measures that can be derived fromour numerical results. Finally, some general conclusionsare drawn in Section VI. t il − t il t il +1 t jm − t jm t jm +1 t il − t il − t il +1 − t il t jm − t jm − t jm +1 − t jm τ ijlm Only used for τ ijlm Only used for τ ijlm Time series j Time series i FIG. 1. Schematic illustration of event synchronization.
II. METHODS
Both ES and ECA provide measures that go beyondsecond-order moments captured by classical correlations.They are based on pairs of event sequences or binaryevent time series as inputs. When being initially giventwo “normal” time series of nonbinary (continuous or dis-crete valued) variables, such sequences are often obtainedby applying a threshold to the underlying time series ata given percentile (other options will be discussed later).We extract the time indices of these threshold exceedanceevents and let t il denote the time of event l in time se-ries i and t jm the time of event m in time series j with l = 1 , , . . . , s i and m = 1 , , . . . , s j , where s i and s j denote the number of events in the respective series. A. Event Synchronization (ES)
As mentioned previously, ES was first introduced byQuian Quiroga et al. [2] as a parameter free method forthe analysis of synchronization phenomena in spiky elec-troencephalography (EEG) data, but has recently beenapplied to other fields of research as well. Fig. 1 schemat-ically illustrates the basic idea behind ES.Two events at t il and t jm are considered synchronizedif they both occur within a certain data adaptive timeinterval of width τ ijlm defined as τ ijlm = 12 min (cid:110) t il +1 − t il , t il − t il − , t jm +1 − t jm , t jm − t jm − (cid:111) , (1)with l = 2 , , . . . , s i − m = 2 , , . . . , s j − τ ijlm is not evaluated for the first and last event in orderto ensure proper consideration of the boundaries.Equation (1) implies that the more rarely events occurin either (or both) of the time series, the larger τ ijlm willbe, so that we refer to it as a dynamic (local) coincidenceinterval . Thus, if events are rare in the vicinity of eitherof the two events, larger deviations from an instantaneouscoincidence might still be considered synchronized. Thefactor of 1 / τ ijlm less or equal than half the minimum inter-eventwaiting time. The dynamic nature of τ ijlm simplifies theseparation of independent events, which in turn resultsin a variety of temporal scales to be captured by a singlemeasure. The trade-off is that, by design, the value of τ ijlm constantly changes between different pairs of events.Counting the number of synchronized event occur-rences in i given an event in j yields c ( i | j ) = s i − (cid:88) l =2 s j − (cid:88) m =2 J ijlm , (2)where J ij is a counting function that incorporates τ ijlm anddepends on whether or not the synchronization condition σ ijlm = (cid:26) < t il − t jm ≤ τ ijlm J ijlm = σ ijlm = 1and σ jim,l − = 0and σ jim +1 ,l = 01 / t il = t jm or σ ijlm = 1 and( σ jim,l − = 1 or σ jim +1 ,l = 1)0 else (4)We note that the counting function in Eq. (4) devi-ates from the original definition of ES and admittedlylooks rather cumbersome. For correct specification, thechanges are inevitable, though, as otherwise erroneousdouble counting might occur. Due to the condition of aninter-event distance that is smaller than, or equal to, thedynamical coincidence interval τ ijlm in both directions, inthe original definition events could be counted twice. Inorder to avoid this, we thus need to check for all eventpairs whether one of the events has already been countedas synchronized in the opposite direction. If this is thecase, a weight of 1/2 is assigned to this pair so that nor-malization is again carried out correctly. Such a situ-ation can only occur if t il − t jm = τ ijlm and the respec-tive events then contribute equally to c ( i | j ) and c ( j | i ).While this correction should always be incorporated, itis especially important for time series with comparablylow temporal resolution (like daily values of some climatevariable). Alternatively, it is possible to exchange ≤ for < in Eq. (3), which has been done in later applications ofthe ES concept [41, 42]. However, this leads to an entirelynew measure, called ‘SPIKE-Synchronization’, with dif-ferent aggregation and normalization. Because our focushere lies on revealing potential shortcomings that resultfrom application of ES in its original form, we leave thesemore recent developments aside and restrict ourselves toa correction of the original ES measure, which has seenextensive use as mentioned earlier. Furthermore, we pre-sume that for high resolution time series the differences will likely be small as the fundamental functioning of thedynamical coincidence interval is left unchanged.By full analogy, we further define c ( j | i ) and infer the strength of event synchronization between i and j as Q ESij = c ( i | j ) + c ( j | i ) (cid:112) ( s i − s j − , (5)which is normalized so that 0 ≤ Q ESij ≤
1, where Q ESij =1 implies complete event synchronization and Q ESij = 0the absence of any synchronized events.For the purpose of generating a functional networkrepresentation of a set of time series, we consider thepairwise ES strength as a generic statistical associationmeasure, the estimated values of which provide the coef-ficients of a matrix Q ES = ( Q ESij ). Since Q ESij as definedabove is symmetric with respect to interchanges between i and j , this matrix is square symmetric and can there-fore be used to construct an undirected network frommultivariate event data (see Section IV A).As a simple example of the proposed modificationsto the original ES definition, consider an alternatingevent time series, e.g., t i ∈ { , , , , , } and t j ∈{ , , , , , } with τ ijlm = 1 ∀ l, m . Using the originaldefinition would yield c ( i | j ) = 3, c ( j | i ) = 4 and there-fore Q ESij = √ · = 1 .
75, which is not normalized as re-quired. Our corrected version instead yields c ( i | j ) = 1 . c ( j | i ) = 2 and thus Q ESij = 0 . Q ESij will always be lessthan 1 because the first event in one time series (herein i ) cannot be synchronized to any event in the othertime series (here j ) in the calculation of c ( i | j ). The ESstrength between two completely synchronized event se-quences, however, always yields 1, e.g., Q ESii = Q ESjj = 1.
B. Event Coincidence Analysis (ECA)
Event Coincidence Analysis (ECA) [3, 4] is a recentlydeveloped method to measure similarities between eventtime series by incorporating static (global) coincidenceintervals (as opposed to the dynamic coincidence intervalof ES) and potentially also mutual time lags (which havebeen rarely considered along with ES in the past, butcould be included here as well). Fig. 2 illustrates thegeneral idea of ECA.In the following, we will explain the estimation of eventcoincidence rates in full detail only for the case whereevents in j precede events in i . Event coincidence ratesfor the other direction are obtained by simply exchang-ing i and j in all formulas. Again, let t il denote the timeindex of event l in series i and t jm the time index of event m in series j , but now again with l = 1 , , . . . , s i and m = 1 , , . . . , s j . All events have been observed duringthe same observation period [ t , t f ]. We define an in-stantaneous event coincidence if two events at t il and t jm t t il − t il t il +1 t f t t jm − t jm t jm +1 t f τ ∆ Tτ + ∆ T τ + ∆ T Truncated for precursor co-incidence rate r p (∆ T, τ, i | j )Truncated for trigger coinci-dence rate r t (∆ T, τ, i | j )Time series j Time series i FIG. 2. Schematic illustration of event coincidence analysis. occur within a certain static coincidence interval ∆ T sothat 0 ≤ t il − t jm ≤ ∆ T . Accordingly, a lagged event coin-cidence is a coincidence between a time shifted event at t il − τ and an event at t jm ≤ t il − τ , with a given time lag τ ≥
0, implying that 0 ≤ ( t il − τ ) − t jm ≤ ∆ T holds.To quantify the degree of synchrony between event inthe two time series i and j , we distinguish between so-called precursor and trigger event coincidence rates [4,31]. This distinction relates to the question whether thenumber of events in i or j are used for normalization: theprecursor rate quantifies the fraction of events in i thathave been preceded by at least one event in j , while thetrigger rate gives the fraction of events in j than havebeen followed by at least one event in i .Formally, the precursor event coincidence rate is thusdefined as r p ( i | j ; ∆ T, τ ) =1 s i − s (cid:48) i s i (cid:88) l =1+ s (cid:48) i Θ (cid:40) s j (cid:88) m =1 [0 , ∆ T ] (cid:0) ( t il − τ ) − t jm (cid:1)(cid:41) . (6)Since it measures the fraction of ( τ -lagged) events in i that are preceded by at least one event in j , multipleevents in j falling within the static coincidence interval∆ T relative to an event in i are only counted once. Here,Θ( • ) denotes the left-continuous Heaviside step functionwith Θ(0) = 0, which rules out double counting. 1 I ( • )represents the indicator function of the interval I definedas 1 I ( x ) = (cid:26) x ∈ I x / ∈ I . (7)Here, the integer s (cid:48) i denotes the number of events in i that occur within the interval [ t , t + τ + ∆ T ], i.e., s (cid:48) i = s i (cid:88) l =1 [ t ,t + τ +∆ T ] ( t il ) . (8)By full analogy, we define r p ( j | i ; ∆ T, τ ) as the precursorevent coincidence rate for the case of events in j beingpreceded by at least one event in i . In turn, by taking the events in j as the basis for nor-malization, we obtain the trigger event coincidence rate r t ( i | j ; ∆ T, τ ) =1 s j − s (cid:48)(cid:48) j s j − s (cid:48)(cid:48) j (cid:88) m =1 Θ (cid:40) s i (cid:88) l =1 [0 , ∆ T ] (cid:0) ( t il − τ ) − t jm (cid:1)(cid:41) , (9)which measures the fraction of events in j that are fol-lowed by at least one event in i . Here, s (cid:48)(cid:48) j counts thenumber of events in j that occur within the interval[ t f − ( τ + ∆ T ) , t f ]: s (cid:48)(cid:48) j = s j (cid:88) m =1 [ t f − ( τ +∆ T ) ,t f ] ( t jm ) (10)The primes in s (cid:48) i and s (cid:48)(cid:48) j are intended to avoid confusiononce the indices are swapped to calculate the same typesof coincidence rates in opposite directions. One primerefers to the number of events in the interval at the be-ginning of the time series, while two primes refer to thenumber of events in the interval at the end. Finally, wecan define r t ( j | i ; ∆ T, τ ) as the trigger coincidence ratefor the opposite case of events in i preceding events in j .Clearly, being defined as fractions, all four possible eventcoincidence rates are confined to values between 0 and 1.The idea to truncate the succeeding time series at thebeginning for the precursor event coincidence rate andthe preceding time series at the end for the trigger eventcoincidence rate provides a necessary correction of theoriginal ECA definition [4]. Ignoring this may again leadto an erroneous normalization especially if τ (cid:54) = 0, be-cause for the precursor event coincidence rate, events in i before t + τ can never coincide with any event in j .Similarly, for the trigger event coincidence rate, eventsin j after t f − τ can never coincide with any event in i . Ifthis is disregarded, the respective event coincidence ratemay end up at a value below 1, even if all events thatcould possibly coincide also do coincide. This might leadto an underestimation, but at least not to values largerthan 1 as for the uncorrected ES definition. The com-mitted error without this correction vanishes with longtime series as the number of events becomes sufficientlylarge. However, for finite event sequences, proper trun-cation should be employed.Altogether, ECA yields four event coincidence ratesfor every pair of event time series, namely the precur-sor and trigger event coincidence rates in both direc-tions. In what follows, we will only use the trigger eventcoincidence rates as the differences with respect to theprecursor event coincidence rates have been found to becommonly very small across all numerical investigationspresented in the following. However, it should not beargued that this is necessarily always the case by con-struction. Consider, as a simple thought example, twoevent sequences with the same number of events. In thefirst series, all events occur with a precise clock after afixed waiting time, and every second of those events isfollowed by two events in the other series that occur inclose succession (within ∆ T ). In this case, each event ofthe second series has a precursor in the first series, whileonly every second event of the first series triggers eventsin the second one.In order to obtain a single statistical association mea-sure Q ECAij as the degree of event synchrony , we can em-ploy either the mean or the maximum value of the twodirected trigger event coincidence rates r t ( i | j ; ∆ T, τ ) and r t ( j | i ; ∆ T, τ ), denoted as Q ECA,meanij or Q ECA,maxij , re-spectively. The former is especially appropriate for bidi-rectional dependencies, whereas the latter emphasizes thestrengths of unidirectional dependencies more strongly,irrespective of the direction.Finally, similar to the ES, we can create a similaritymatrix from pairwise values of event coincidence ratesif more than two simultaneously measured time seriesare available. The resulting similarity matrix Q ECA =( Q ECAij ) of either the mean or maximum values is againnormalized to 0 ≤ Q ECAij ≤ c ( i | j ) and c ( j | i ) [2, 16, 17]. In turn, for ECA either a similar dif-ference between any of the two event coincidence ratesunder an exchange of i and j or the distinction betweentrigger and precursor event coincidence rates could beused for a similar purpose. We outline correspondingfurther theoretical studies as a subject of future work. III. NUMERICAL STUDY: BIVARIATE AR(1)PROCESSES
We first compare ES and ECA for artificial data stem-ming from a simple bivariate first order vector autoregres-sive (VAR(1)) process. A comprehensive, yet not com-pletely exhaustive, comparison of different synchroniza-tion measures for dynamical systems has been providedin Kreuz et al. [8], but only included ES and not ECA,and only compared the different approaches with regardto coupling and noise strength, but not to serial depen-dency and particularly event clustering, which remainsour focus in the following.Fig. 3 shows a stylized example where events tend tooccur in pairs at subsequent time steps. Evidently, thereis some form of lagged synchronization since a sequenceof one to three consecutive events in j is always followedby an event pair in i . However, ES is unable to detectthis type of synchronization, because the dynamical coin-cidence interval τ ijlm always remains small due to the shortinter-event distances. The illustrative result of Q ESij = 0is unsatisfactory and highlights a major caveat of ES: itfails to properly unravel different temporal scales if eventsare clustered in time. By contrast, given appropriate pa-rameter values for ∆ T and τ , ECA would clearly detect Time series j Time series i FIG. 3. Schematic illustration of an example of clusteredevent time series. a directional relationship so that Q ECAij = 1.Although, for a juxtaposition of two nonlinear meth-ods, it might appear sensible to analyze coupled nonlin-ear systems, a simple linear VAR(1) process already suf-fices to demonstrate the fundamental differences betweenES and ECA with regard to event clustering. This ap-proach is also much easier to comprehend, whereas non-linear systems often defy straightforward imagination.Therefore, we use a VAR(1) model x t = ϕ x t − + κ y t − + ε ,t (11) y t = ϕ y t − + κ x t − + ε ,t (12)to create artificial time series that depend on autoregres-sive parameters ϕ and ϕ (modeling serial correlations)and coupling parameters κ and κ (modeling cross cor-relations) for two variables x t and y t . The error terms ε ,t and ε ,t follow two independent standard normal dis-tributions with mean µ = 0 and variance σ = 1, i.e., ε ,t ∼ N (0 ,
1) and ε ,t ∼ N (0 , x and y are also drawn from two independentstandard normal distributions.For a given realization of this bivariate stochastic pro-cess, the time steps where the associated values exceeda given percentile threshold yield two event time series,which can be used as an input for both ES and ECA.As an illustrative example, for ECA we set the param-eters to τ = 0 and ∆ T = 3. We simulate with a timeseries length of 1,000 points and threshold at the 90thpercentile, yielding s i = s j = 100 events per time se-ries. We consider 1,000 independent runs and calculatethe averages of Q ESij , Q ECA,meanij and Q ECA,maxij over i, j = 1 , . . . , , Q ES , Q ECA,mean and Q ECA,max , respec-tively, all of which are functions of ϕ , ϕ , κ and κ . Fol-lowing our previous considerations (see Fig. 3), we expectthat ES will show an undesirable negative dependencyon ϕ and ϕ . However, simulating the time series fora discrete set of varying parameters with ϕ , ϕ , κ , κ ∈{ , . , . , . , . . . , . , } together with the choice ofusing either Q ECA,meanij or Q ECA,maxij entails (too) manydegrees of freedom for a meaningful analysis. In orderto reduce computational efforts and focus on the mostrelevant aspects, we consider here only two illustrativesettings, one simplified bivariate and one extreme uni-variate case.Fig. 4 shows four 3D surface plots of the ensemblemeans Q ES and Q ECA,mean in dependence of the AR(1)and coupling parameters. The top row with subplots(a) and (b) contains the simplified bidirectional case, forwhich we set ϕ := ϕ = ϕ and κ := κ = κ . The bot-tom row with subplots (c) and (d) contains the extremeunidirectional case, for which we again set ϕ := ϕ = ϕ ,but now κ = 0 so that κ = κ remains the only freecoupling parameter. For brevity, we leave out the resultsof Q ECA,max as they turned out to be qualitatively in-distinguishable from Q ECA,mean . Note that the plots donot share the same z axes as the absolute values are notdirectly comparable among the different measures.Looking at the simplified bidirectional case (Fig. 4a,b),we see that for ϕ + κ ≥ x t and y t diverge and all events occur atsubsequent time steps. This leads to a perfect, albeitmeaningless, synchronization. Much more interesting isthe behavior of ES and ECA in dependence on the AR(1)parameter ϕ , which controls the serial dependency and,hence, temporal clustering of events. For any κ (cid:54) = 0,we observe that Q ES first decreases with ϕ , whereas Q ECA,mean continuously increases with ϕ . This is de-picted more clearly in Fig. 5, where a transect of bothpanels at κ = 0 . ϕ as justified below. On the otherhand, ECA values show a strictly positive dependenceon ϕ . This is a much more understandable and mean-ingful behaviour as increasing ϕ also increases statisticalpersistence, which makes events occur less erratic, butrather in temporal clusters. While this does not justifya stronger synchronization as such, it does lead to thefact that if a synchronization occurs it is more likely toinclude several events at once. As normalization in ECAis carried out over the number of individual events in thetime series, increasing serial dependencies through ϕ ul-timately increases the event coincidence rate (as it is alsocommon in other classical statistical association measureslike Pearson correlation). Such reasoning would also beplausible for ES, but clearly does not become apparentfrom our numerical results.Turning to the extreme unidirectional case (Fig. 4c,d),the strong bias of ES in the presence of serial dependen-cies stands out even more prominently as Q ES decreaseswith respect to ϕ for all values of κ up to ϕ ≈ .
7. Quitecontrarily, the values of Q ECA,mean increase slowly butsteadily for almost all values of κ up to ϕ ≈ .
7. Thesepatterns match the expected behavior, underpinning ourconceptual concerns regarding the potential caveats ofES and providing indications in support of ECA as amore reliable measure of event synchronicity. Yet, for ϕ > ∼ .
7, the results show intriguing patterns. The ESvalues now increase again with rising ϕ , only to dropabruptly at ϕ = 1. The initial increase is likely a re- sult of increasing persistence that eventually leads to theexpected pattern that is commonly observed for ECA,but only commences for very high ϕ in ES. The abruptdrop appears to be a statistical artifact stemming fromthe nonstationarity of the process at ϕ = 1. In turn,for ECA, the obtained values above ϕ > ∼ . ϕ , events in both x t and y t form more andmore clusters in a certain part of the underlying timeseries. For x t , which is entirely independent of y t as κ = 0, this means that the probability of event clustersfalling into the same period as in y t accordingly declines.For y t , even a large unidirectional coupling parameter κ does not curtail this outcome substantially so thatthe overall mean Q ECA,mean also declines. Altogether,the region above ϕ ≈ . ϕ stipulatedabove on the inter-event distances that ultimately leadto the described behavior in both considered cases, wedefine a simple measure for event clustering, the pair-ing coefficient P i , as the fraction of events occurring onsubsequent time steps, P i = 1 s i − s i − (cid:88) l =1 δ (cid:2)(cid:0) t il +1 − t il (cid:1) − (cid:3) , (13)so that 0 ≤ P i ≤
1, where P i = 0 means that no eventpairs form at all and P i = 1 that all events occur on con-secutive time steps. Note that P i is a univariate measurecharacterizing the properties of a single time series i thatessentially entails the first value (∆ t = 1) of the respec-tive empirical inter-event distance distribution. Thus, P i allows us to scrutinize the characteristics of the inputdata used for estimating ES and ECA in a simple man-ner without any notion of coupling or synchronization.Fig. 6 shows the dependency of the ensemble mean of P i , denoted as P , on the AR(1) strength ϕ , with calcula-tions based on 1,000 independent realizations as before.Clearly, P is strictly monotonically increasing with ϕ .At high values of ϕ , events exhibit strong serial depen-dencies, thereby increasing statistical persistence. Thisconfirms our reasoning of using ϕ as a proxy for serial de-pendency in the resulting event time series by providingthe missing, but expected, link. The pairing coefficientwill be drawn upon later again.Even though the absolute difference between ES andECA results might not appear overwhelming for the bi-variate case, it should be noted that the respective quan-tities should only be compared to each other in relativeterms as they are usually ranked internally, before beingreferred to. In functional network analysis, for instance,the values of ES and ECA would be ranked so that onlythe strongest links are included in the network represen-tation [29]. This implies that even small changes in thecorresponding association values might entail large con- A R ( ) s tr e n g t h ϕ . . . . . . C o u p li n g s tr e n g t h κ . . . . . . E n s e m b l e m e a n Q E S . . . . . . . . . (a) A R ( ) s tr e n g t h ϕ . . . . . . C o u p li n g s tr e n g t h κ . . . . . . E n s e m b l e m e a n Q E C A , m e a n . . . . . . . . . (b) A R ( ) s tr e n g t h ϕ . . . . . . C o u p li n g s tr e n g t h κ . . . . . . E n s e m b l e m e a n Q E S . . . . . . . (c) A R ( ) s tr e n g t h ϕ . . . . . . C o u p li n g s tr e n g t h κ . . . . . . E n s e m b l e m e a n Q E C A , m e a n . . . . . . (d) FIG. 4. Ensemble means of the statistical association measures provided by ES and ECA for simplified bivariate and univariatecases of the bivariate VAR(1) model. (a) : Bidirectional coupling, ϕ := ϕ = ϕ , κ := κ = κ , ES. (b) : Bidirectional coupling,ECA. (c) : Unidirectional coupling, ϕ := ϕ = ϕ , κ = 0, ES. (d) : Unidirectional coupling, ECA. sequences for the inferred network structure. Thus, itis not the absolute values that should be compared, butrather the patterns in response to changes in ϕ and κ asjust described.Taken together, our simulation results confirm our ini-tial expectation that the dynamical coincidence interval τ ijlm unambiguously renders ES insensitive to properly de-tect synchronization if the events in a time series arestrongly clustered, which is a common property of cli-mate extremes [43]. This undesirable outcome may ham-per the reliability of results and interpretations obtainedfrom such networks, as we will show in the following. IV. REAL WORLD EXAMPLES
Following our numerical results for the artificial data inSection III, we next attempt to shed some more light onthe real world implications of those findings. Consider-ing the extensive previous research on functional networkanalysis based on ES as a statistical similarity measure,Section IV A concisely reviews the key elements of net-work construction and analysis. Then, in Section IV B,we demonstrate that ES yields biased results for climatenetwork representations of heavy rainfall events, since cli-mate time series commonly exhibit serial dependenciesand clustering among extreme events. Along with repro-ducing some results of previous studies based on ES, we . . . . . . ϕ . . . . . . . . E n s e m b l e m e a n Q E S , Q E C A , m e a n ESECA
FIG. 5. Transect of the simulation results for the symmetri-cally bidirectionally coupled AR(1) processes (Fig. 4a,b) with κ = 0 .
2. Solid lines and shadings indicate the ensemble meansand standard deviations, respectively. . . . . . . ϕ . . . . . . E n s e m b l e m e a n P FIG. 6. Dependency of pairing coefficient P i on AR(1)strength ϕ , shown as ensemble mean (solid line) and stan-dard deviation (shaded band). present substantially different results based on ECA thatare not only more robust in the presence of event cluster-ing, but also allow us to analyze the temporal evolutionof extreme events in a functional precipitation network,enabling worthwhile customized analyses on a more de-tailed level than when using ES.As an illustrative counterexample, Section IV C ana-lyzes five sets of bivariate rat EEG signals, includingone non-epileptic example and four epileptic spike trains,for which ES and ECA yield qualitatively very similarresults. This highlights that due to a relatively nar-row waiting time distribution of clearly discernible EEGspikes (i.e., the existence of a relatively regular inter-nal pacemaker), temporal clustering among events is nota major issue here (i.e., there exist serial dependencies Consider 6 timeseries on a spatialgrid time Extract event timeseries at each location,e.g. by thresholding Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Calculate similaritymatrix Q using e.g.ES or ECA1 2 34 5 6 Threshold similaritymatrix: Link if Q ij > θ Symmetric adjacencymatrix A shows links k k k k k k Calculate and visu-alize network measures,e.g. degree
FIG. 7. Flow chart of functional network analysis for a genericexample of spatially embedded time series. among events, but of an entirely different form than inthe previous numerical example). Furthermore, we com-pare our results in view of a different event definition thatis arguably more directly applicable to the detection oflocal spikes in noisy EEG signals.
A. Functional networks
The combination of nonlinear time series analysis withcomplex network theory is a rapidly growing area of re-search as it allows to disentangle and visualize spatiotem-poral patterns from large scale datasets. Fig. 7 shows aflow chart of how to incorporate ES and ECA into theconstruction of a spatially embedded functional network.While this bears the hallmarks of climate network analy-sis [29, 30], it is straightforward to extend this approachto other applications.A network G = ( V, E ) is defined by a set of vertices (ornodes) V = { , . . . , N } with N = | V | and a set of edges (or links) E ⊆ V × V . The edges E with K = | E | arerepresented in the adjacency matrix A , in which self con-nections conventionally do not exist so that A ii = 0 ∀ i .Having calculated the similarity matrix Q from either ESor ECA, we link i and j if Q ij is above a certain threshold θ . Thus, we obtain a binary square symmetric adjacencymatrix A of an undirected network, where A ij = A ji = 1indicates a link between i and j and A ij = A ji = 0 thelack thereof. As a valuable alternative to choosing a par-ticular value of θ , it is common practice to instead prede-fine a global edge density ρ , which extracts the strongeststatistical associations and thereby indirectly defines θ .Subsequently, the resulting adjacency matrix can beanalyzed by means of complex network theory. Out ofthe abundance of existing network measures [44], for thesake of brevity, in this work we only consider the degree density ρ i = 1 N − N (cid:88) j =1 A ij (14)as the simplest and most prominent vertex measure,yielding the number of connections associated with eachnode, normalized to 0 ≤ ρ i ≤ B. Precipitation data
Following upon previous works on ES based functionalclimate networks for the Indian subcontinent [16, 17, 24],we use gauge adjusted satellite data from the TropicalRainfall Measuring Mission (TRMM)[45] to construct aclimate network representation of extreme precipitationduring the Indian Summer Monsoon (ISM). We select re-sampled daily precipitation sums on a square grid with aspatial resolution of 0 . ◦ × . ◦ ( ∼
28 km) for the period1998–2015 (TMPA 3B42 V7), from which we extract thesummer monsoon season of June to September. For thechosen area from 60 . − . ◦ E and 5 . − . ◦ N,we thus have 22,400 (160 × · et al. [24] we also threshold at the 90th percentile andselect a global edge density of ρ = 0 .
05 as a convenientvalue [29].Although the capability of ES to dynamically incor-porate different time scales at once through τ ijlm appearsworthwhile at first, it entails a major disadvantage byoveremphasizing the node degree of the resulting func-tional climate network in regions where events are tem-porally isolated. This is a consequence of the systematicunderestimation of the degree values in regions whereevents tend to occur temporally clustered and, hence,an immediate manifestation of the adverse effect of tem-poral event clustering on the results of ES as demon-strated in Section III. In order to quantify and visualizethis undesirable property, we calculate the previously in-troduced pairing coefficient P i for all grid points and plotthis alongside the ES degree density ρ i in Fig. 8.Notably, Fig. 8b reproduces the results of Stolbova et al. [24, Fig. 3] almost exactly, with minor differencesoriginating from the described corrections of the ES algo-rithm and a slightly longer time period covered. Further-more, our results are also very similar to those of Malik et al. [16, 17], who used a different data source withoutocean coverage.Although Fig. 8a solely contains local information ofthe dynamical characteristics of events at each individualgrid point without any notion of coupling to other gridpoints whatsoever, we observe interesting similarities incomparison with Fig. 8b. Specifically, in areas wherethe degree density is high, the pairing coefficient very of-ten has low values, and vice versa. In many areas, thetwo figures resemble negative images of each other. This ◦ N ◦ N ◦ N ◦ E ◦ E ◦ E (a) Pairing coefficient P i . . . . . ◦ E ◦ E ◦ E (b) ES degree density ρ i .
000 0 .
025 0 .
050 0 .
075 0 . FIG. 8. Spatial pattern of the pairing coefficient and thedegree density of the ES based functional climate network. . . . . . . . P i . . . . E Sd e g r ee d e n s i t y ρ i P r o b a b ili t y d e n s i t y FIG. 9. Estimate of the joint probability density of the pairingcoefficient and the degree density of the ES based functionalclimate network. holds especially true for regions that had been reported inprevious studies as important for governing monsoon dy-namics at large spatial scales, such as Northern Pakistan,the Tibetan Plateau, or the Eastern Ghats [17, 24]. Notethat even minor artifacts such as interspersed grid pointswith low ρ i values on the Tibetan Plateau can exhibit ahigh pairing coefficient P i . The effect of P i on ρ i is furtherdisplayed in a two-dimensional density plot, see Fig. 9,which reveals a negative relationship between ρ i and pair-ing coefficient up to P i ≈ .
2. Fundamentally, this meansthat a trivial property of the local time series often prede-termines the event synchronization strength to other gridpoints in many areas, especially those deemed pivotal formonsoon dynamics. This suggests that ES may not bean ideal similarity measure for time series of extreme cli-mate events, which frequently entail serial dependenciesand temporal event clustering [43]. Furthermore, our ob-servations are compatible with the artificial data resultsfrom Section III and highlight the practical implicationsof the described weaknesses of ES.On the other hand, the ECA based networks do notexhibit such adverse dependencies of the node degree onthe pairing coefficient. Fig. 10 shows the degree densityfield for four different parameter settings with ∆ T = 5and τ ∈ { , , , } days. Varying τ while fixing ∆ T enables us to clearly isolate lagged synchrony and thus0 ◦ N ◦ N ◦ N ◦ E ◦ E ◦ E (a) ECA ρ i , ∆ T = 5, τ = 00 .
00 0 .
05 0 .
10 0 .
15 0 . ◦ E ◦ E ◦ E (b) ECA ρ i , ∆ T = 5, τ = 20 .
00 0 .
05 0 .
10 0 .
15 0 . ◦ N ◦ N ◦ N ◦ E ◦ E ◦ E (c) ECA ρ i , ∆ T = 5, τ = 50 .
00 0 .
05 0 .
10 0 .
15 0 . ◦ E ◦ E ◦ E (d) ECA ρ i , ∆ T = 5, τ = 100 .
00 0 .
05 0 .
10 0 .
15 0 . FIG. 10. Spatial pattern of the degree density of the ECAbased functional climate network for four different parametersettings. to analyze the spatiotemporal evolution of synchronousextreme precipitation events. Effectively, this moves thefixed-length synchrony detection window back in time(see Fig. 2). This stands in contrast to fixing τ whilevarying ∆ T , which would not have allowed us to clearlydisentangle lagged from instantaneous synchrony (within∆ T ). The obtained spatial patterns vastly differ fromthose in Fig. 8b, and the additional parameters of ECAeven allow to isolate distinct time scales, thereby enablingus to assess the temporal evolution of heavy precipitationpatterns along the sequence of climate networks. In ouropinion, this provides a valuable extension of the ES ap-proach, where the specific underlying time scales remainunknown, rendering the outcome of functional climatenetwork analysis comparatively opaque. Similarly to theES case, Fig. 11 shows the influence of P i on ρ i for thefour considered parameter settings of the ECA based net-work. In contrast to Fig. 9, the relationship between bothcharacteristics is far less pronounced as the points scattermuch more. The slightly positive dependency in Fig. 11agradually evolves into a slightly negative one in Fig. 11d,but is subject to strong dispersion.Based on our findings as reported above, we tenta-tively propose ECA as the preferable similarity mea-sure for event based functional climate network construc-tion and analysis and stress that previous results basedon ES should be interpreted with caution and carefullyre-examined wherever possible. Apart from enabling aproper disentanglement of synchrony from serial depen- dency, we further advocate ECA’s ability to precisely an-alyze the temporal scales encoded in a given network.Yet, even if such a refined analysis is not desired, theadverse effects of event clustering on ES can still over-shadow the potential benefits of incorporating multipletime scales at once to a large extent. C. EEG data
As a second rather common example for the applica-tion of bivariate event based statistics, we analyze fivepairs of rat EEG signals obtained from electrodes at theleft and right frontal cortex of male adult rats, includ-ing one non-epileptic case and four epileptic cases. Thisdataset was also included in Quian Quiroga et al. [2]and is publicly available [46] with more detailed infor-mation on experiment settings, recordings and results tobe found in [47]. Each signal was recorded for 5 s anddigitized at 200 Hz, yielding a time series of 1,000 values.Fig. 12 shows plots of the five EEG pairs. Example (a)displays the normal non-epileptic signal, while examples(b) to (e) all exhibit clear epileptic activity as apparentfrom regular spike discharges.Under the plausible assumption that neither the spikeshape nor the background noise carry valuable informa-tion, we may again extract events from the depicted timeseries, for which we compare two options. Firstly, similarto the climate data case, we simply impose a threshold atthe 97th percentile to obtain a well-defined discrimina-tion of extreme events. Secondly, we follow the approachof Quian Quiroga et al. [2] by declaring events at t as localmaxima fulfilling the following conditions: (1) x t > x t + k for k = − K, . . . , − , , . . . , K (also including k = 0 asstated in the original reference appears incorrect since x t > x t would never be satisfied) and (2) x t > x t ± K + h ,with free parameters K and h . In agreement with theoriginal definition, we set K = 3 and h = 0 .
1. Subse-quently, we calculate Q ESij and Q ECA,meanij for both typesof event definitions.Our results do not exactly replicate those inQuian Quiroga et al. [2] as they replaced the dynamic lo-cal ES coincidence interval τ ijlm with a fixed global valueof 2 time steps, i.e., 10 ms. In our opinion, this prac-tice of educated guessing should be handled with careas it might void the normalization by enabling uninten-tional double counting of event pairs. However, to allowfor some degree of comparison, we set ∆ T = 10 ms and τ = 0 ms for the ECA parameters, which can be donewithout jeopardizing normalization.Table I displays the results of ES and ECA for bothevent definitions and all five EEG pairs. Since ES andECA merely share a normalization to the same intervalbut differ substantially between these boundaries, abso-lute values are not strictly comparable as also mentionedpreviously. In order to allow for a fair comparison, wethus also provide the percentage values in one columnas compared with the maximum value of that column.1 . . . . . . . . E C A D e g r ee d e n s i t y ρ i (a) ∆ T = 5, τ = 0 . . . . . . . . (b) ∆ T = 5, τ = 2 . . . . P i . . . . E C A D e g r ee d e n s i t y ρ i (c) ∆ T = 5, τ = 5 . . . . P i . . . . (d) ∆ T = 5, τ = 10 P r o b a b ili t y d e n s i t y FIG. 11. Estimate of the joint probability density of the pairing coefficient and the degree density of the ECA based functionalclimate networks. − − m V (a) − − − (b) − − (c) − − m V (d) − − (e) FIG. 12. The five considered examples of rat EEGs, with (a) a normal non-epileptic signal and (b)-(e) epileptic spike trains.Left hemisphere signals (green) are plotted with a vertical offset ( − . − − TABLE I. ES and ECA results for five selected rat EEGs.Threshold exceedance Method from [2]Case Q ESij Q ECA,meanij Q ESij Q ECA,meanij (a) 0.22 (31%) 0.32 (43%) 0.70 (78%) 0.55 (79%)(b) 0.50 (70%) 0.63 (85%) 0.77 (85%) 0.53 (76%)(c) 0.21 (30%) 0.23 (31%) 0.77 (85%) 0.51 (74%)(d) 0.43 (60%) 0.45 (60%) 0.85 (94%) 0.57 (81%)(e) 0.71 (100%) 0.75 (100%) 0.91 (100%) 0.70 (100%)
Note that this relative approach is also consistent withour strategy used when studying the artificial simulationand the real world climate data. In the first case rela-tive patterns rather than absolute values were analyzed,while in the second case the network construction processextracted only the strongest links on a relative basis (seeFig. 7).In relative terms, the differences between ES and ECAturn out to be small for both event definitions, with amaximum difference of 15% for example (b). All otherexamples exhibit even smaller differences. Irrespectiveof either the event definition or the similarity measure,within the epileptic signals (e) is always ranked first(most strongly synchronized) and (c) last (least synchro-nized). Within each event definition method and againexcept for example (a), the ranking order remains con-sistent across both ES and ECA. This confirms that theresults of both approaches resemble each other fairly wellfor time series that are characterized by relatively con-stant event rates. However, this observation even holdstrue for the non-epileptic case (a), where results are alsocomparable for both event definitions. Yet, we do observepronounced differences between the two event definitions,which are discussed in more detail below. Of course, theECA results depend markedly on the parameter valuesand similar values can be obtained by setting the delay τ sufficiently close to the mean inter-event distance. How-ever, admittedly, the potential to analyze the temporalevolution of event synchrony might not be regarded as anequally important feature for EEG as for climate data.Thus, if the two different event definitions are con-sidered separately, ES and ECA yield very similar re-sults. This is a direct and plausible consequence of asufficiently narrow inter-event distance distribution forepileptic EEG spike trains in examples (b)-(e), whichstands in marked contrast to the precipitation data usedpreviously. Since ES was originally designed with EEGapplications in mind, our findings are conceptually justi-fied and not surprising. While the observed consistencyof ES and ECA was also facilitated by very distinct eventsdue to the recorded epileptic activity in (b)-(e), notablythis also holds true for the non-epileptic example (a),most likely because the probability distribution of wait-ing times between subsequent events hardly shows anyvalues close to zero (not shown), which in turn wouldcorrespond to a regime where we might expect devia-tions between the two methods. This finding further un- derlines the versatility of ECA.The importance of unambiguous event extraction ismoreover revealed in examples (a) and (c). In exam-ple (a), the threshold method only yields results of 0.22(31%) and 0.32 (43%), while the event definition from [2]leads to values of 0.70 (78%) and 0.55 (79%). Similarly,in example (c) we obtain 0.21 (30%) and 0.23 (31%) us-ing threshold events versus 0.77 (85%) and 0.51 (74%)using [2]. In both examples, these differences in bothabsolute, but more importantly also relative values, arevery likely caused by less pronounced peaks over a dy-namic background in (a) and (c) as compared with theother examples, which probably led to error-prone eventdefinitions for a simple threshold method. This alreadyhints to the overarching issue of statistically disentan-gling events from underlying phenomena, which will befurther discussed below in conjunction with the subse-quent choice of a proper similarity measure. V. DISCUSSION
Within the scope of (extreme) event analysis, the prob-lem of serial dependencies in time can principally be tack-led in two ways: either by choosing a robust analysismethod or by defining events in different ways. For thefirst approach, we have provided here tentative explana-tions why ECA may outperform ES. However, it appearsalso viable to address the issue already in the previousstep of event definition. This necessarily requires moreinvolved preprocessing techniques than a mere thresholdexceedance strategy.For EEG signals, the temporal resolution relative tothe number of spikes is usually much higher than for cli-matological data. Thus, several time values that clearlybelong to the same peak might fulfill the threshold ex-ceedance criterion, even for very high percentiles. Sincethe focus in EEG spike train analysis is on the extrac-tion of singular events that are local (or relative) maxima,which may well have different amplitudes among them-selves, the event definition method by Quian Quiroga et al. [2] is a sensible approach for this delicate task. Forclimatological extremes, we are however not interested inlocal, but indeed rather in global (or absolute) maximawith respect to some threshold (i.e., taking a peaks-over-threshold approach as common in extreme value theory)because these are the type of events with potentially dev-astating impacts. Thus, applying said EEG event defi-nition method also to climate data appears unreason-able. Yet another possibility would be the utilization ofsophisticated clustering algorithms. However, we rejectthis as unnecessarily complex for the sake of the presentwork, since ECA fulfills the same purpose in a much morestraightforward and easily understandable manner.In a broader context, the task of disentangling eventsynchrony from serial dependency therefore transitionsinto the more profound endeavour of disentangling sta-tistical events from underlying phenomena, which we3are eventually interested in. In our opinion, this mustbe either informed by a priori knowledge of the systemor guided by specific research questions targeted to thegiven dataset under study. A shared feature of the twoaforementioned options is that they both require someexternal parameters, which determine the expected timescales of serial dependencies and which cannot be set in-dependent of data and context. In essence, the values K and h for the event definition in [2] or ∆ T and τ forECA constitute different parameterizations of just thisissue. Thus, in our opinion, the introduction of a sophis-ticated clustering algorithm [19], which merges severalpreviously defined global extremes into one, may onlyshift the choice of these time scale parameters elsewhere,without relieving us of the actual task of determiningthem.In this line of argument, there seems to be no uni-versally optimal method. While the inclusion of multi-ple time scales into one single association measure forevent sequences might constitute a noteworthy advan-tage of ES over ECA, this can only be truly exploited ifthe data has been preprocessed diligently. Succumbingto the tempting pitfall of using ES without preprocessingmakes it vulnerable to the identified weaknesses originat-ing from serial dependencies among events. ECA, on theother hand, offers a more refined analysis of time scalesthrough ∆ T (and τ ) with the considerable advantage ofdealing properly with both serially dependent and inde-pendent event time series, albeit without the possibilityto dynamically incorporate changing event rates into asingle measure. Using a sliding windows approach couldhowever alleviate this alleged restriction.For future research on functional climate networks, wetherefore reinforce our suggestion to use ECA instead ofES. Even if serial dependencies were dealt with beforeso that ES worked as intended, the algorithm would stilllack a clear declaration on the involved time scales aswell as the possibility to scrutinize the temporal evolu-tion, which we perceive as a valuable advantage in its ownright. Additionally, using ECA elegantly circumvents theneed for declustering along with its ensuing parameteri-zation difficulties. Only if the inclusion of varying timescales into one single measure is essential, if events areclearly delineated and if refined temporal analyses areundesired, ES appears to be the method of choice. Thisis probably the case for EEG data, even though ECA canalso be adjusted to fit such applications as well.At a final note, we emphasize that many other mea-sures have been used in the past for quantifying differentaspects of statistical interdependence between two timeseries, particularly in the context of functional networkanalysis. For example, in the context of spatiotemporalfields of climate data, Pearson correlation and mutualinformation have been often used for this purpose. Tounderstand the differences with respect to ES and ECAused in the present work, one should note that those (aswell as many other) measures take all existing data points(i.e., values from both the bulk and the tails of the prob- ability distributions of the variable of interest) into ac-count and attempt to quantify the strength of a linear(Pearson correlation) or arbitrary functional relationship(mutual information) between two series. Hence, thosemeasures are necessarily dominated by statistical associ-ations among the bulk values due to their by far largernumber than that of the values in the tails. In turn,there are important applications where statistical associ-ations among those bulk values are not of primary inter-est, since the relationship between extraordinary values(or extremes) can be believed to differ from that under“normal conditions” [33] (for example, in the precipita-tion example discussed in Sec. IV B). On purpose, onlythe latter aspect has been addressed in the present work.Due to their conceptual difference (continuum-based ver-sus event-based association measures), inferred statisti-cal associations among the given data sets may cruciallychange when either considering all data or focusing onlyon the extremes (or even more, in our present work, onlythe timing of extremes and not their magnitudes). Inter-comparisons between measures of both types have beenpublished elsewhere (e.g., in [48] for neuroscience appli-cations or [49] for climate data), and we did not attemptto repeat such studies here to maintain the topical focusof the present work. VI. CONCLUSIONS
In this paper, we have explored the key differences oftwo statistical association measures for event time series,event synchronization (ES) and event coincidence analy-sis (ECA). Both measures have been used extensively indifferent disciplines, yet had hardly been systematicallycompared or applied to the same data before [50].First of all, building on identified ambiguities in thetheoretical definitions of both measures, we introducedand implemented corrected versions of the original pro-posals. We then created artificial data from a coupledautoregressive process, by which we were able to provideevidence that ES does not allow to unambiguously disen-tangle event synchrony from serial dependencies, whereasECA is significantly less susceptible to such issues. Re-producing the results of previous studies, we demon-strated ensuing implications for real-world data that re-inforce our conceptual concerns. We specifically providedindications that results from several past climate networkstudies, which rely on ES as a similarity metric, need tobe interpreted with caution as they might exhibit severebiases originating from unaddressed serial dependenciesbetween events. On the other hand, for epileptic EEGdata, we showed that both ES and ECA yield qualita-tively similar results as the characteristic spike trains en-tail relatively constant event rates without notable tem-poral clustering.Furthermore, we discussed the nexus of event defini-tions and appropriate similarity measures conceptually.We argued that disentangling event synchrony from se-4rial dependency is on a lower level tantamount to disen-tangling statistical events from underlying phenomena.While methods that extract local extremes prove to besensible for data with clear spikes of varying amplitudesuch as EEG signals, they are not applicable for caseswhere the focus lies on global maxima such as extremeclimate events. Because ES only works properly for pre-processed data with a priori knowledge, we propose ECAas an alternative measure that can handle both seriallydependent and independent data. Furthermore, ECAallows to explicitly analyze temporal evolution and ele-gantly bypasses the need for unnecessarily complex clus-tering algorithms that would be required if we wanted toanalyze extremes with ES.While both ES and ECA have strength and weak-nesses, the nonparametric nature of ES makes it all tooeasy to succumb to the temptation of omitting the defini-tion of a time scale, within which multiple events belongto the same phenomenon. However, it appears impossibleto leap over this step for synchronization to be attributedbetween meaningfully defined events. Whether it is bestto stipulate this time scale via parameters in the eventdefinition or rather in the subsequent analysis, remainsopen for debate at this point. Even though the questfor a universally optimal method in our view thus consti-tutes an ill-posed problem, we advocate in favor of ECAas providing a straightforward event based statistical as-sociation measure to analyze event time series that mayor may not include serial dependencies, without caveats due to temporal event clustering.Ultimately, the question which event definition andwhich similarity measure is most appropriate remains amatter of choice. But that choice should be made well-informed, imperatively embedded into the research con-text and data characteristics.
ACKNOWLEDGMENTS
This work has been financially supported by the Ger-man Federal Ministry for Education and Research viathe Young Investigators Group CoSy-CC2: Complex Sys-tems Approaches to Understanding Causes and Conse-quences of Past, Present and Future Climate Change(grant no. 01LN1306A) and the Belmont Forum/JPIClimate project GOTHAM (Globally Observed Telecon-nections and their Representation in Hierarchies of At-mospheric Models, grant no. 01LP16MA). AO has beenpartially funded by the Foundation of German Business(sdw) and received travel support from the GraduateSchool of Integrated Climate System Sciences (SICSS)of the Cluster of Excellence CLICCS at the Univer-sity of Hamburg. AO would like to express his grati-tude to Joachim Krug for supporting and co-supervisingparts of the presented work at the University of Cologne.Functional network analysis has been performed usingthe Python package pyunicorn [51] freely available atGitHub. [1] S. Albeverio, V. Jentsch, and H. Kantz,
Extreme eventsin nature and society (Springer, Berlin, 2006).[2] R. Quian Quiroga, T. Kreuz, and P. Grassberger, Phys-ical Review E , 041904 (2002).[3] J. F. Donges, R. V. Donner, M. H. Trauth, N. Marwan,H.-J. Schellnhuber, and J. Kurths, Proceedings of theNational Academy of Sciences , 20422 (2011).[4] J. Donges, C.-F. Schleussner, J. Siegmund, and R. Don-ner, European Physical Journal Special Topics , 471(2016).[5] T. Kreuz, R. G. Andrzejak, F. Mormann, A. Kraskov,H. St¨ogbauer, C. E. Elger, K. Lehnertz, and P. Grass-berger, Physical Review E , 061915 (2004).[6] E. Pereda, R. Q. Quiroga, and J. Bhattacharya, Progressin Neurobiology , 1 (2005).[7] T. Kreuz, J. S. Haas, A. Morelli, H. D. Abarbanel, andA. Politi, Journal of Neuroscience Methods , 151(2007).[8] T. Kreuz, F. Mormann, R. G. Andrzejak, A. Kraskov,K. Lehnertz, and P. Grassberger, Physica D: NonlinearPhenomena , 29 (2007).[9] C. Stam, Clinical Neurophysiology , 2266 (2005).[10] J. Dauwels, F. Vialatte, T. Weber, and A. Cichocki, in International Conference on Neural Information Process-ing (Springer, 2008) pp. 177–185.[11] G. N. Angotzi, F. Boi, S. Zordan, A. Bonfanti, andA. Vato, Scientific Reports , 5963 (2014). [12] R. D. Singh, S. J. Gibbons, S. A. Saravanaperumal,P. Du, G. W. Hennig, S. T. Eisenman, A. Mazzone,Y. Hayashi, C. Cao, G. J. Stoltz, et al. , Journal of Phys-iology , 4051 (2014).[13] G. Varni, G. Volpe, and A. Camurri, IEEE Transactionson Multimedia , 576 (2010).[14] S. Butail, V. Mwaffo, and M. Porfiri, Physical Review E , 042411 (2016).[15] W.-X. Zhou and D. Sornette, Physica A: Statistical Me-chanics and its Applications , 543 (2003).[16] N. Malik, N. Marwan, and J. Kurths, Nonlinear Pro-cesses in Geophysics , 371 (2010).[17] N. Malik, B. Bookhagen, N. Marwan, and J. Kurths,Climate Dynamics , 971 (2012).[18] N. Boers, B. Bookhagen, N. Marwan, J. Kurths, andJ. Marengo, Geophysical Research Letters , 4386(2013).[19] N. Boers, B. Bookhagen, H. M. J. Barbosa, N. Marwan,J. Kurths, and J. A. Marengo, Nature Communications , 5199 (2014).[20] N. Boers, A. Rheinwalt, B. Bookhagen, H. M. Barbosa,N. Marwan, J. Marengo, and J. Kurths, GeophysicalResearch Letters , 7397 (2014).[21] N. Boers, R. V. Donner, B. Bookhagen, and J. Kurths,Climate Dynamics , 619 (2015).[22] N. Boers, B. Bookhagen, J. Marengo, N. Marwan, J.-S.von Storch, and J. Kurths, Journal of Climate , 1031 (2015).[23] N. Boers, B. Bookhagen, N. Marwan, and J. Kurths,Climate Dynamics , 601 (2016).[24] V. Stolbova, P. Martin, B. Bookhagen, N. Marwan, andJ. Kurths, Nonlinear Processes in Geophysics , 901(2014).[25] S.-H. He, T.-C. Feng, Y.-C. Gong, Y.-H. Huang, C.-G.Wu, and Z.-Q. Gong, Chinese Physics B , 059202(2014).[26] K. Rehfeld and J. Kurths, Climate of the Past , 107(2014).[27] A. Rheinwalt, N. Marwan, J. Kurths, P. Werner, andF.-W. Gerstengarbe, in High Performance Computing,Networking, Storage and Analysis (SCC), 2012 SC Com-panion: (IEEE, 2012) pp. 500–505.[28] N. Marwan and J. Kurths, Chaos: An InterdisciplinaryJournal of Nonlinear Science , 097609 (2015).[29] R. V. Donner, M. Wiedermann, and J. F. Donges, in Nonlinear and Stochastic Climate Dynamics , edited byC. Franzke and T. O’Kane (Cambridge University Press,Cambridge, 2017) 1st ed., pp. 159–183.[30] H. A. Dijkstra, E. Hern´andez-Garc´ıa, C. Masoller, andM. Barreiro,
Networks in climate (Cambridge UniversityPress, Cambridge, 2019).[31] C.-F. Schleussner, J. F. Donges, R. V. Donner, and H. J.Schellnhuber, Proceedings of the National Academy ofSciences , 9216 (2016).[32] J. F. Siegmund, T. G. M. Sanders, I. Heinrich, E. van derMaaten, S. Simard, G. Helle, and R. V. Donner, Fron-tiers in Plant Science , 733 (2016).[33] J. F. Siegmund, M. Wiedermann, J. F. Donges, andR. V. Donner, Biogeosciences , 5541 (2016).[34] J. F. Siegmund, N. Siegmund, and R. V. Donner, Com-puters & Geosciences , 64 (2017).[35] A. Rammig, M. Wiedermann, J. F. Donges, F. Babst,W. von Bloh, D. Frank, K. Thonicke, and M. D. Ma-hecha, Biogeosciences , 373 (2015).[36] L. Baumbach, J. F. Siegmund, M. Mittermeier, andR. V. Donner, Biogeosciences , 4891 (2017).[37] M. Wiedermann, J. F. Siegmund, J. F. Donges,J. Kurths, and R. V. Donner, “Differential imprints of distinct enso flavors in global extreme precipitation pat-terns,” (2017), arXiv:1702.00218 [physics.ao-ph].[38] S. Sippel, J. Zscheischler, M. D. Mahecha, R. Orth,M. Reichstein, M. Vogel, and S. I. Seneviratne, EarthSystem Dynamics , 387 (2017).[39] N. Sarlis, Entropy , 561 (2018).[40] F. Wolf, J. Bauer, N. Boers, and R. V. Donner, Chaos , 033102 (2020).[41] T. Kreuz, M. Mulansky, and N. Bozanic, Journal of Neu-rophysiology , 3432 (2015).[42] M. Mulansky, N. Bozanic, A. Sburlea, and T. Kreuz,in (IEEE, 2015) pp. 1–8.[43] J. Kropp and H.-J. Schellnhuber, eds., In Extremis -Disruptive Events and Trends in Climate and Hydrology (Springer, Berlin, 2011).[44] L. d. F. Costa, F. A. Rodrigues, G. Travieso, and P. R. V.Boas, Advances in Physics , 167 (2007).[45] Goddard Earth Sciences Data and Information ServicesCenter, TRMM (TMPA) Precipitation L3 1 day 0.25 de-gree x 0.25 degree V7 , Goddard Earth Sciences Data andInformation Services Center (GES DISC) (2016).[46] The rat EEG data may be downloaded from as dataset 4.[47] G. van Luijtelaar,
The WAG/Rij Rat Model of AbsenceEpilepsy: Ten Years of Research: a Compilation of Pa-pers (Nijmegen University Press, 1997).[48] R. Quian Quiroga, A. Kraskov, T. Kreuz, and P. Grass-berger, Phys. Rev. E , 041903 (2002).[49] L. N. Ferreira, N. C. R. Ferreira, M. L. L. M. Gava,L. Zhao, and E. E. N. Macau, “The influence of timeseries distance functions on climate networks,” (2019),arXiv:1902.03298 [physics.data-an].[50] F. Hassanibesheli and R. V. Donner, Chaos , 083125(2019).[51] J. F. Donges, J. Heitzig, B. Beronov, M. Wiedermann,J. Runge, Q. Y. Feng, L. Tupikina, V. Stolbova, R. V.Donner, N. Marwan, H. A. Dijkstra, and J. Kurths,Chaos25