Studying Complex Adaptive Systems With Internal States: A Recurrence Network Approach to the Analysis of Multivariate Time-Series Data Representing Self-Reports of Human Experience

We discuss formal, theoretical, and practical issues with the statistical analysis of multivariate time-series data that represent self-reports of human experience, often referred to as Ecological Momentary Assessment (EMA) data or Experience Sampling Method (ESM) data. We argue that such time series likely violate the assumptions required for valid statistical inference, such as the memoryless-ness property (due to the presence of long-range temporal correlations) and ergodicity (due to non-stationarity and non-homogeneity of central moments). Moreover, we consider the common practice of interpreting outcomes of self-reports as if they were outcomes of classical physical measurements as extremely problematic and suggest to consider them as records of the temporal evolution of observables of a complex adaptive system with internal state dynamics. We propose to address some of these issues by analyzing the change profile instead of the observed time series, using recurrence-based analyses, specifically (multiplex) recurrence networks. We analyze a publicly available dataset in which four participants rated six questions about their self-esteem and physical self, twice a day over a period of 512 days and introduce the concept of recurrence networks weighted by recurrence time. The edge weights represent either recurrence times or recurrence time frequencies and results show that the scaling relation between vertex degree and vertex strength (the weighted variant of vertex degree) is associated to the scaling relation between frequency and spectral power based on the “raw” time series. We present a new spiral layout for recurrence networks that might be more appropriate when the detection of critical periods, regime shifts, and tipping points requires insight into the temporal order in which those events occur. We conclude that a complex systems approach to analyzing multivariate time series of self-reports of human experience is preferred over and above fitting statistical models like the Gaussian Graphical Model or its derivatives.

We discuss formal, theoretical, and practical issues with the statistical analysis of multivariate time-series data that represent self-reports of human experience, often referred to as Ecological Momentary Assessment (EMA) data or Experience Sampling Method (ESM) data. We argue that such time series likely violate the assumptions required for valid statistical inference, such as the memoryless-ness property (due to the presence of long-range temporal correlations) and ergodicity (due to non-stationarity and non-homogeneity of central moments). Moreover, we consider the common practice of interpreting outcomes of self-reports as if they were outcomes of classical physical measurements as extremely problematic and suggest to consider them as records of the temporal evolution of observables of a complex adaptive system with internal state dynamics. We propose to address some of these issues by analyzing the change profile instead of the observed time series, using recurrence-based analyses, specifically (multiplex) recurrence networks. We analyze a publicly available dataset in which four participants rated six questions about their self-esteem and physical self, twice a day over a period of 512 days and introduce the concept of recurrence networks weighted by recurrence time. The edge weights represent either recurrence times or recurrence time frequencies and results show that the scaling relation between vertex degree and vertex strength (the weighted variant of vertex degree) is associated to the scaling relation between frequency and spectral power based on the "raw" time series. We present a new spiral layout for recurrence networks that might be more appropriate when the detection of critical periods, regime shifts, and tipping points requires insight into the temporal order in which those events occur. We conclude that a complex systems approach to analyzing multivariate time series of self-reports of human experience is preferred over and above fitting statistical models like the Gaussian Graphical Model or its derivatives.

TIME SERIES OF SELF-REPORTS OF HUMAN EXPERIENCE: A COMPLEX SYSTEMS APPROACH
Technological advances have made collecting densely sampled physiological, behavioral, and psychological data from individuals relatively easy and inexpensive (e.g., sensors and apps on "wearables" and personal electronic devices). The availability of such large, multivariate time-series data of human physiology and well-being is of great interest to the social, behavioral, and medical sciences; however, there is at present no consensus on how to analyze such datasets or how to interpret the results obtained by different analytic approaches [1][2][3].
The challenges faced by researchers in these disciplines entail more than expanding their analytical toolbox of population statistics to include multivariate time-series analyses. Measurements of the temporal evolution of mental or emotional states are known as Experience Sampling or Ecological Momentary Assessment (e.g., daily ratings of "Today I feel happy" on an ordered categorical scale from 1 = "completely disagree" to 5 = "completely agree"). Analysis strategies of such datasets expose some well-known concerns about the validity of assumptions pertaining to the theoretical object of measurement and the nature of the data generating processes underlying human behavior and cognition [4][5][6].
Our main concerns are about the validity of the following assumptions: 1. The assumption that the ergodic theorems apply [7] to the theoretical objects of measurement and data generating processes: Ensemble averages of variables observed in samples of sufficiently many individuals are expected to be arbitrarily similar to the time averages of variables evolving over a sufficiently long interval of time, from any single initial condition. 2. The assumption that the interpretation of outcomes of psychological measurement is, or, should be, equivalent to classical physical measurement [8]: It is considered unproblematic to interpret a measurement outcome as a property of the theoretical object of measurement confounded by sources of random, additive, measurement noise, or sampling error.
The assumptions related to ergodicity (i.e., stationarity and homogeneity of central moments) are obviously important for making valid statistical inferences and generalizations. A common practice in psychopathology studies in which the Gaussian Graphical Model (GGM) is fitted to multivariate self-report data to estimate so-called "symptom networks" [9] is the use of rather ad hoc methods to achieve data reduction. For example, subsets of the multivariate time-series data are combined by simply taking their average or by using techniques like principal component analysis and exploratory factor analysis to achieve dimension reduction by estimating the network on the extracted component or factor structure. Such aggregation methods are known to be problematic if the assumptions of independence, stationarity, and homogeneity are violated [7,10].
Even if the core assumptions for an ergodic data generating process formally apply, one cannot rely on parameter estimates to converge on a characteristic expected value within the time scale of observation or magnitude of fluctuation. This occurs for example when the process samples from a stable distribution with one or more undefined central moments like the Cauchy distribution. This has led some scholars to suggest that "the very notion of probability may not make sense" [11] when studying complex systems with internal state dynamics. Recent observations of discrepancies between inferred statistical properties at the ensemble level (inter-individual) and the individual level (intra-individual), have been suggested as a cause of the so-called reproducibility crisis in the social and life sciences [11][12][13]. A study that observed a lack of "group-to-individual generalizability" in the context of psychopathology described the phenomenon as a threat to human subjects research: "In clinical research, diagnostic tests may be systematically biased and our classification systems may be at least partially invalid. In terms of theory development, we may have a misleading impression about the nature of psychological variables and their interactions" [13]. A study of the neuroanatomical phenotypes of schizophrenia and bipolar disorder [14] concluded: "This study found that group-level differences disguised biological heterogeneity and interindividual differences among patients with the same diagnosis. This finding suggests that the idea of the average patient is a noninformative construct in psychiatry that falls apart when mapping abnormalities at the level of the individual patient." The second concern is about the lack of a clear notion of how to incorporate the measurement context and the act of measurement of psychological variables into the description of a phenomenon [15]. It seems uncontroversial to suggest that the very act of asking someone to project their current internal state of happiness onto an arbitrary ordinal scale will interfere with their "true" state of happiness (if such a thing even exists without the measurement context). However, this actually implies that psychological measurement cannot be regarded as classical physical measurement, because the outcomes of selfreports emerge as the result of an interaction between the theoretical object of measurement and the measurement context. This causes measurement problems that are at least conceptually related to contemporary discussions about the ontological status of measurement in physics [16]. Resolutions have been proposed for measurement of psychological variables, for example, various types of conjoint measurement [8,17] or measurement models based on contemporary particle physics [15,18,19]. When the temporal evolution of internal states is concerned, additional issues arise due to the fact that living systems are subject to aging (loss of identity over time) and appear to be able to coordinate their current behavior relative to some record of previously experienced events (i.e., memory). In more general terms, the behavior of a complex adaptive system will display after-effects of interactions with its internal or external environment that extend far beyond any time scale that might be understood as a simple stochastic process with autoregressive components [10,20]. Time series of observables of living systems will often lack the memoryless-ness property [21,22], suggesting that, at least, anomalous rather than normal diffusion processes should be considered as a model for the data generating process [23].

Addressing (Some) Concerns: Recurrence Networks of Change Profiles
In the present paper, we discuss the benefits and challenges of considering multivariate time-series data that represent selfreports of human experience as the observables of a complex dynamical system with internal state dynamics (internal degrees of freedom). More specifically, we will consider the framework of (multiplex) recurrence networks [24,25] as the preferred analytic toolbox for describing the dynamics of multivariate time series collected using the Experience Sampling Method. Recurrence networks have been used to analyze multivariate physiological data, like human EEG [26]. A publicly available dataset will be used in which four participants evaluated six questions about their self-esteem and attitudes toward their physical self, twice a day, during 512 consecutive days [27]. Figure 1 displays the time series for each participant and each variable; the ratings were evaluated on a visual analog scale, and participants moved a slider between the values 0 and 10. The analyses in the original study revealed the presence of long-range dependence, power-law scaling in virtually all cases close to 1/f noise [27]. Most time series are obviously non-stationary with respect to level and trend; some appear to reveal clear transitions between dynamical regimes, such as periods of low and high variability in responses.
There are several advantages of using recurrence-based methods for time-series analysis over statistical models that require stationarity and (conditional) independence like the GGM [9], but also over models that can handle non-stationary signals like Time Varying Auto Regressive models (TV-AR) [28]. These models generally limit estimating temporal dynamics to the linear domain, and to very short time scales (e.g., lag-1 vector autoregression [9]), which severely limits the range of potential data generating processes that can be considered to underlie the observed series.
In general, recurrence analytic approaches to time-series analysis are as follows: • Essentially model-free and make few assumptions about the data [29]. • Describe linear as well as non-linear dynamical phenomena [25]. • Describe (transitions between) dynamical regimes, even in exceptionally noisy environments [30]. • Quantify recurrent dynamics across all available time scales [25,31]. • Quantify structural similarities between different time series represented as multiplex recurrence networks [24], suspending the need for potentially problematic aggregation and/or dimension reduction of multivariate time-series data.
By using recurrence-based methods, most of the concerns related to the assumption of ergodicity are addressed to the extent that one is able to quantify dynamical behavior associated with local ergodicity breaking. However, the concerns about the interpretation of measurement outcomes of EMA data remain, which is why we propose to transform the data prior to analysis.
In the context of the current paper, we can expect the temporal evolution of self-reported ratings of internal states projected onto an arbitrary ordinal scale to be affected by at least three factors. First, current ratings are likely evaluated relative to some moving average of previous ratings (memory). People might recall their evaluations of a couple of days ago, perhaps a week, but it is less likely they will recall the previous month and use it as a reference level. Second, changes in the function that projects the internal state onto the ordinal scale (observed as sudden shifts in level, trend, or variance) are expected to be unpredictable and nonlinear. For example, a sudden re-evaluation of the relation between internal state and values of the rating scale might be triggered by the occurrence of an event that is extreme in the sense that the scale boundaries no longer adequately represent the current state relative to previously experienced states. The projection of "Today I feel sad" onto a scale of 1-5 might change dramatically relative to previous evaluations due to the sudden loss of a loved one. Finally, it is unlikely that using the same arbitrary scale to simultaneously measure two or more different internal states within the same individual will yield projection functions that are approximately equivalent for each distinct state or that such projections will be affected in an equivalent way by the state history and experience of extreme events. Based on these assertions, we assume that the way in which the ratings change from one moment to the next relative to some (non-stationary) level is more informative than the sequence of absolute values of the scale itself. For example, in Figure 1, Participants 1, 2, and 3 stay within a very limited range of the values allowed for by the scale, but analyses reveal that these bounded fluctuations are not random at all, but represent a complex pattern of longrange dependence.
We suggest to use a change profile (CP) whenever the raw data consist of time series of self-reports on bounded ordinal scales. This transformation of the data takes the cumulative sum of deviations of observed values from an estimate of the reference level used to evaluate the current state. One could, for example, use a simple moving average in a sliding window (a linear filter) to represent the reference level. The window size should be a reasonable estimate of the lag of time across which a person's current ratings are still affected by their previous ratings. Such a lag could be estimated using a partial autocorrelation function (first lag with non-significant autocorrelations), a lagged mutual information function (e.g., lag at the first or global minimum), or a lag that makes sense based on theoretical and empirical knowledge of the process in question. Here, we consider a CP based on a running average in a right-aligned window.
The CP is a simple extension of the global profile (GP), which is obtained by first subtracting the mean of the time series from each observed value and subsequently taking the cumulative sum: In Equation 1, GP j is the global profile of x i , a time series of i = 1, . . . , N equidistant observations, x represents the timeseries mean 1 The main difference between the GP and CP is that instead of subtracting the constant x, a running average series w j calculated in a window of size W is used: is the cumulative sum of the differences between x i and w j , but because w j is a shorter sequence than x i , one can choose how to align the two series. Using a leftaligned window refers to the differences taken starting at x i=1 , in a centered window at x i=W/2 , and in a right-aligned window at x i=W . Equation 3 represents a right-aligned window: In subsequent analyses that quantify recurrent states, values that are similar in the observed bounded series will be dissimilar in the unbounded CP, which we believe to be more realistic given the context of self-reports of emotional and psychological states. The remainder of this paper is structured as follows: In Section Weighted Recurrence Network Analysis: Recurrence Times as Edge Weights, we construct recurrence networks from the CPs of the time series in Figure 1 and discuss a method for extracting a scaling exponent from a weighted recurrence network that is similar to scaling exponents obtained from fractal analysis techniques. In Section Swing on a Spiral: A Sequential Layout for Recurrence Networks, we present novel layouts for recurrence networks designed to provide insight into the occurrence of regime shifts and network measures within certain epochs. In Section Multiplex Weighted Recurrence Networks, we construct multiplex recurrence networks and evaluate multilayer mutual information based on similarities between the strength distributions of the layers in the network.

WEIGHTED RECURRENCE NETWORK ANALYSIS: RECURRENCE TIMES AS EDGE WEIGHTS
The goal of the original study reporting on the multivariate time series of the four participants in Figure 1 was to quantify the temporal structure of fluctuations in self-esteem ratings in terms of scaling exponents [27]. The results showed that all subjects produced time series that resemble long-memory processes that resemble 1/f noise. In this section, we explore whether recurrence networks can be used to extract similar information about the data, by creating networks that represent information about recurrent states as well as information about the frequency with which they occur.
To the best of our knowledge, this approach is different from previous studies using, e.g., the degree distribution (kspectrum [32]), the transitivity, or clustering dimension to quantify dimension characteristics [24]. Weighted recurrence networks have been constructed based on multivariate EEG data [26], by creating a joint recurrence matrix (from the recurrence matrices of simultaneously recorded EEG channels). The joint recurrence rates (RRs) served as the edge weights. We propose using recurrence times as edge weights. This allows for a comparison of the scaling relation between the vertex degree k i and vertex strength s i (sum of the number of edges times edge weights) [33]. It has been shown that, absent any correlations between edge weights and vertex degree, the average strength of vertices of degree k scales with k as s(k) ∼ k β with β = 1. Any deviations from β = 1 can be interpreted to arise from systematic dependencies between the edge weights and vertex degree. In the case of a recurrence network weighted by recurrence times, a scaling relation β > 1 represents the case in which the mean recurrence time of an observed state increases exponentially with the frequency with which an observed state will recur, relative to a randomly connected network. We will refer to the slope β as the recurrence rate-recurrence time dependency (irrespective of whether frequency or duration is used), abbreviated as β RR∼RT .

Weighted Recurrence Network Construction
For each participant, we constructed distance matrices by delay embedding [34] the CPs of the six self-esteem variables using a window size of 14 (7 days). We chose a common delay τ and embedding dimension m for all six time series observed within a participant, using the median of the six τ values (found as first local minima of the lagged mutual information function) and the maximum of the six embedding dimensions (found by false nearest neighbor analysis [25,35,36]). Settling on one set of parameters for all variables observed within an individual is necessary for the construction of a multiplex network for each participant (see Section Multiplex Weighted Recurrence Networks), which requires the same number of vertices for each layer of the network [24]. Figure 3A (created with R package casnet [37]) displays the distance matrix of the delayembedded CP of Perceived Strength for Participant 1. To turn the distance matrix into a recurrence plot (RP, Figure 3B), which is a visualization of the recurrence matrix R i,j , a threshold value ε is chosen. The color bar in Figure 3A represents the range of distance values in the matrix. On the right side of the color bar, distance values are reported that would result in the RR reported on the left side of the color bar, should they be chosen as the threshold value. In the present paper, we chose ε such that RR would be 5% for all R i,j , denoted as ε .05 in Equation 4: The RP in Figure 3B was constructed with ε .05 = 2.08, with N as the length of the state space vectors (series length − (m − 1) · τ ), Θ as the Heaviside function (returning 0 if the distance exceeds ε .05 , 1 otherwise), and · as the Euclidean norm. To create a version of R i,j weighted by recurrence time, the result of the Heaviside function (0 or 1) is multiplied by the temporal separation between two state coordinates ( Figure 3C):    comparing the strength-degree scaling relation s(k) ∼ k β to the scaling relation between frequency and power obtained from a power spectrum of the time series, we represent the weights as a frequency 1 |i−j | .

Power-Law Scaling in Time Series and Weighted Networks
First, we study how the scaling relation s(k) ∼ k β (in the present context with β RR∼RT ) behaves when we create weighted recurrence networks of simulated ideal colored noise series. Time series were generated using a procedure based on the inverse Fourier transform that allows constructing a specific scaling relation between frequency and power in the frequency domain (a log-log slope in Power Spectral Density, see [38]). Figure 4A displays a set of 51 simulated time series of which the PSD slope varies from −2.5 to 2.5 in 0.1 increments; that is, noises include Brownian (red) noise (−2) to pink (−1), white (0), blue noise (1), and violet noise (2). For each of the 51 scaling relations, we generated 100 noise series (each with different random seeds) and created recurrence networks, choosing the delay τ and embedding dimension m as described in the previous paragraph (see Table 1).
A threshold ε was selected that yielded a 5% RR. This means the average vertex degree of each network will be approximately the same, but will vary with the length of the embedded series. Figure 4B shows the relation between β RR∼RT of each recurrence network and the log-log PSD slope from the simulated series from which it was constructed. The plot suggests a logistic relation, which has an almost perfect linear relationship (r S = 0.93) for noise series with log-log PSD slopes in the range of [−2, 0]. Outside of this range, the β RR∼RT slope does not seem to be able to distinguish very well between log-log PSD slopes; the Spearman correlation for the full range is r S = TABLE 1 | Mean and standard deviation (SD) of the embedding parameters and recurrence network measures for the simulated series in Figure 4. 0.75. Another interesting pattern that can be seen in Figure 4B is that the highest mean vertex strengths are associated with negative log-log PSD slopes (r S = −0.94), whereas the highest average vertex degrees are associated with positive log-log PSD slope (r S = 0.85). Compared to the observed series in Figure 1, the simulated series are relatively stationary, homogeneous, and void of noise. When the scaling relation β RR∼RT is calculated based on the observed data (Figure 5), we arrive at the same conclusion as the original study [27]: These observed time-series display a structure of temporal dependencies that resembles the category of colored noises between pink noise (PSD slope −1) and red noise (PSD slope −2). Table 2 displays the values of β RR∼RT and the estimated log-log PSD slope for each participant and variable. In addition, we use the linear relationship for log-log PSD slope ranges of [−2, 0] to predict the log-log PSD slope from β RR∼RT , which is approximately: PSD slope ≈ −1.94 + 1.38 · β RR∼RT . As can be seen in the final row of Table 2, the deviation between the predicted and estimated PSD slope is rather small on average.

SWING ON A SPIRAL: A SEQUENTIAL LAYOUT FOR RECURRENCE NETWORKS
Graphical representations of recurrence networks often use layouts that emphasize certain structural properties (e.g., "spring" layouts) and discard information about the temporal order, the vertex sequence. However, in a clinical setting, e.g., if multivariate time-series data are used to detect imminent regime shifts, such as sudden gains and losses in symptom severity [39,40], it is important to retain temporal information. Here, we present a simple layout for recurrence networks based on the spiral.  In Figure 6, each arc corresponds to an epoch, which is color coded, but epochs could also describe months, years, phases in an experiment, or other qualitative descriptors of segments of time. If edges connect vertices within the same epoch, they are represented in the vertices; otherwise, they are gray (see the Supplementary Material for spiral RNs of each participant and each variable in the dataset). This provides a rough indication of differences between longer and shorter recurrence times.
Although the estimated scaling exponents of the time series for each participant fall within a similar range [27], Figure 7 shows that the recurrence network on a spiral suggests that the power laws arise from rather different patterns of recurrent states, occurring at different moments in time. The RN of Participant 3 reveals lots of recurrent states in the first two arcs, and very few later on in the time series, whereas the RN of Participant 4 indicates most recurrent states occurring in the final arc and across different arcs (longer recurrence times). It is possible to change the spiral layout to highlight different structural aspects of the RN. Figure 8 shows a number of different layout options for the CP of Participant 3 in Figure 7. The Archimedean spiral (top left) in Figure 8 is similar to the layout in Figures 6, 7, except for the start and end of the time series, which are reversed. In Archimedean spirals, points on the spiral that intersect a line going through the origin will all be equidistant. This means the first part instead of the last part of the time series will be emphasized if the order is reversed, because it is now displayed on the longest arcs of the spiral. In the Bernoulli and Fermat spiral layout (Figure 8, top right and bottom left), the length of the arcs increases with each FIGURE 7 | Spiral graphs with four arcs of the variable Global Self-Esteem for each participant. Edges with recurrence times of 10 or smaller (5 days) were removed to improve visibility. turn, which emphasizes the part of the series represented on the outer arcs. The Euler spiral layout can emphasize the middle part of the time series.

MULTIPLEX WEIGHTED RECURRENCE NETWORKS
Multiplex recurrence networks provide a framework for studying the temporal structure in multivariate time series. To obtain a multiplex weighted recurrence network (MwRN), we simply consider the six recurrence time frequency weighted RNs for each participant as the layers of a larger network, the MwRN. The structural association between layers α and β in an unweighted MRN can be quantified by the interlayer mutual information measure I α,β that is based on the degree distributions P κ [α] , P κ [β] of the vertices in each layer [41,42]: (6) The term P κ [α] , κ [β] denotes the joint probability of observing the same vertex with degree κ [α] in layer α and degree κ [β] in layer β. P κ [α] and P κ [β] represent the marginal probabilities of observing the degree distribution in each layer.
Taking I α,β as the weight between two layers in a multiplex network can be interpreted as the quantification of structural FIGURE 8 | A showcase of spiral graphs using the variable Global Self-Esteem of Participant 3 in Figure 7. Each spiral type (Archimedean, Bernoulli, Fermat, and Euler) can highlight different aspects of the temporal structure of recurrent states. Edges with recurrence times of 10 or smaller (5 days) were removed to improve visibility. regularities between the time series in the multivariate dataset. This is of course a more abstract-level interpretation of the mutual information quantity: The averaged mutual information across all layers, I α,β , in a multiplex network can be interpreted as a quantification of typical information flow between the dimensions of the multivariate system [41,42].
For the weighted RNs, we construct a multilayer mutual information based on the strength distributions in each layer, e.g., P s(κ) [α] , P s(κ) [β] , by simply exchanging vertex degree for vertex strength: (7) Figure 9 displays the MwRNs for each participant; the edges connecting a layer α and β represent the mutual information I α,β based on the vertex strength distribution in those layers. I α,β can be said to quantify the typical information flow between the layers of the MwRN; in the present case, this can be interpreted as the characteristic way in which self-reports of internal states related to self-esteem and physical self-worth are associated across time scales spanning from half a day to 512 days. Different patterns of similar structure can be identified among participants. The MwRN of Participant 3 reveals that associations between Global Self-Esteem, Physical Strength, and Perceived Fitness are highest, whereas the strongest structural similarities in the MwRN of Participant 4 occur between the layers Global Self-Esteem, Sport Competence, and Attractive Body. Participants 1 and 2 show more diverse flow patterns.
In Table 3, several MwRN measures were reported: The Layer  Strength for each layer in the MwRN, which in the present case is just the number of layers times the value of I α,β (the vertex degree is always 6, because we did not impose a threshold to remove edges). The Average Edge Overlap measure represents the proportion of edges that are shared between any two vertices across all layers [42]. The average I α,β is simply the mean of the edge weights.

CONCLUSION AND DISCUSSION
A primary goal of the present paper was to address some of the issues that arise when multivariate time series of self-reports of human experience are analyzed and interpreted under the assumption that: 1. The data generating process is essentially stationary, homogeneous, and void of long-range dependencies. 2. The measurement process can be considered classical physical measurement.
We explored a publicly available dataset of which it is known that the multivariate time series represent non-stationary data with long-range dependence [27]. First, we showed that focusing on how internal states change relative to a non-stationary reference level, the CP seems to be a more sensible observable than the sequence of values allowed for by the arbitrary ordinal rating scale. A secondary goal was to show that analytic methods developed to study complex adaptive systems can be applied to time series of self-reports of human experience (EMA data) and that the results can be sensibly interpreted and visualized. We were able to construct recurrence networks using regular methods for delay-embedding time series (mutual information, false nearest neighbor analysis). Moreover, we study recurrence networks weighted by recurrence time (or recurrence time frequency) and show that a scaling relation β RR∼RT can be estimated that is very similar to the scaling relation estimated from the power spectrum, for a specific range of scaling exponents [−2,0].
Many of the potential applications of recurrence networks to EMA data would benefit from a visualization of the network in which the temporal order of observed values is retained and we suggest that the spiral layout might be beneficial in this respect. Finally, we constructed multilayer recurrence networks that revealed that characteristic patterns of fluctuations in self-reports about self-esteem vary across participants and reveal a limited subset of variables sharing similar structure over time in some participants, whereas other participants appear to share structure in temporal dynamics that include almost all observed variables. Of course, the present study is in fact "merely" a multiple case study, and more research will be required in order to decide whether the approach presented here generalizes to other datasets as well. We do believe that the present paper provides at least a proof of principle with respect to the applicability of the recurrence network approach to EMA data, which can potentially be a valuable analytic method for idiographic research, especially in psychopathology.
To summarize, we suggest that using a recurrence-based approach to analyze multivariate time-series data that can be considered to emerge from a complex adaptive system with internal state dynamics should be preferred over analytic methods that attempt to fit the parameters of a statistical model. The data generating processes underlying the multivariate series will almost certainly violate the key assumptions of stationary and homogeneity, and likely do not possess the memorylessness property.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study can be found on the website of the first authors of the original publication [27]: https:// didierdelignieresblog.wordpress.com/recherche/databank/.

AUTHOR CONTRIBUTIONS
FH conceived of hypotheses, methods, analyses, and structure. FH and AB wrote the article.

FUNDING
This article was written while FH was employed by the School of Pedagogical and Educational Sciences at the Radboud University.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams. 2020.00009/full#supplementary-material