Front. Appl. Math. Stat.Frontiers in Applied Mathematics and StatisticsFront. Appl. Math. Stat.2297-4687Frontiers Media S.A.10.3389/fams.2020.00009Applied Mathematics and StatisticsMethodsStudying Complex Adaptive Systems With Internal States: A Recurrence Network Approach to the Analysis of Multivariate Time-Series Data Representing Self-Reports of Human ExperienceHasselmanFred*BosmanAnna M. T.School of Pedagogical and Educational Sciences, Radboud University, Nijmegen, Netherlands
Edited by: Norbert Marwan, Potsdam-Institut für Klimafolgenforschung (PIK), Germany
Reviewed by: Miguel Pineda, University College London, United Kingdom; Zhong-Ke Gao, Tianjin University, China
*Correspondence: Fred Hasselman f.hasselman@pwo.ru.nl
This article was submitted to Dynamical Systems, a section of the journal Frontiers in Applied Mathematics and Statistics
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
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.
weighted recurrence networksmultiplex recurrence networkscomplex systems approachecological momentary assessmentexperience sampling methodidiographic analysisfractal scalingTime 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–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–6].
Our main concerns are about the validity of the following assumptions:
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.
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–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 non-informative 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 self-reports 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 self-reports 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.
Bi-daily ratings of self-esteem and perception of the physical self by four participants over the course of 512 days. The data were originally reported by Delignieres et al. [27].
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 long-range 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:
GPj=∑i=1jxi-x¯, for j=1,2,…,N
In Equation 1, GPj is the global profile of xi, a time series of i = 1, …, N equidistant observations, x¯ represents the time-series mean 1N∑i=1Nxi. 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:
w¯j=1W∑i=0W-1xi+j, with j =1,2,…,N-W
CPj in Equation 3 is the cumulative sum of the differences between xi and w¯j, but because w¯j is a shorter sequence than xi, one can choose how to align the two series. Using a left-aligned window refers to the differences taken starting at xi = 1, in a centered window at xi = W/2, and in a right-aligned window at xi = W. Equation 3 represents a right-aligned window:
CPj=∑i=0j-1xi+W-w¯j, with j =1,2,…, N-W
Figure 2 displays the CPs for the observed variable Perceived Fitness for each participant, in which the window size was W = 14 (corresponding to 7 days). The effect of calculating the CP is that any changes in the state evaluation relative to the rolling average are “exaggerated”; the description of the state dynamics is no longer bounded by the range of the ordinal rating scale (see the Supplementary Material for a more detailed description of the effect of window size on the CP). The observed series of Participant 1, for example, oscillate around a rather stationary level in the first half of the observation period, but displays a downward trend in the second half. This can also be seen in the series of Participants 3 and 4. The CPs indicate that the dynamics driving these downward trends are likely more diverse than what might be expected based on the fluctuations of observed values on the bounded ordinal scale. 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.
Change Profiles (blue lines) of the variable Perceived Fitness for each participant. The profiles were calculated using a right aligned window of 14 data points, representing 7 days (see Equation 1).
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 (k-spectrum [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 ki and vertex strength si (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 delay-embedded 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 Ri,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 Ri,j, denoted as ε.05 in Equation 4:
Ri,j(ε.05)=Θ(ε.05-‖x⃗i-x⃗j‖),i,j=1,…,N
(A) shows a distance matrix for the Change Profile of Perceived Strength of Participant 1. The series was embedded with parameters, m = 3 and τ = 37. (B) is a recurrence matrix, a threshold value ε = 2.08 was chosen to yield RR = 0.05. A weighted recurrence matrix is shown in (C), created by eliminating all values above the threshold of ε = 2.08 and replacing the distance values with recurrence times.
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 Ri,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):
Ri,j(ε.05)=Θ(ε.05-‖x⃗i-x⃗j‖).|i-j|,i,j=1,…,N
If a time stamp is available, the weights |i − j| can be expressed as a duration in the observed units of time. For the purpose of 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) shows one set of 51 (out of 100) simulated ideal noise series with log–log PSD slopes varying between −2.5 and 2.5. (B) represents the relation between the log–log PSD slopes and βRR~RT. Each of the 5,100 circles represents a recurrence network constructed from a simulated noise series. Circles are color coded according to the mean vertex strength s(k), and their size represents the mean vertex degree k. The Spearman correlation between log–log PSD slopes in the range [−2,0] and βRR~RT is 0.93.
Mean and standard deviation (SD) of the embedding parameters and recurrence network measures for the simulated series in Figure 4.
PSD slope
τ
SD
m
SD
βRR-RT
SD
s(k)
SD
k
SD
(−2.5, −1.5)
26.59
6.03
3.73
0.47
0.18
0.02
5.37
0.73
47.59
0.34
(−1.5, −0.5)
8.85
3.97
4.90
0.32
0.71
0.32
1.89
0.94
49.43
0.66
(−0.5, 0.5)
2.15
0.47
5.00
0.00
1.20
0.02
0.70
0.06
50.67
0.10
(0.5, 1.5)
2.20
0.42
4.75
0.42
1.21
0.02
0.64
0.02
50.69
0.02
(1.5, 2.5)
1.85
1.11
4.25
0.35
1.19
0.01
0.62
0.02
50.80
0.20
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 (rS = 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 rS = 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 (rS = −0.94), whereas the highest average vertex degrees are associated with positive log–log PSD slope (rS = 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.
The relation between the log–log PSD slopes and βRR~RT for each participant and each observed variable.
Relation between βRR−RT, predicted PSD slope, and estimated PSD slope for each participant and variable.
ID
Variable
βRR-RT
Predicted PSD slope
Estimated PSD slope
Difference
Participant 1
Attractive body
0.27
−1.56
−2.09
0.53
Participant 1
Global self-esteem
0.26
−1.57
−1.80
0.22
Participant 1
Perceived fitness
0.16
−1.71
−1.30
−0.41
Participant 1
Perceived strength
0.24
−1.61
−1.18
−0.42
Participant 1
Physical self-worth
0.23
−1.62
−1.83
0.21
Participant 1
Sport competence
0.19
−1.67
−0.91
−0.76
Participant 2
Attractive body
0.29
−1.53
−1.51
−0.02
Participant 2
Global self-esteem
0.30
−1.52
−0.85
−0.67
Participant 2
Perceived fitness
0.20
−1.66
−1.68
0.02
Participant 2
Perceived strength
0.25
−1.59
−2.22
0.63
Participant 2
Physical self-worth
0.25
−1.59
−1.49
−0.10
Participant 2
Sport competence
0.23
−1.62
−1.07
−0.55
Participant 3
Attractive body
0.26
−1.58
−1.58
−0.01
Participant 3
Global self-esteem
0.28
−1.56
−2.40
0.84
Participant 3
Perceived fitness
0.28
−1.54
−1.15
−0.39
Participant 3
Perceived strength
0.31
−1.51
−0.98
−0.53
Participant 3
Physical self-worth
0.30
−1.53
−2.12
0.59
Participant 3
Sport competence
0.30
−1.53
−2.03
0.50
Participant 4
Attractive body
0.25
−1.59
−1.43
−0.16
Participant 4
Global self-esteem
0.26
−1.58
−0.83
−0.75
Participant 4
Perceived fitness
0.34
−1.46
−1.89
0.43
Participant 4
Perceived strength
0.36
−1.44
−2.09
0.65
Participant 4
Physical self-worth
0.34
−1.46
−1.22
−0.24
Participant 4
Sport competence
0.24
−1.61
−1.27
−0.33
Mean
0.27
−1.57
−1.54
−0.03
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. Figure 6 displays a basic Archimedean spiral layout in which the vertices are equally divided over six arcs. The edge weights are set to the recurrence time frequency. The RN represents the CP of the variable Perceived Strength produced by Participant 1 and weighted by recurrence time.
A spiral layout of the recurrence network of the change profile of variable Perceived Strength produced by Participant 1. Edges with recurrence times of 10 or smaller (5 days) were removed to improve visibility.
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.
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.
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 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.
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.
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]:
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 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(k)[α]),P(s(k)[β]), by simply exchanging vertex degree for vertex strength:
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.
Multiplex weighted recurrence networks for each multivariate dataset. The layers of the network represent the RNs of the observed change profiles (AB, Attractive Body; PS, Physical Strength; GS, Global Self-Esteem; SW, Physical Self-Worth; SC, Sport Competence; PF, Perceived Fitness). The edge weights represent the multilayer mutual information, Iα,β, based on the strength distribution in each layer.
Layer strength, average edge overlap, and multilayer mutual information of the MwRNs in Figure 9.
Layer strength
Average
Average
ID
GS
PS
AB
PF
SC
SW
Edge overlap
Iα,β
Participant 1
10.44
10.44
9.67
9.67
10.44
10.44
0.024
2.04
Participant 2
9.60
10.03
8.44
10.13
10.36
10.15
0.028
1.96
Participant 3
11.61
11.61
11.36
11.61
11.36
11.36
0.035
2.30
Participant 4
10.99
10.64
10.99
9.91
10.60
10.04
0.017
2.11
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:
The data generating process is essentially stationary, homogeneous, and void of long-range dependencies.
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 memoryless-ness 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.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
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
ReferencesShiffmanS. Conceptualizing analyses of ecological momentary assessment data. Nicotine Tob Res. (2014) 16 (Suppl. 2):S76–87. 10.1093/ntr/ntt19524323571NezlekJ. Ecological momentary assessment: statistical issues. Eur Psychiatry. (2015) 30:121. 10.1016/S0924-9338(15)31837-XConnerTSTennenHFleesonWBarrettLF. Experience sampling methods: a modern idiographic approach to personality research. Soc Personal Psychol Compass. (2009) 3:292–313. 10.1111/j.1751-9004.2009.00170.x19898679BringmannLFElmerTEpskampSKrauseRWSchochDWichersM. What do centrality measures measure in psychological networks?J Abnorm Psychol. (2019) 128:892–903. 10.1037/abn000044631318245MolenaarPCM. Consequences of the ergodic theorems for classical test theory, factor analysis, and the analysis of developmental processes. In: HoferSMAlwinDF, editors. Handbook of Cognitive Aging: Interdisciplinary Perspectives. Thousand Oaks, CA: SAGE Publications, Inc. (2008). p. 90–104.van BorkRRhemtullaMWaldorpLJKruisJRezvanifarSBorsboomD. Latent variable models and networks: statistical equivalence and testability. Multivariate Behav Res. (2019) 16:1–24. 10.1080/00273171.2019.1672515MolenaarPC. On the implications of the classical ergodic theorems: analysis of developmental processes has to focus on intra-individual variation. Dev Psychobiol. (2008) 50:60–9. 10.1002/dev.2026218085558LuceRDNarensL. Symmetry, scale, types, and generalizations of classical physical measurment. J Math Psychol. (1983) 27:44–85. 10.1016/0022-2496(83)90026-3EpskampSCramerAOJWaldorpLJSchmittmannVDBorsboomD. qgraph: Network visualizations of relationships in psychometric data. J Stat Softw. (2012) 48:1–18. 10.18637/jss.v048.i04BeranJ. Statistical methods for data with long-range dependence. Stat Sci. (1992) 7:404–16. 10.1214/ss/1177011122PressWH. Reproducibility Now at Risk?Symposium on Evidence in the Natural Sciences 30-05-2014.New York, NY: Simons Foundation (2014).WallotSKelty-StephenDG. Interaction-dominant causation in mind and brain, and its implication for questions of generalization and replication. Minds Machines. (2017) 28:353–74. 10.1007/s11023-017-9455-0FisherAJMedagliaJDJeronimusBF. Lack of group-to-individual generalizability is a threat to human subjects research. Proc Natl Acad Sci USA. (2018) 115:E6106–E15. 10.1073/pnas.171197811529915059WolfersTDoanNTKaufmannTAlnaesDMobergetTAgartzI. Mapping the heterogeneous phenotype of schizophrenia and bipolar disorder using normative models. JAMA Psychiatry. (2018) 75:1146–55. 10.1001/jamapsychiatry.2018.246730304337CoxRFAHasselmanFSeevinckMP editors. A process approach to psychological measurement: lessons from quantum mechanics. In: 16th International Conference on Perception–Action. Ouro Preto (2011).WallaceD. Decoherence and its role in the modern measurement problem. Philos Trans A Math Phys Eng Sci. (2012) 370:4576–93. 10.1098/rsta.2011.049022908343TverskyA. A general theory of polynomial conjoint measurement. J Math Psychol. (1967) 4:1–20. 10.1016/0022-2496(67)90039-9BruzaPBusemeyerJRGaboraL. Introduction to the special issue on quantum cognition. J Math Psychol. (2009) 53:303–5. 10.1016/j.jmp.2009.06.002Van OrdenGCKelloCTHoldenJG. Situated behavior and the place of measurement in psychological theory. Ecol Psychol. (2010) 22:24–43. 10.1080/10407410903493145Van OrdenGCHoldenJGTurveyMT. Human cognition and 1/f scaling. J Exp Psychol Gen. (2005) 134:117–23. 10.1037/0096-3445.134.1.11715702967RamachandranB. On the “Strong Memorylessness Property” of the exponential and geometric probability laws. Sankhya: Indian J Stat. (1979) 41:244–51.WijnantsMLHasselmanFCoxRFBosmanAMVan OrdenG. An interaction-dominant perspective on reading fluency and dyslexia. Ann Dyslexia. (2012) 62:100–19. 10.1007/s11881-012-0067-322460607MetzlerRJeonJHCherstvyAGBarkaiE. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys Chem Chem Phys. (2014) 16:24128–64. 10.1039/C4CP03465A25297814ZouYDonnerRVMarwanNDongesJFKurthsJ. Complex network approaches to nonlinear time series analysis. Phys Rep. (2019) 787:1–97. 10.1016/j.physrep.2018.10.005MarwanNRomanoMCThielMKurthsJ. Recurrence plots for the analysis of complex systems. Phys Rep. (2007) 438:237–329. 10.1016/j.physrep.2006.11.001GaoZKLiuCYYangYXCaiQDangWDDuXL. Multivariate weighted recurrence network analysis of EEG signals from ERP-based smart home system. Chaos. (2018) 28:085713. 10.1063/1.501882430180616DelignieresDFortesMNinotG. The fractal dynamics of self-esteem and physical self. Nonlinear Dynamics Psychol Life Sci. (2004) 8:479–510. Available online at: http://www.societyforchaostheory.org/ndpls/askFILE.cgi?vol=08&iss=04&art=02&desc=ABSTRACT15473949BringmannLFHamakerELVigoDEAubertABorsboomDTuerlinckxF. Changing dynamics: time-varying autoregressive models using generalized additive modeling. Psychol Methods. (2017) 22:409–25. 10.1037/met000008527668421BianciardiMSirabellaPHagbergGEGiulianiAZbilutJPColosimoA. Model-free analysis of brain fMRI data by recurrence quantification. Neuroimage. (2007) 37:489–503. 10.1016/j.neuroimage.2007.05.02517600730ZbilutJPGiulianiAWebberCL. Detecting deterministic signals in exceptionally noisy environments using cross-recurrence quantification. Phys Lett A. (1998) 246:122–8. 10.1016/S0375-9601(98)00457-5DonnerRVSmallMDongesJFMarwanNZouYXiangRX. Recurrence-based time series analysis by means of complex network methods. Int J Bifurcat Chaos. (2011) 21:1019–46. 10.1142/S0218127411029021JacobRHarikrishnanKPMisraRAmbikaG. Characterization of chaotic attractors under noise: a recurrence network perspective. Commun Nonlinear Sci Numer Simul. (2016) 41:32–47. 10.1016/j.cnsns.2016.04.028BarratABarthelemyMPastor-SatorrasRVespignaniA. The architecture of complex weighted networks. Proc Natl Acad Sci USA. (2004) 101:3747–52. 10.1073/pnas.040008710115007165PackardNHCrutchfieldJPFarmerJDShawRS. Geometry from a time-series. Phys Rev Lett. (1980) 45:712–6. 10.1103/PhysRevLett.45.712MarwanN. How to avoid potential pitfalls in recurrence plot based data analysis. Int J Bifurcat Chaos. (2011) 21:1003–17. 10.1142/S0218127411029008WebberCLJrZbilutJP. Dynamical assessment of physiological systems and states using recurrence plot strategies. J Appl Physiol. (1994) 76:965–73. 10.1152/jappl.1994.76.2.9658175612HasselmanF. casnet: A toolbox for studying Complex Adaptive Systems and NETworks. R package version 0.1.3 ed2017. 10.17605/osf.io/nauhrLittleMAMcSharryPERobertsSJCostelloDAMorozIM. Exploiting nonlinear recurrence and fractal scaling properties for voice disorder detection. Biomed Eng Online. (2007) 6:23. 10.1186/1475-925X-6-2317594480OlthofMHasselmanFStrunkGAasBSchiepekGLichtwarck-AschoffA. Destabilization in self-ratings of the psychotherapeutic process is associated with better treatment outcome in patients with mood disorders. Psychother Res. (2019) 1:1–12. 10.1080/10503307.2019.1633484OlthofMHasselmanFStrunkGvan RooijMAasBHelmichMA. Critical fluctuations as an early-warning signal for sudden gains and losses in patients receiving psychotherapy for mood disorders. Clin Psychol Sci. (2019) 8:25–35. 10.1177/2167702619865969LacasaLNicosiaVLatoraV. Network structure of multivariate time series. Sci Rep. (2015) 5:15508. 10.1038/srep1550826487040ErogluDMarwanNStebichMKurthsJ. Multiplex recurrence networks. Phys Rev E. (2018) 97:012312. 10.1103/PhysRevE.97.01231229448424
Funding. This article was written while FH was employed by the School of Pedagogical and Educational Sciences at the Radboud University.