Unidimensional and Multidimensional Methods for Recurrence Quantification Analysis with crqa
Moreno I. Coco, Dan Mønster, Giuseppe Leonardi, Rick Dale, Sebastian Wallot
UUnidimensional and Multidimensional Methods forRecurrence Quantification Analysis with crqa
Moreno I. Coco a,b , Dan Mønster c , Giuseppe Leonardi d , Rick Dale e , SebastianWallot f a School of Psychology, University of East London, Water Lane, E15 4LZ, London, UK b CICPSI, Faculdade de Psicologia, Universidade de Lisboa, Alameda da Universidade,1649-013 Lisboa, Portugal c School of Business and Social Sciences, Aarhus University, DK-8000, Aarhus C, Denmark d Institute of Psychology, University of Economics and Human Sciences in Warsaw, 01-043,Warsaw, Poland e Department of Communication, University of California, Los Angeles, CA 90005, LosAngeles, USA f Department of Language and Literature, Max Planck Institute for Empirical Aesthetics,GR-60322, Frankfurt a.M., Germany
Abstract
Recurrence quantification analysis is a widely used method for characterizingpatterns in time series. This article presents a comprehensive survey for conductinga wide range of recurrence-based analyses to quantify the dynamical structure ofsingle and multivariate time series, and to capture coupling properties underlyingleader-follower relationships. The basics of recurrence quantification analysis (RQA)and all its variants are formally introduced step-by-step from the simplest auto-recurrence to the most advanced multivariate case. Importantly, we show howsuch RQA methods can be deployed under a single computational framework in Rusing a substantially renewed version our crqa
Keywords: recurrence quantification analysis, unidimensional andmultidimensional time series, non-linear dynamics, R-package ∗ Corresponding author
Email address: [email protected] (Moreno I. Coco) a r X i v : . [ phy s i c s . d a t a - a n ] M a y . Introduction In the current article, we present the updated 2.0 version of the R package crqa to perform many variants of recurrence-based analyses [1] including somevery recent developments for the treatment of multivariate and categoricaldata. Recurrence-based techniques allow the quantification of temporalstructure and generalized autocorrelation properties of individual time series[2, 3], the quantification of bivariate relationships and coupling between twotime series [4, 5], as well as the quantification of multidimensional dynamicsof multivariate time series [6, 7]. Recurrence-based techniques originatefrom the description and analysis of dynamical systems [8, 5], and havebeen widely applied to data from physics [9, 10, 11, 12, 13], physiology[14, 15, 16, 17, 18, 19], and psychology [20, 21, 22, 23, 24, 25], to name a fewfields.The real success of recurrence-based analyses has revolved around theirpower of capturing the dynamics of complex and non-stationary time seriesdata, and of time series exhibiting qualitatively different patterns along theirtemporal evolution [8]. This is because recurrence-based analyses are model-free techniques that make few assumptions and hence are well suited forthe analysis of complex systems. Moreover, recurrence-based analyses areversatile and can be applied to interval-scale data as well as nominal data,continuously sampled data and inter-event data alike [26, 4].The new version of the crqa package features the integration of majordevelopments in recurrence analysis, such as its extension to multidimensionaldata, as well as a key simplification of its design and a marked improvementof the underlying computational procedures. It includes useful new functions,including a tool for mining the parameter settings for continuous data, andpiece-wise computation of recurrence plots to mitigate computational cost oflong time series.The remainder of this article is divided into two broad sections. In the firstsection, we provide a concise introduction to the core concepts of recurrenceanalysis from the simplest case of auto-recurrence of a unidimensional timeseries, to the most complex case of multidimensional cross-recurrence, whichis now integrated in the new version of the crqa package. In the secondsection, we showcase example applications of the different analysis methodsusing empirical data, hence providing a hands-on tutorial for how to use thedifferent functions of the package.
2. Methodological background
In section 2.1, we briefly introduce the framework of recurrence quantifi-cation analysis (RQA hereafter) with the simplest case of a unidimensional2ime series. Then, in section 2.2, we discuss how RQA can be applied totwo different unidimensional times series. Finally, in section 2.3 and 2.4 weexplain how RQA methods can be extended to multidimensional data.
The concept of recurrence is at the heart of all recurrence-based analyses[5, 27], which mostly apply to time series or sequenced data (but see [28]). Aswe will see, recurrences are often defined in terms of phase space coordinates,not directly in terms of the values of a single time series, but the concept iseasily demonstrated using this case—a single time series or sequence—as astarting point. Loosely speaking, a recurrence is the repetition of a value in asequence of data points. More precisely, recurrence in a time series x , with n data points x , x , . . . , x n is defined as: R ij = (cid:26) x i = x j x i (cid:54) = x j i, j = 1 , , . . . n (1)The above equation defines the elements R ij of the recurrence matrix R in terms of identical repetitions. Such a definition is useful for a nominalsequence, where there are categorical elements that are either identical or not,and so no meaningful ‘distance’ norm among categories can be defined (seesection 3.3.1 for an application to text data). In order to define recurrencesfor continuous data, we need to establish a threshold parameter (or radiusparameter) ε , which provides the width of a tolerance band in the chosendistance norm within which similar, but not identical values in a time seriesare counted as recurrent: R ij = (cid:26) | x i − x j | ≤ ε | x i − x j | > ε i, j = 1 , , . . . n (2)Setting a threshold is necessary in most cases for empirical data, becausesuch data feature intrinsic fluctuations as well as measurement error [8]. Byusing recurrences as defined above, we can convert any unidimensional timeseries x into a recurrence plot (RP), which is a two-dimensional portrait ofits dynamics expressed through its recurrence characteristics.Figure 1 shows some examples of time series of various complexity (i.e.,a sinusoidal, a chaotic attractor and white noise) and their associated RPs,given some value for the threshold parameter ε .If a time series x constitutes the one-dimensional measurement of anunderlying multidimensional system, and the dynamics of the underlyingdimensions are co-dependent, then these underlying dimensions can be re-covered via the method of time-delayed embedding from the unidimensional3 igure 1: The three time series in the first row (A) — a periodical sine wave, one of thedimensions of the chaotic Lorenz attractor and a white noise signal — were subjected torecurrence analysis without dimensional embedding, that is the recurrence plots depictrecurrences based on the values of the one-dimensional time series (B). The third row (C)shows their associated recurrence plots with dimensional embedding (the sine wave wasembedded in 2 dimensions, the Lorenz attractor and the white noise signal in 3 dimensions). x is delayed (or lagged)by a certain number of data points, τ , and the number of times such delaysare applied to x depend on an embedding dimension parameter m . Thetime-shifted copies of x , x τ , x τ ,. . . , x ( m − τ , can now be integrated into a sin-gle m -dimensional phase space, which shows the recovered multidimensionaldynamics behind the measured unidimensional time series (Figure 2).If the time series data comes from a multidimensional system, embeddingthe data into a higher dimensional phase space before the computation of anRP will improve the quantification of the dynamics of the systems from whichthat time series was recorded [8]. Now, with embedded data, recurrences aredefined not in terms of the individual values of the original, unidimensionaltime series x , but in terms of coordinates in m -dimensional phase space. It ispossible to construct a total of N = n − ( m − τ points in phase space, withthe following coordinates: X = ( x , x τ , x τ , . . . , x m − τ ) X = ( x , x τ , x τ , . . . , x m − τ ) ... (3) X k = ( x k , x k + τ , x k +2 τ , . . . , x k +( m − τ ) ... X N = ( x N , x N + τ , x N +2 τ , . . . , x n ) In terms of points in the m -dimensional phase space, the elements of therecurrence matrix are then given by R ij = (cid:26) || X i − X j || ≤ ε || X i − X j || > ε i, j = 1 , , . . . N (4)where || · || is a distance norm in the m -dimensional phase space and ε isthe threshold used to determine whether the points are close enough in phasespace to be considered recurring or not. Examples of embedding and theresulting RPs are shown in Figure 1, lower panel, and in Figure 2, panel D.Note, however, that in the process of phase space reconstruction, the numberof coordinates available in phase space is less than the number of data pointsin the original time series, the difference being equal to n − N = τ ( m − .Moreover, the resulting phase space portraits do not exactly reflect the trueunderlying multidimensional dynamics, but are isomorphic to them [31]. Forcontinuous time series (i.e., not categorical), the delay τ , the embeddingdimension m and the radius ε are usually unknown, and have to be estimated5 igure 2: Graphical exemplification of the embedding method. A unidimensional timeseries (B)— here represented with time running along the y-axis—is embedded in a two-dimensional phase space (A) by choosing as the second dimension (C) of every point inpanel B the τ -lagged value (here τ = 3) of the very same original unidimensional timeseries (B). After selecting an appropriate value of the threshold (or radius) parameter(here ε = 0.25) recurrence analysis can be applied to the trajectory in the two-dimensionalphase space (A), and a recurrence plot is generated (D). In the recurrence plot (D), therecurrence points around the coordinates of ( X ; X τ ) are highlighted in red color. τ can be obtained by examining the average mutual informationfunction (AMI) of x [32], and the false-nearest-neighbor function (FNN) of x ,given some value for τ , can in turn be used to determine m [33]. The radius ε is then chosen to achieve a desired proportion of recurrence points that isexpected given the type of data at hand ([34]; see [7] for practical informationon parameter estimation). The crqa package implements functions to carryout parameter estimation semi-automatically.RPs are a very useful visualisation to qualitatively explore the dynamicsof a time series. However, their main advantage is that they provide the basisto quantify the dynamics of a time series based on the patterns of recurrencepoints found in the plot [3]. There are various measures that can be computedfrom RPs. Here, we briefly describe a selection of measures, including themost common ones, that are implemented in the new version of the crqa package.The most basic measure is the recurrence rate (RR) or percentage recur-rence, which is defined as the sum of all recurrence points in an RP divided bythe area of that RP. RR provides a measure for how many individual valuesof a time series—or its phase space coordinates—recur over time: RR = 1 N N (cid:88) i,j =1 R ij (5)All other measures characterise dynamics by exploiting the patterns ofrecurrences along the vertical and diagonal structures of the RP (see Table 1for a concise summary of the measures). In particular, measures based ondiagonal-line structures reflect repetitions of the trajectories of the time series,whereas measures based on vertical-line structures focus on the states duringwhich a time series slows down its dynamics. The entropy of the time seriescan also be computed on the basis of diagonal and vertical line structures ofthe RP. In version 2.0 of the crqa package, we include a novel entropy measure,which is based on the distribution of the areas of the rectangular structuresin an RP generated from categorical time series. This measure provides amore accurate estimation of entropy over the classic diagonal-line entropyfor categorical time series that predominantly evolve in terms of changes ofstates [35]. Until now, we focused on quantifying the dynamics of a system by wayof recurrences of a single time series. But, the concept of recurrence can beextended to that of cross-recurrence, which extends the univariate recurrenceanalysis to a bivariate analysis technique that allows quantification of the7emporal coupling properties or similarity of two time series [4, 8]. In otherwords, cross-recurrence is the recurrence of a value in a time series x , withdata points x , x , . . . , x n with the values of a time series y , with data points y , y , . . . , y n . Formally: CR ij = (cid:26) | x i − y j | ≤ ε | x i − y j | > ε i, j = 1 , , . . . n (6)As for Equation 2, a threshold parameter ε is applied to identify similar,but not necessarily identical, values that are recurrent across the two timeseries. As in the univariate case, this parameter can be set to values closer to 0,forcing cross-recurrences to be identical values, such as between two nominal(categorical) sequences. Also in the case of cross-recurrence analysis, the twotime series x and y can be embedded before computing the cross-recurrenceplot (CRP): CR ij = (cid:26) || X i − Y j || ≤ ε || X i − Y j || > ε i, j = 1 , , . . . N (7)Commonly it is expected that x and y have the same number of data points,and that the delay and embedding dimension parameters, τ and m , have to bethe same too (see [7, 36] for practical aspects of parameter estimation) . Therecurrence measures obtained from cross-recurrence quantification analysisare calculated in the same way as for the univariate recurrence quantificationanalysis (see Table 1). However, the values for cross-recurrence now reflectthe coupled dynamics of the two time series [37], rather than the dynamics ofan individual time series in univariate RQA. There is a key difference betweenRPs and CRPs. RPs always have recurrence points all along the main diagonalline of the plot (so-called line of identity, LOI), because a time series is bydefinition recurrent with itself at lag 0, which is what recurrences along themain diagonal reflect (i.e., x i = x j when i = j ). This is not necessarily thecase for CRPs as two time series are not necessarily synchronized ( x i neednot be the same as y j when i = j ). If the time series are not synchronized,cross-recurrences around the main diagonal are absent or sparse. The presenceof a full LOI in a CRP implies that the dynamics of the two time series areidentical, or that they have a strict linear relationship to each other. It is possible for x and y to be different lengths, and so to produce a rectangularCRP, but this is very rare in practice and it introduces some complications in computingsynchrony measures .3. Multidimensional recurrence quantification analysis (MdRQA) Multidimensional recurrence quantification analysis (MdRQA) is one ofthe recent extensions of RQA [6] that is now also available in the crqa package.MdRQA allows the analysis of multidimensional time series z with samples z , z , . . . , z n , where each point (sample) in z has d dimensions: z k = ( z k, , z k, , . . . , z k,d ) (8)Here, recurrences are defined on the d -dimensional coordinate space madeup of the points of z : R ij = (cid:26) || z i − z j || ≤ ε || z i − z j || > ε i, j = 1 , , . . . n (9)Like their unidimensional counterparts, multidimensional time series canbe embedded into a higher dimensional space. The logic of estimating thedelay and embedding dimension parameters, τ and m , are the same as withunivariate RQA [38, 36]. However, if one has multivariate time series andwants to estimate embedding parameters for those time series, one shoulduse the multivariate embedding functions which are also provided with thenew version of the crqa package, because they provide superior estimates ofembedding parameters for multidimensional time series [7].Depending on the underlying data, MdRQA can, in principle, be used fortwo different purposes. On the one hand, MdRQA can be used to quantify amultidimensional construct, such as physiological arousal, by simultaneouslylooking at its different measurable dimensions (e.g., heart rate, breathing andbody temperature). This would be the multivariate version of how univariateRQA quantifies the dynamics of a single unidimensional time series (e.g.,breathing alone). On the other hand, MdRQA can be used to examine theshared dynamics of multiple individual time series, such as, for example, threeelectrodermal signals measured from three members of a team performing acollaborative task. Here, MdRQA variables would be interpreted as capturinghigher-order inter-correlative properties between the three signals at the levelof the group [6].In either case, one has to make a decision about whether to normalizethe different dimensions of the time series or not. If each dimension of themultidimensional time series is normalized, for example z -scored, it wouldeffectively give each time series equal weight for the definition of recurrence.In particular, if one does not know how the different time series interact, or ifthey are measured on different scales without regard to their potential effectson one another, then normalization is recommended, otherwise the risk is toassign a greater weight to the time series bearing greater variance. If one9s certain that the values of each dimension of the multidimensional timeseries are already properly scaled with regard to each other, such as whensimultaneously analyzing the three dimensions of the Lorenz system [39], thenone needs not—and perhaps should not—normalize the dimensions. Multidimensional cross-recurrence quantification analysis (MdCRQA) ex-tends MdRQA in the same way that CRQA extends RQA. Effectively, Md-CRQA allows for the computation of cross-recurrences between two multi-dimensional time series x and y [40], where cross-recurrences are definedbetween the two d -dimensional coordinate spaces between the points of x and y : x i = ( z i, , z i, , . . . , z i,d ) (10) y j = ( z j, , z j, , . . . , z j,d ) (11) CR ij = (cid:26) || x i − y j || ≤ ε || x i − y j || > ε i, j = 1 , , . . . n (12)It is important that the different dimensions of the two multivariate timeseries enter the analysis in the same order. For the example of physiologicalarousal, if heart rate is the first dimension in x , breathing is the seconddimension in x , and body temperature is the third dimension in x , then heartrate also needs to be the first dimension in y , breathing needs to be thesecond dimension in y , and accordingly body temperature needs to be thethird dimension in y . Otherwise, the resulting MdCRQA measure will not beinterpretable.Unidimensional and multidimensional cross-recurrence analysis can alsobe performed in a time-dependent manner. This so-called windowed cross-recurrence analysis is conceptually very similar to windowed cross-correlationanalysis as in Boker et al. [41] and it can be used to track how cross-recurrencechanges over the time course. To that end, one simply partitions the timeseries of interest into a number of overlapping or non-overlapping sub-seriesand calculates CRPs for each of the sub-series. For each CRP, i.e., eachsub-series of the original time series, the recurrence measures are calculatedwhich allows tracking changes in cross-recurrence over time. Finally, from cross-recurrence plots of either two unidimensional timeseries (i.e., CRQA) or two multidimensional time series (i.e., MdCRQA), it is10ossible to extract the diagonal cross-recurrence profiles (DCRPs), and useit to capture leader-follower-relationships [26, 8]. To that end, one has tospecify a window size w for the number of lags on the recurrence plot thatone wants to investigate. For example, a window-size of 10 data-points, i.e., w = 10 would span recurrences of ±
10 diagonals from the LOI. Specifically,this procedure determines the cross-recurrence rates of each diagonal, CR w ,by summing up all cross-recurrence points that fall along such diagonal, anddivide them by its length: CR w = 1 N − w N − w (cid:88) j − i = w R i,j (13)The DCRP permits quantification of the recurrence rate over differentrelative lags between two time series. If the peak of cross-recurrences fall alongthe central diagonal, the line of synchronization (LOS), then this suggestsstrong coupling at lag 0 between the two time series. If the peak of cross-recurrences falls instead on one of the diagonals off the LOS, it indicates thatthe dynamics of one time series follow the dynamics of the other time seriesby some lag equal to that diagonal position (Figure 3).Note, however that interpreting the lags in terms of the sampling rate of theunderlying measured time series only applies to CRPs based on unembedded(often categorical) time series. If the time series are embedded, this meansthat the observed lag is based on coordinates that are made up of severaldata points from the respective original time series, and hence introducesa degree of uncertainty with regard to the precise time interval of the lagwhen one tries to map particular recurrence points back to data points ofthe original time series. In addition, the leader-follower interpretation of thelags cannot be granted the status of a causal interpretation. For example, aparent can deliberately ‘lag’ their behavior behind that of their child. In sucha case, it would be odd to say definitively that the child’s behavior is causingthe parent’s behavior. For this reason, the DCRP has to be interpretedwith caution and is best treated as a general description of relative temporalrelationships.
3. Package design
The crqa is available from the Comprehensive R Archive Network (
CRAN )at https://cran.r-project.org/web/packages/crqa/index.html. The software iswritten in R with the exception of the spdiags function, which contains asection written in Fortran . The main function is crqa , which takes its namefrom the package. This function performs all types of recurrence quantification11 igure 3: Illustration of DCRP. Time series of two of the three dimensions from the Lorenzsystem, one in red, the other one in black (A). After computation of their CRP, one candefine their diagonal cross-recurrence profile (DCRP) in order to quantify their time-laggedcoupling properties (B). As can be seen, the two time series are most strongly coupled ataround the lags of 5 to 9, where a peak in the DCRP can be observed. analyses discussed in the previous section, and it returns the actual recurrenceplot, along with the measures listed in Table 1 extracted from it. Here is howto call this primary function, with arguments and defaults: crqa ( ts1 , ts2 , delay , embed , rescale = 0 , radius ,normalize = 0 , mindiagline = 2 , minvertline = 2 ,tw = 0 , whiteline = FALSE , recpt = NULL ,side = ’ upper ’ , method = ’ crqa ’ , metric = ’ euclidean ’ ,datatype = ’ continuous ’)
Several arguments in this function allow the user to refine aspects of thecomputation. For example, the user can modify the metric to obtain thedistance matrix when estimating recurrence, or the settings of thresholds toaccept contiguous points as recurring (see Table 2 for the list of argumentsof crqa() with a brief explanation of each). drpfromts() can be usedto obtain the diagonal cross-recurrence profile, and this function is builtusing the crqa() at its core. In fact, most arguments stay exactly thesame, and we will illustrate the additional arguments that are specific tothis function when describing it. Likewise, the functions wincrqa() and windowdrp() are built on the main function crqa() and are used to computewindowed cross-recurrence. The package also contains functions, such as optimizeParam() , to automatically estimate the settings for the three mainparameters of RQA analyses, i.e., radius ε , delay τ and embedding dimension m , for continuous measures. As recurrence quantification analysis is heavyon memory requirements, the package features a function ( piecewiseRQA() )which can be used to break down the analysis of long time series into smaller12nd more manageable chunks that are computationally more tractable. Finally,the package provides the user with functions to simulate data from classicexamples of dynamical systems, such as the Lorenz ( lorenzattractor() ) orcategorical series with different distributions ( simts() ). In what follows, webriefly describe the data available with the package. These data are used toillustrate the different functionalities of the crqa() package. A search of CRAN shows very few other packages offering alternativemethods to compute recurrence quantification analysis. In particular, tseri-esChaos , nonlinearTseries and RHRV . The function recurr in tseriesChaos computes a recurrence plot for a (section of) a single unidimensional timeseries,but it does not compute any derived measures from the plot characterisingthe dynamics of the system, nor does it handle cross recurrence plots orany of the many extensions to simple recurrence plots provided by our crqa package. Going a little further, the function rqa in nonlinearTseries returnskey measures from the recurrence plots (e.g., recurrence rate), as well as, aquick visualisation of the recurrence plot (available also with the function RecurrencePlot ). Finally, the function
RecurrencePlot in RHRV is simplya wrapper built on the function in nonlinearTseries with the same name, butspecifically tailored to electrocardiogram data. The functions available in non-linearTseries provide very basic functionality in terms of the range of metricsavailable and optimization routines. They do not allow the user to examinediagonal structures for leader-follower analyses, nor explore the evolution ofrecurrence rate using windowed methods. All such features are integratedinto crqa , which, to the best of our knowledge, is the most comprehensivestatistical package to perform recurrence quantification analysis in R.
Different types of categorical and continuous time series, both unidimen-sional and multidimensional, are available with the package. The command load(crqa) will load the data into the R workspace. In particular, we includea nursery rhyme “The wheels on the bus” by Verna Hills to illustrate the mostbasic recurrence quantification analysis. This text is a vector of 120 strings(i.e., the words of the song), and as it is extremely simple and highly repetitive,it makes it a very good example to illustrate the core concept of recurrence.Then, we move on to cross-recurrence quantification analysis and include inthe data object of the package eye-tracking data from the study by [42]. Inthis study, a narrator describes the characters of a TV series (Friends) to alistener, who will have to later answer some comprehension questions aboutthem, while their eye-movement are co-registered. Here, we use a single trial of13his study, which is stored as a dataframe of 2,000 observations of six possiblescreen locations that are looked at by the narrator and the listener. These arenumerically coded from 1 to 6, representing a 2x3 visual grid ). This data willalso be used to illustrate diagonal and windowed cross-recurrence. Finally, weillustrate multidimensional cross-recurrence analysis using hand movementdata from the study by [43]. In this study, dyads were instructed to cooperatein a complex LEGO joint construction task, under different conditions, whiletheir hand movements and heart-rates were co-registered. Again, we selectonly a single trial of hand-movement from the turn-taking condition. Thedataframe comprises of 5,799 observations from two participants (P1 and P2)for the dominant (_d) and non-dominant (_n) hand. crqa package3.3.1. RQA As explained in section 2.1, RQA entails computing the auto-recurrenceof a unidimensional time series. In the context of the nursery rhyme, weexpect clear phases of recurrence to emerge because the words greatly repeat.In order to run this analysis, we use the main function crqa , and specifyin the argument method that we are running an RQA analysis ( method ="rqa" ). As the data that we use is categorical, we need to specify in theargument datatype that the nature of the data is categorical ( datatype ="categorical" ) this will automatically recode the categorical states of theseries (i.e., the words) into unique numerical integers so that recurrence can becomputed using a radius that has to be smaller than 1 (e.g., radius = 0.01 ),so that we can capture the recurrence of identical words. For categoricalRQA, the delay and embedding dimension have to be set to 1 ( delay = 1;embed = 1 . ). We also need to set the Theiler window parameter to 1 ( tw = 1), so that we can exclude the LOI from all recurrence measures. Finally,the same unidimensional time series has to be input both as ts1 and ts2 toobtain its auto-recurrence. res ← crqa ( text , text , delay = 1 , embed = 1 , rescale ,radius = 0 .01 , normalize , mindiagline ,minvertline ,tw = 1 , whiteline , recpt , side , method = " rqa " ,metric = " euclidean " , datatype = " categorical " ) There are two more states, 10 and 11, to indicate when the listener or the narratorblinked or looked outside of the screen, or to identify possible blinks. They are coded witha different number so that these two states will not recur when radius is set near 0. If embedding dimension is set to higher values, this becomes equivalent to doingrecurrence on n -grams, where m = n . In fact this interpretation of categorical recurrencecreates bridges to traditional natural language processing, summarized in [44]
14e can use the plotting function plotRP to visualise the resulting recur-rence plot. This function provides some basic arguments to change the sizeof the points in the plot ( pcex ), the color ( cols ), or their type ( pch ) whichare taken verbatim from the generic plot function. RP ← res $ RPparC ← list ( unit = 10 , labelx = " Time " , labely = " Time " ,cols = " black " , pcex = .5 , pch = 15 , las = 0 ,labax = seq (0 , nrow ( RP ) , 10) ,labay = seq (0 , nrow ( RP ) , 10) )plotRP ( RP , parC ) In order to get a closer understanding of recurrences, we zoom into asegment of the text, re-run the crqa() function and visualise it (Figure 4).We add the labels of the axes (x,y), print the words vertically using the las argument, and decrease the temporal unit argument to print each individualword on the axes. text_zoom ← text [81:110]ans_zoom ← crqa ( text_zoom , text_zoom , delay , embed ,rescale , radius , normalize , mindiagline ,minvertline , tw , whiteline , recpt , side ,method , metric , datatype )RP ← ans_zoom $ RPparC $ labay ← parC $ labax ← text_zoomparC $ las ←
2; parC $ unit ← ← parC $ labely ← " Words "plotRP ( RP , parC ) As it can be clearly seen in Figure 4, the bulk of recurrence in that portionis driven by the repeated use of the single word wah . When looking at afew measures associated with the RP, we observe an overall determinism of85.4%, which implies that the system is fairly repetitive, an average diagonalline length of 3.88, which implies that on average there are sequences of fourwords that repeat, and a maximum diagonal length of 9, which means thatthe longest sequence repeating is made of 9 words. Some measures such asdeterminism or the average diagonal line length depend on the setting of theargument mindiagline (equivalent to l min in Table 1). The default value ofthis argument is 2, as two contiguous points form a line, but it can depend onthe type of data (e.g., words vs. eye-movement), or the sample rate at whichit is acquired (55 Hz vs. 1,000 Hz). For example, if we have acquired data at1,000 Hz, we would practically have one data-point every 2 ms. This meansthat if we use a default mindiagline of 2, we would be considering as lines,any states that contiguously repeat over a 4 ms window. This value wouldcertainly be unrealistic for some type of responses that unfolds over longer15 igure 4: Recurrence quantification of the nursery rhyme ‘wheels on the bus’. On theleft panel, we show the full recurrence and the red box indicates the region of RP that iszoomed in on the right panel. period of time (e.g., an eye-movement fixation lasts for an average of 200 ms). crqa() As already explained in section 2.2, cross-recurrence is an extension ofauto-recurrence to two different unidimensional time series. In order to run across-recurrence analysis with the crqa package, we simply need to change the method argument to method = "crqa" . Here, we illustrate its use throughtwo time series of eye-movement data explained in section 3.2. Also, we inputtwo different time series of eye-movement data, narrator and listener ,rather than just one. An optional argument that is available in the crqa function is the side (upper, lower or both) of the recurrence plot on whichrecurrence measures are computed. This may be useful, for example, forresearchers interested in leader-follower dynamics (see Figure 5, left panel,for the visualisation of the cross recurrence plot). drpfromts
In section 2.2, we explained what diagonal cross-recurrence is, and how itcan be used. In the crqa package, this measure is computed by the function drpfromts , which utilises the same arguments of the main crqa functionplus an additional argument, windowsize , to define the number of lags (ordiagonals) around the line of synchronization (LOS) of the CRP over whichrecurrence rate is computed. In the example visualised in Figure 5, centerpanel, we have chosen a window of 100 lags, which spans about ± secondsaround the LOS. We can clearly see that the peak recurrence is shifted by16 igure 5: Cross-recurrence using a single trial of eye-movement data measured on alistener and a narrator. On the left panel, we visualise the full cross-recurrence plot.The transparent red band represents the lags around the line of coincidence that havebeen used to calculate the diagonal cross-recurrence plot (center panel), whereas the bluetransparent squares within it are the overlapping windows used to compute the windowedcross-recurrence profile (right panel). ≈ second from the LOS, i.e., lag 0. This reflects the time taken by thelistener to look at the same panel the narrator was looking at—namely, about1 second for a listener to “catch up” to the speaker. res ← drpfromts ( narrator , listener , windowsize = 100 ,radius = 0 .001 , delay = 1 , embed = 1 ,rescale = 0 , normalize = 0 ,mindiagline = 2 , minvertline = 2 , tw = 0 ,whiteline = F , recpt = F , side = " both " ,method = " crqa " , metric = " euclidean " ,datatype = " continuous " ) windowdrp Windowed cross-recurrence captures the evolution of recurrence rate overtime. In the context of the eye-movement data, this measure reflects howconsistently listener and narrator are looking at the same scene location at anygiven point in time. More importantly, the windowed methodology revealshow recurrence changes over time (refer to section 2.5 for more details). Weuse the function windowdrp to compute this measure, which again shares thesame arguments of crqa , plus three more that are specific to it. In particular,we have to set: (a) the size of the window that slides over the time-course,e.g., windowsize = 100 , (b) the step that we want this window to move, e.g., windowstep = 20 , and (c) the number of lags within the window of interest Note, that the number of lags cannot be greater than the size of the window. lagwidth = 50 . In Figure 5,right panel, we observe that recurrence rate grows over time, which meansthat the narrator and the listener tend to look more and more at the samepanels as the trial progresses. res ← windowdrp ( narrator , listener , windowstep = 20 ,windowsize = 100 , lagwidth = 50 ,radius = 0 .001 , delay = 1 , embed = 1 ,rescale = 0 , normalize = 0 , mindiagline = 2 ,minvertline = 2 , tw = 0 , whiteline = F ,side = " both " , method = " crqa " ,metric = " euclidean " , datatype = " continuous ") Finally, we compute multidimensional cross-recurrence by setting the method argument, i.e., method = "mdcrqa" . Note, in order to computemultidimensional recurrence, the user needs to provide the same dataframe asinput to both ts1 and ts2 . This method is also available for drpfromts and windowdrp . Applied to the hand movement data, we first restructure the dataof the two participants into independent sets. We re-use the same parametersettings of [43] to compute MdCRQA, and leave all other arguments withtheir default values. P1 ← cbind ( handset $ P1_TT_d , handset $ P1_TT_n )P2 ← cbind ( handset $ P2_TT_d , handset $ P2_TT_n )res ← crqa ( P1 , P2 , delay = 5 , embed = 2 ,rescale = 0 , radius = 0 .1 , normalize = 0 ,mindiagline = 10 , minvertline = 10 , tw = 0 ,whiteline , recpt , side , method = " mdcrqa " ,metric , datatype ) piecewiseRQA Often, researchers are interested in time series which contain severalthousands of observations. Sometimes the dimensionality of these time seriescan be reduced without losing too much information, such as by down-sampling. This strategy may not always be possible or serve the researcher’spurpose. Recurrence quantification analysis can require more RAM thanis available in standard laptops or personal workstations, hence making itnearly impossible to run. In the new version of the crqa package, we providethe user with the piecewiseRQA function, which can be used to compute alldifferent variants of recurrence quantification analysis described above on longtime series. Conceptually, this function divides the time series into blocks, it18btains a recurrence plot for each individual block, and then fills the originalrecurrence plot in with all such sub-blocks, before computing the measures. For example, if we are handling a time series of 10,000 observations,we could divide it into 10 blocks of 1,000 observation each. piecewiseRQA has exactly the same arguments that we have already encountered in themain crqa function, but has an additional two that are used to control thesize of the block, e.g., blockSize = 100 , and the argument typeRQA , whichcan take two options, either full or diagonal . If the value for typeRQA is diagonal only the diagonal cross-recurrence will be computed; if full , therecurrence measures will be obtained out of the full plot . In Figure 6, wevisualise the computational speed (left panel) and memory demand (rightpanel) on simulated time series of sinusoids of increasing size, and comparewhat happens if we run the piecewiseRQA , also with blocks of increasingsize, as compared to running the main crqa function. We can clearly seethat for time series of increasing length, memory demands are kept lower bythe piecewiseRQA function as compared to the crqa . However, we can alsosee that there is a wide variance for blocks of different sizes. Therefore, itmay be wise to explore the block sizes to find the one that can optimize thecomputational performance over a single trial before running the piecewiserecurrence analysis on an entire dataset. optimizeParam The last function we showcase is optimizeParam , which helps the userexploring the space of values for the parameters of delay, embedding dimensionand radius to compute recurrence quantification on continuous valued timeseries. In particular, optimizeParam first estimates the average mutualinformation (AMI) of either the unidimensional or multidimensional timeseries, and chooses the value that minimizes it. Then, it takes such a delayvalue, and evaluates the embedding dimensions maximizing the false-nearestneighbors (FNN) . As a last step, optimizeParam applies the values of delayand embedding dimension obtained, to find a radius which returns a recurrencerate within a minimum and maximum value established by the user. Weapply optimizeParam() to simulated sinusoids to show the unidimensionalcase, and to the hand movement data of [43] to show the multidimensional This function is similar to the crp_big function in the long-standing CRP-toolbox forMATLAB by Marwan and colleagues: http://tocsy.pik-potsdam.de. Currently, the windowed cross-recurrence is not implemented to work with the piecewiseRQA . Users can also access the functions to compute delay (
MdDelay ) and embeddingdimensions (
MdFnn ) independently of optimizeParam . igure 6: Evaluating the speed (time in seconds, left panel) and memory (peak RAMin MB, right panel) performance of crqa() and piecewiseRQA() for simulated data ofincreasing size (from 3000 to 7000 data points). We compare the use of blocks of differentsizes (from 1000 to 6500 in increments of 500, coded using color and point type) with thecase of running crqa() on the entire sequence of data points. case.In order to set up the estimation of the delay and embedding dimensionsfor unidimensional time series, the user has to decide how to choose thevalue of average mutual information (i.e., typeami = mindip , the lag atwhich minimal information is observed, or typeami = maxlag , the maximumlag at which minimal information is observed) and the relative percentageof information gained in FNN, relative to the first embedding dimension,when higher embeddings are considered (i.e., fnnpercent ). Then, as crqa is integrated in the optimizeParam to estimate the radius, most of thearguments are the same (e.g., mindiagline or tw ), except the number ofvalues of that are considered (i.e., radiusspan = 100 ). ts1 ← seq (0 .1 , 200 , .1 )ts1 ← sin ( ts1 ) + linspace (0 , 1 , length ( ts1 ) )ts2 ← ts1par ← list ( method = " rqa " , metric = " euclidean " ,maxlag = 20 , radiusspan = 100 , normalize = 0 ,rescale = 4 , mindiagline = 10 , minvertline = 10 ,tw = 0 , whiteline = FALSE , recpt = FALSE ,side = " both " , datatype = " continuous " ,fnnpercent = 10 , typeami = " mindip " )results ← optimizeParam ( ts1 , ts2 , par ,min.rec = 2 , max.rec = 5)print ( unlist ( results ) )radius emddim delay .17 2 18 For multidimensional series, the user needs to specify the right RQAmethod (i.e., method = "mdcrqa" ). Then, for the estimation of the delay viaAMI: (1) nbins the number of breaks used to define the bins within whichthe two-dimensional histogram (or frequency distribution) of the original anddelayed time series are computed and (2) the criterion to select the delay( firstBelow to use the lowest delay at which the AMI function drops belowthe value set by the threshold argument, and localMin to use the positionof the first local AMI minimum). The estimation of the embedding dimensionsinstead needs the following arguments: (1) maxEmb , which is the maximumnumber of embedding dimensions considered, (2) noSamples , which is thenumber of randomly drawn coordinates from phase space used to estimatethe percentage of false-nearest neighbors, (3)
Rtol , which is the first distancecriterion for separating false neighbors, and (4)
Atol , which is the seconddistance criterion for separating false neighbors. The radius is estimated asbefore. par $ method = " mdcrqa " ; par $ nbins = 50par $ criterion = " firstBelow " ; par $ threshold = 1 .6par $ maxEmb = 20; par $ numSamples = 500; par $ Rtol = 10;par $ Atol = 2results ← optimizeParam ( P1 , P2 , par ,min.rec = 2 , max.rec = 5)print ( unlist ( results ) )radius emddim delay0 .032 11 2
4. Conclusion
This paper describes recurrence quantification analysis, a statistical methodto characterise the nonlinear dynamics of a system. It has received a grow-ing interest by researchers across very different fields from physiology topsychology because of its flexibility, ease of application, and explanatorypower. In particular, we explain recurrence analysis from the simplest caseof auto-recurrence of a unidimensional time series to the most complex caseof multidimensional cross-recurrence. More importantly, we presented a sig-nificantly updated version of the crqa to perform all different variants ofrecurrence analysis described in the theoretical section of this manuscript. Weshowcased the different functions available in crqa with real data of categoricaland continuous nature, and illustrated how starting parameters for continuousdata can be obtained (i.e., radius, embedding dimension and delay) as well as21andling long time series in a memory efficient way using additional functionsavailable in crqa .It is useful to end on some observations regarding the broader relevanceof the package presented here. The RQA methodology and the updated crqa tap into a number of evolving problems in data analysis across variousdisciplines. For example, there is a drive to improve measures and modelsof multimodal data [45, 46, 47], including multi-person measures [48, 49,50, 51, 43]. RQA and its multidimensional counterpart implemented in ourpackage constitute an extraordinarily expansive analysis tool for exploringvaried kinds of complex and multidimensional data. In addition, a demandfor more dynamic quantitative analyses has now also penetrated into thesocial sciences [52, 53, 54, 55]. rqa is designed to be a comprehensive analysispackage for studying the dynamics of diverse systems, especially systems thatexhibit high degrees of interdependence, and that show signatures in theirdynamics that are critical for understanding them.
SW acknowledges support from the Deutsche Forschungsgemeinschaft(DFG, German Research Foundation)–WA 3538/4-1. MIC acknowledgessupport from the Fundaçâo para a Ciência e Tecnologia under grant agreementPTDC/PSI-ESP/30958/2017 22 eferences [1] M. I. Coco, R. Dale, Cross-recurrence quantification analysis of categori-cal and continuous time series: an R package, Frontiers in psychology 5(2014) 510. doi: .[2] C. L. Webber, J. P. Zbilut, Dynamical assessment of physiologicalsystems and states using recurrence plot strategies, Journal of AppliedPhysiology 76 (1994) 965–973. doi: .[3] J. P. Zbilut, C. L. Webber, Embeddings and delays as derived fromquantification of recurrence plots, Physics Letters A 171 (1992) 199 –203. doi: .[4] J. P. Zbilut, A. Giuliani, C. L. Webber, Detecting deterministic signalsin exceptionally noisy environments using cross-recurrence quantification,Physics Letters A 246 (1998) 122 – 128. doi: .[5] N. Marwan, J. Kurths, Nonlinear analysis of bivariate data with crossrecurrence plots, Physics Letters A 302 (2002) 299–307.[6] S. Wallot, A. Roepstorff, D. Mønster, Multidimensional recurrencequantification analysis (mdrqa) for the analysis of multidimensional time-series: A software implementation in MATLAB and its application togroup-level data in joint action, Frontiers in psychology 7 (2016) 1835.doi: .[7] S. Wallot, D. Mønster, Calculation of average mutual information(ami) and false-nearest neighbors (fnn) for the estimation of embeddingparameters of multidimensional time series in
MATLAB , Frontiers inpsychology 9 (2018). doi: .[8] N. Marwan, M. C. Romano, M. Thiel, J. Kurths, Recurrence plots forthe analysis of complex systems, Physics Reports 438 (2007) 237 – 329.doi: .[9] P. Alex, S. Arumugam, K. Jayaprakash, K. Suraj, Order–chaos–order–chaos transition and evolution of multiple anodic double layers in glowdischarge plasma, Results in Physics 5 (2015) 235–240.[10] B. Ambrożkiewicz, Y. Guo, G. Litak, P. Wolszczak, Dynamical responseof a planetary gear system with faults using recurrence statistics, in:Topics in Nonlinear Mechanics and Physics, Springer, 2019, pp. 177–185.2311] R. Donner, M. Thiel, Scale-resolved phase coherence analysis of hemi-spheric sunspot activity: a new look at the north-south asymmetry,Astronomy & Astrophysics 475 (2007) L33–L36.[12] V. Hilarov, Detection of the fracture zone by the method of recurrenceplot, Physics of the Solid State 59 (2017) 2401–2406.[13] N. Zolotova, D. Ponyavin, Phase asynchrony of the north-south sunspotactivity, Astronomy & Astrophysics 449 (2006) L1–L4.[14] N. Marwan, N. Wessel, U. Meyerfeldt, A. Schirdewan, J. Kurths,Recurrence-plot-based Measures of Complexity and their Applicationto Heart-rate-variability Data, Physical Review E 66 (2002) 026702.doi: .[15] J. Langbein, G. Nürnberg, G. Manteuffel, Visual discrimination learningin dwarf goats and associated changes in heart rate and heart ratevariability, Physiology & behavior 82 (2004) 601–609.[16] N. Thomasson, T. J. Hoeppner, C. L. Webber Jr, J. P. Zbilut, Recurrencequantification in epileptic eegs, Physics Letters A 279 (2001) 94–101.[17] D. Mestivier, H. Dabiré, N. P. Chau, Effects of autonomic blockers onlinear and nonlinear indexes of blood pressure and heart rate in shr,American Journal of Physiology-Heart and Circulatory Physiology 281(2001) H1113–H1121.[18] D. Mønster, D. D. Håkonsson, J. K. Eskildsen, S. Wallot, Physiologicalevidence of interpersonal dynamics in a cooperative production task,Physiology & Behavior 156 (2016) 24 – 34. doi: https://doi.org/10.1016/j.physbeh.2016.01.004 .[19] L. T. Timothy, B. M. Krishna, U. Nair, Classification of mild cognitiveimpairment eeg using combined recurrence and cross recurrence quantifi-cation analysis, International Journal of Psychophysiology 120 (2017)86–95.[20] D. H. Abney, A. S. Warlaumont, A. Haussman, J. M. Ross, S. Wallot,Using nonlinear methods to quantify changes in infant limb movementsand vocalizations, Frontiers in psychology 5 (2014) 771.[21] M. I. Coco, R. Dale, F. Keller, Performance in a collaborative searchtask: The role of feedback and alignment, Topics in cognitive science 10(2018) 55–79. 2422] M. I. Coco, L. Badino, P. Cipresso, A. Chirico, E. Ferrari, G. Riva,A. Gaggioli, A. D’Ausilio, Multilevel behavioral synchronization in a jointtower-building task, IEEE Transactions on Cognitive and DevelopmentalSystems 9 (2016) 223–233.[23] K. Shockley, M.-V. Santana, C. A. Fowler, Mutual interpersonal pos-tural constraints are involved in cooperative conversation., Journal ofExperimental Psychology: Human Perception and Performance 29 (2003)326.[24] S. Wallot, J. T. Lee, D. G. Kelty-Stephen, Switching between readingtasks leads to phase-transitions in reading times in l1 and l2 readers,PloS one 14 (2019) e0211502.[25] M. Wijnants, F. Hasselman, R. Cox, A. Bosman, G. Van Orden, Aninteraction-dominant perspective on reading fluency and dyslexia, Annalsof dyslexia 62 (2012) 100–119.[26] R. Dale, A. S. Warlaumont, D. C. Richardson, Nominal cross recur-rence as a generalized lag sequential analysis for behavioral streams,International Journal of Bifurcation and Chaos 21 (2011) 1153–1161.doi: .[27] L. Trulla, A. Giuliani, J. Zbilut, C. Webber Jr, Recurrence quantificationanalysis of the logistic equation with transients, Physics Letters A 223(1996) 255–260.[28] S. Wallot, G. Leonardi, Deriving inferential statistics from recurrenceplots: A recurrence-based test of differences between sample distributionsand its comparison to the two-sample kolmogorov-smirnov test, Chaos:An Interdisciplinary Journal of Nonlinear Science 28 (2018) 085712.doi: .[29] N. H. Packard, J. P. Crutchfield, J. D. Farmer, R. S. Shaw, Geometryfrom a time series, Physical review letters 45 (1980) 712. doi: .[30] F. Takens, Detecting strange attractors in turbulence, in: D. Rand,L.-S. Young (Eds.), Dynamical Systems and Turbulence, Warwick 1980,Springer Berlin Heidelberg, Berlin, Heidelberg, 1981, pp. 366–381. doi: .[31] J. Garland, E. Bradley, J. D. Meiss, Exploring the topology of dynamicalreconstructions, Physica D: Nonlinear Phenomena 334 (2016) 49 – 59.25oi: , topology in Dynamics, DifferentialEquations, and Data.[32] A. M. Fraser, H. L. Swinney, Independent coordinates for strangeattractors from mutual information, Phys. Rev. A 33 (1986) 1134–1140.doi: .[33] M. B. Kennel, R. Brown, H. D. I. Abarbanel, Determining embedding di-mension for phase-space reconstruction using a geometrical construction,Phys. Rev. A 45 (1992) 3403–3411. doi: .[34] C. Webber Jr, J. Zbilut, Recurrence quantification analysis of nonlineardynamical systems. national science foundation, washington dc, usa.chapter 2, methods for the behavioral sciences, eds. m. riley and g.van orden, 2005. URL: https://nsf.gov/pubs/2005/nsf05057/nmbs/nmbs.jsp .[35] G. Leonardi, A method for the computation of entropy in the recurrencequantification analysis of categorical time series, Physica A: StatisticalMechanics and its Applications 512 (2018) 824 – 836. doi: .[36] S. Wallot, G. Leonardi, Analyzing multivariate dynamics using cross-recurrence quantification analysis (crqa), diagonal-cross-recurrence pro-files (dcrp), and multidimensional recurrence quantification analysis(mdrqa) – a tutorial in R , Frontiers in Psychology 9 (2018) 2232.doi: .[37] K. Shockley, M. Butwill, J. P. Zbilut, C. L. J. Webber, Cross recurrencequantification of coupled oscillators, Physics Letters A 305 (2002) 59–69.[38] S. Wallot, Recurrence quantification analysis of processes and productsof discourse: A tutorial in R , Discourse Processes 54 (2017) 382–405.doi: .[39] E. N. Lorenz, Deterministic Nonperiodic Flow, Journal of the Atmo-spheric Sciences 20 (1963) 130–141.[40] S. Wallot, Multidimensional cross-recurrence quantification analysis(mdcrqa) – a method for quantifying correlation between multivariatetime-series, Multivariate Behavioral Research 54 (2019) 173–191. doi: . 2641] S. M. Boker, J. L. Rotondo, M. Xu, K. King, Windowed cross-correlationand peak picking for the analysis of variability in the association betweenbehavioral time series., Psychological methods 7 (2002) 338.[42] D. C. Richardson, R. Dale, Looking to understand: The coupling betweenspeakers’ and listeners’ eye movements and its relationship to discoursecomprehension, Cognitive science 29 (2005) 1045–1060.[43] S. Wallot, P. Mitkidis, J. J. McGraw, A. Roepstorff, Beyond synchrony:joint action in a complex production task reveals beneficial effects ofdecreased interpersonal synchrony, PloS one 11 (2016) e0168306.[44] R. Dale, N. D. Duran, M. Coco, Dynamic natural language processingwith recurrence quantification analysis, arXiv preprint arXiv:1803.07136(2018).[45] J. Abawajy, Comprehensive analysis of big data variety landscape,International journal of parallel, emergent and distributed systems 30(2015) 5–14.[46] R. Dale, An integrative research strategy for exploring synergies innatural language performance, Ecological Psychology 27 (2015) 190–201.[47] D. Lahat, T. Adali, C. Jutten, Multimodal data fusion: an overview ofmethods, challenges, and prospects, Proceedings of the IEEE 103 (2015)1449–1477.[48] N. J. Cooke, J. C. Gorman, C. W. Myers, J. L. Duran, Interactive teamcognition, Cognitive science 37 (2013) 255–285.[49] D. López Pérez, G. Leonardi, A. Niedźwiecka, A. Radkowska,J. Rączaszek-Leonardi, P. Tomalski, Combining recurrence analysisand automatic movement extraction from video recordings to study be-havioral coupling in face-to-face parent-child interactions, Frontiers inpsychology 8 (2017) 2228.[50] L. Schilbach, B. Timmermans, V. Reddy, A. Costall, G. Bente, T. Schlicht,K. Vogeley, Toward a second-person neuroscience 1, Behavioral andbrain sciences 36 (2013) 393–414.[51] J. von Zimmermann, D. C. Richardson, Verbal synchrony and actiondynamics in large groups, Frontiers in Psychology 7 (2016) 2034.[52] A. Chemero, Radical embodied cognitive science, MIT press, 2011.2753] J. Friedenberg, Dynamical psychology: Complexity, self-organization andmind., ISCE Publishing, 2009.[54] M. Spivey, The continuity of mind, Oxford University Press, 2008.[55] L. M. Ward, Dynamical cognitive science, MIT press, 2002.28 easure Abbreviation Definition Recurrence Rate RR N N (cid:88) i,j =1 R ij Determinism DET N (cid:88) l = l min lP ( l ) (cid:44) N (cid:88) l =1 lP ( l ) Average Diagonal Line Length L N (cid:88) l = l min lP ( l ) (cid:44) N (cid:88) l = l min P ( l ) Maximum Diagonal Line Length maxL max( { l i } N l i =1 ) , N l = (cid:88) l ≥ l min P ( l ) Diagonal Line Entropy ENTR − N (cid:88) l = l min p ( l ) log p ( l ) Laminarity LAM N (cid:88) v = v min vP ( v ) (cid:44) N (cid:88) v =1 vP ( v ) Trapping Time TT N (cid:88) v = v min vP ( v ) (cid:44) N (cid:88) v = v min P ( v ) Categorical Area-based Entropy catH − N a (cid:88) a> p ( a ) log p ( a ) Table 1: Summary and definition of RQA measures. Here, l is some diagonal line length onthe recurrence plot, i.e. the number of diagonally adjacent recurrence points; P ( l ) is thehistogram or frequency distribution of such diagonal line lengths; p ( l ) is probability of somediagonal line length; v is some vertical line length on the recurrence plot, i.e., the number ofvertical adjacent recurrence points; P ( v ) is the histogram or frequency distribution of suchdiagonal line lengths; l min and v min are the minimum diagonal and vertical line lengthsincluded in the measures ( ≥ ); a is the value of the area of a rectangular recurrence blockgenerated in a categorical recurrence analysis; p ( a ) is the probability of the recurrenceblocks of area a . rgument Description Usage Notes ts1 and ts2 The input time series, either unidimen-sional or multidimensional time series If auto-recurrence is used (see method = rqa , then t1and t2 should be the same time series. delay
A constant indicating the number of timepoints used to lag the time series It corresponds to the τ of the equations. embed A constant indicating the number of em-bedding dimensions applied to the timeseries It corresponds to the m of the equations. rescale Whether the distance matrix on whichrecurrence is evaluated should be rescaledand how If rescale = 0 , keep the distance matrix as is; if rescale = 1 , rescale the distance matrix to its mean;if rescale = 2 , rescale to distance matrix to its maxi-mum value. radius
A constant used to decide whether the dis-tance between two points is small enoughto be considered recurrent For categorical time series, the radius needs to be setat values smaller than 1. For continuous time series,the value of the radius needs to tailored to the type ofdata and its range. normalize
Normalize the time series if normalize = 0 , keep the time series at their originalscale; if normalize = 1 , normalize the time series tounit interval; if normalize = 2 , z-score the time series. mindiagline
The minimum number of contiguouspoints along the diagonals to consider thesystem into a recurrent state The default value is usually 2, as it takes a minimumof two points to define any line. minvertline
The minimum number of contiguouspoints along the vertical lines to considerthe system into a recurrent state The default value is usually 2, as it takes a minimumof two points to define any line. tw The Theiler window parameter It defines the number of diagonals off the line of identitythat are excluded from recurrence quantification. whiteline
A logical flag to calculate (TRUE) or not(FALSE) empty vertical lines. The default is FALSE as the calculation of such linesadds on the time of computation. recpt
A logical flag indicating whether mea-sures of cross-recurrence are calculateddirectly from a recurrent plot (TRUE) ornot (FALSE) It is mostly useful if the user wants to compute joint-recurrence analysis. The RP or CRP is supplied inplace of ts1, and ts2 needs to be assigned as NA. side
A string indicating the side of the recur-rence plot on which recurrence measuresshould be calculated For side = upper , recurrence measures are calculatedon the upper triangle of the RP, for side = lower on the lower triangle of the RP, for size = both onthe full RP. Note, the line of identity is automaticallyexcluded for upper and lower setting. method
A string to indicate the type of recurrenceanalysis to perform For method = rqa , Auto-recurrence is calculated, i.e.,a unidimensional series; for method = crqa cross-recurrence is calculated; for method = mdcrqa multidi-mensional recurrence is calculated. Note, the defaultvalue is crqa . metric A string to indicate the type of distancemetric used To see the list of all other possible metrics see thehelp for the rdist() function. Note, the default is euclidean . datatype A string (continuous or categorical) toindicate the nature of the data type If the time series contain categorical information, itwill automatically be recoded into a continuous integer-based time series with a warning sent to the user.
Table 2: Overview of the arguments that can be used in crqa to set up the recurrencequantification analysis.to set up the recurrencequantification analysis.