# Recurrence Networks

A recurrence network is a representation of a recurrence matrix as a graph (a network). The nodes in the network represent time points, and if a value at any point in time will recur at some later point in time, an edge is drawn between the time points.

## Recurrence Networks

Recurrence networks are graphs created from a recurrence matrix. This means the nodes of the graph represent time points and the connections between nodes represent a recurrence relation betwen the values observed at those time points. That is, often the matrix represents recurrences in a reconstructed state space, the values are coordinates and therefore we would say the edges of a recurrence network represent a temporal relation between recurring states. The ultimate reference for learning everything about recurrence networks is:

Zou, Y., Donner, R. V., Marwan, N., Donges, J. F., & Kurths, J. (2018). Complex network approaches to nonlinear time series analysis. Physics Reports

Package casnet has some functions to create recurrence networks, they are similar to the functions used for CRQA: * rn() is very similar to rp(), it will create a matrix based on embedding parameters. One difference is the option to create a weighted matrix. This is a matrix in which non-recurring values are set to 0, but the recurring values are not replaced by a 1, the distance value is retained and acts as an edge-weight * rn_plot() will produce the same as rp_plot()

We can turn the recurrence matrix into an adjecency matrix, an igraph object. This means we can use all the igraph functions to calculate network measures.

library(igraph)
#library(qgraph)
#library(survival)
# Reload the data we used earlier
series <- rio::import("https://github.com/complexity-methods/CSA-assignments/raw/master/assignment_data/BasicTSA_arma/series.xlsx")

# Lets use a shorter dataset to speed things up
series <- series[1:500,]

We’ll analyse the three time series as a recurrence network: * Compare the results by comparing the different measures - Remember: TS_1 was white noise, TS_2 was a sine with added noise, TS_3 was the logistic map in the chaotic regime. * Note that some of the RQA measures can be exactly calculated from the measures of the network representation. - Try to understand why the Recurrence is represented as the degree centrality of the network (igraph::centr_degree())

### TS 1

#----------------------
# Adjacency matrix TS_1
#----------------------
plot(ts(series$TS_1))  arcs <- 6 # Because these are generated signals, look for a drop in FNN below 1%. p1 <- est_parameters(y = series$TS_1, nnThres = 1)


# By passing emRad = NA, a radius will be calculated
RN1 <- casnet::rn(y1 = series$TS_1, emDim = p1$optimDim, emLag = p1$optimLag, emRad = NA, targetValue = 0.05) > > Auto-recurrence: Setting diagonal to (1 + max. distance) for analyses > > Searching for a radius that will yield 0.05 for RR casnet::rn_plot(RN1)  # Get RQA measures rqa1 <- casnet::rp_measures(RN1, silent = FALSE) > > ~~~o~~o~~casnet~~o~~o~~~ > > Global Measures > Global Max.rec.points N.rec.points Recurrence.Rate Singular.points > 1 Recurrence Matrix 244530 12228 0.05000613 11140 > Divergence Repetitiveness Anisotropy > 1 0.1428571 1.58364 1 > > > Line-based Measures > Line.based N.lines N.points.on.lines Measure Rate Mean Max > 1 Diagonal 420 1088 Determinism 0.08897612 2.590476 7 > 2 Vertical 774 1723 V Laminarity 0.14090612 2.226098 6 > 3 Horizontal 774 1723 H Laminarity 0.14090612 2.226098 6 > Entropy.of.lengths Relative.entropy CoV.of.lengths > 1 1.0355136 0.16695004 0.3890404 > 2 0.5764252 0.09293379 0.2592375 > 3 0.5764252 0.09293379 0.2592375 > > ~~~o~~o~~casnet~~o~~o~~~ # Create RN graph g1 <- igraph::graph_from_adjacency_matrix(RN1, mode="undirected", diag = FALSE) igraph::V(g1)$size <- igraph::degree(g1)
g1r <- casnet::make_spiral_graph(g1,arcs = arcs, epochColours = getColours(arcs), markTimeBy = TRUE)


# Get RN measures
rn1 <- rn_measures(g1, silent = FALSE)
>
> ~~~o~~o~~casnet~~o~~o~~~
>
> Global Network Measures
>
>   EdgeDensity MeanStrengthDensity GlobalClustering NetworkTransitivity
> 1  0.05000613                  NA        0.5104217           0.5161686
>   AveragePathLength GlobalEfficiency
> 1          51.95484         3.212035
>
> ~~~o~~o~~casnet~~o~~o~~~

# Should be the same
rqa1$RR==rn1$graph_measures$EdgeDensity > [1] TRUE ### TS 2 #---------------------- # Adjacency matrix TS_2 #---------------------- plot(ts(series$TS_2))


# Because these are generated signals, look for a drop in FNN below 1%.
p2 <- est_parameters(y = series$TS_2, nnThres = 1)  RN2 <- rn(y1 = series$TS_2, emDim = p2$optimDim, emLag = p2$optimLag, emRad = NA,  targetValue = 0.05)
>
> Auto-recurrence: Setting diagonal to (1 + max. distance) for analyses
>
> Searching for a radius that will yield 0.05 for RR
rn_plot(RN2)


# Get RQA measures
rqa2 <- rp_measures(RN2, silent = FALSE)
>
> ~~~o~~o~~casnet~~o~~o~~~
>
>  Global Measures
>              Global Max.rec.points N.rec.points Recurrence.Rate Singular.points
> 1 Recurrence Matrix         244530        12228      0.05000613            9874
>   Divergence Repetitiveness Anisotropy
> 1  0.1111111        1.52294          1
>
>
>  Line-based Measures
>   Line.based N.lines N.points.on.lines      Measure      Rate     Mean Max
> 1   Diagonal     804              2354  Determinism 0.1925090 2.927861   9
> 2   Vertical    1552              3585 V Laminarity 0.2931796 2.309923   8
> 3 Horizontal    1552              3585 H Laminarity 0.2931796 2.309923   8
>   Entropy.of.lengths Relative.entropy CoV.of.lengths
> 1          1.3217260        0.2130945      0.4538389
> 2          0.7052368        0.1137014      0.3101225
> 3          0.7052368        0.1137014      0.3101225
>
> ~~~o~~o~~casnet~~o~~o~~~

# Create RN graph
g2 <- igraph::graph_from_adjacency_matrix(RN2, mode="undirected", diag = FALSE)
V(g2)$size <- degree(g2) g2r <- make_spiral_graph(g2,arcs = arcs ,epochColours = getColours(arcs), markTimeBy = TRUE)  # Get RN measures rn2 <- rn_measures(g2, silent = FALSE) > > ~~~o~~o~~casnet~~o~~o~~~ > > Global Network Measures > > EdgeDensity MeanStrengthDensity GlobalClustering NetworkTransitivity > 1 0.05000613 NA 0.5781826 0.5369779 > AveragePathLength GlobalEfficiency > 1 7.020586 3.510234 > > ~~~o~~o~~casnet~~o~~o~~~ # Should be the same rqa2$RR==rn2$graph_measures$EdgeDensity
> [1] TRUE

### TS 3

#----------------------
# Adjacency matrix TS_3
#----------------------
plot(ts(series$TS_3))  # Because these are generated signals, look for a drop in FNN below 1%. p3 <- est_parameters(y = series$TS_3, nnThres = 1)


RN3 <- rn(y1 = series$TS_3, emDim = p3$optimDim, emLag = p3$optimLag, emRad = NA, targetValue = 0.05) > > Auto-recurrence: Setting diagonal to (1 + max. distance) for analyses > > Searching for a radius that will yield 0.05 for RR rn_plot(RN3)  # Get RQA measures rqa3 <- rp_measures(RN3, silent = FALSE) > > ~~~o~~o~~casnet~~o~~o~~~ > > Global Measures > Global Max.rec.points N.rec.points Recurrence.Rate Singular.points > 1 Recurrence Matrix 232806 11642 0.0500073 7662 > Divergence Repetitiveness Anisotropy > 1 0.1666667 0.170603 1 > > > Line-based Measures > Line.based N.lines N.points.on.lines Measure Rate Mean Max > 1 Diagonal 1776 3980 Determinism 0.34186566 2.240991 6 > 2 Vertical 324 679 V Laminarity 0.05832331 2.095679 3 > 3 Horizontal 324 679 H Laminarity 0.05832331 2.095679 3 > Entropy.of.lengths Relative.entropy CoV.of.lengths > 1 0.6092768 0.09862128 0.2508340 > 2 0.3154837 0.05106613 0.1405776 > 3 0.3154837 0.05106613 0.1405776 > > ~~~o~~o~~casnet~~o~~o~~~ # Create RN graph g3 <- igraph::graph_from_adjacency_matrix(RN3, mode="undirected", diag = FALSE) V(g3)$size <- degree(g3)
g3r <- make_spiral_graph(g3,arcs = arcs ,epochColours = getColours(arcs), markTimeBy = TRUE)


# Get RN measures
rn3 <- rn_measures(g3, silent = FALSE)
>
> ~~~o~~o~~casnet~~o~~o~~~
>
> Global Network Measures
>
>   EdgeDensity MeanStrengthDensity GlobalClustering NetworkTransitivity
> 1   0.0500073                  NA        0.6030504            0.597613
>   AveragePathLength GlobalEfficiency
> 1          4.010532         3.260776
>
> ~~~o~~o~~casnet~~o~~o~~~

# Should be the same
rqa3$RR==rn3$graph_measures$EdgeDensity > [1] TRUE ## Multiplex Recurrence Networks Consider the three time series to be part of a multi-layer recurrence network. Common properties of the multiplex network are inter-layer mutual information and edge overlap can be calculated using function casnet::mrn(). One problem, the networks have to be all of the same size (same number of nodes, a multivariate time series), but here we have reconstructed the phase space using different embedding parameters… let’s choose one set of parameters for all time series. emDim <- mean(c(p1$optimDim,p2$optimDim,p3$optimDim))
emLag <- median(c(p1$optimLag,p2$optimLag,p3$optimLag)) RNs <- plyr::llply(1:3, function(r) rn(y1 = series[,r], emDim = emDim, emLag = emLag, emRad = NA, targetValue = 0.05)) > > Auto-recurrence: Setting diagonal to (1 + max. distance) for analyses > > Searching for a radius that will yield 0.05 for RR > > Auto-recurrence: Setting diagonal to (1 + max. distance) for analyses > > Searching for a radius that will yield 0.05 for RR > > Auto-recurrence: Setting diagonal to (1 + max. distance) for analyses > > Searching for a radius that will yield 0.05 for RR layers <- plyr::llply(RNs, function(r) igraph::graph_from_adjacency_matrix(r, mode="undirected", diag = FALSE)) names(layers) <- c("g1","g2","g3") mrn(layers = layers) > > ~~~o~~o~~casnet~~o~~o~~~ > > Welcome to the multiplex... in layer similarity mode! > > > ~~~o~~o~~casnet~~o~~o~~~ >$interlayerMI
> $interlayerMI$window: 1 | start: 1 | stop: 495
>    g1        g2        g3
> g1 NA 0.2609506 0.3133915
> g2 NA        NA 1.4864990
> g3 NA        NA        NA
>
>
> $edgeOverlap >$edgeOverlap$window: 1 | start: 1 | stop: 495 > g1 g2 g3 > g1 NA 0.02568631 0.06124579 > g2 NA NA 0.02625123 > g3 NA NA NA > > >$meanValues
> $meanValues$window: 1 | start: 1 | stop: 495
>     mi_mean     mi_sd    eo_mean      eo_sd eo_sum eo_all    eo_joint
> 1 0.6869471 0.6929286 0.03772777 0.02036916  18572     34 0.001830713
>
>
> $Nsize > [1] 495 A variety of plots can be created using casnet::mrn_plot()  # Simple mrn_plot(layers = layers, showEdgeColourLegend =TRUE) > > ~~~o~~o~~casnet~~o~~o~~~ > > Welcome to the multiplex... in layer similarity mode! > > > ~~~o~~o~~casnet~~o~~o~~~ >$MRN
> $MRN$window: 1 | start: 1 | stop: 495

>
>
> $interlayerMI >$interlayerMI$window: 1 | start: 1 | stop: 495 > g1 g2 g3 > g1 NA 0.2609506 0.3133915 > g2 NA NA 1.4864990 > g3 NA NA NA > > >$edgeOverlap
> $edgeOverlap$window: 1 | start: 1 | stop: 495
>    g1         g2         g3
> g1 NA 0.02568631 0.06124579
> g2 NA         NA 0.02625123
> g3 NA         NA         NA
>
>
> $meanValues > NULL mrn_plot(layers = layers, MRNweightedBy = "EdgeOverlap", showEdgeColourLegend =TRUE) > > ~~~o~~o~~casnet~~o~~o~~~ > > Welcome to the multiplex... in layer similarity mode! > > > ~~~o~~o~~casnet~~o~~o~~~ >$MRN
> $MRN$window: 1 | start: 1 | stop: 495

>
>
> $interlayerMI >$interlayerMI$window: 1 | start: 1 | stop: 495 > g1 g2 g3 > g1 NA 0.2609506 0.3133915 > g2 NA NA 1.4864990 > g3 NA NA NA > > >$edgeOverlap
> $edgeOverlap$window: 1 | start: 1 | stop: 495
>    g1         g2         g3
> g1 NA 0.02568631 0.06124579
> g2 NA         NA 0.02625123
> g3 NA         NA         NA
>
>
> $meanValues > NULL # Include picture of Layers mrn_plot(layers = layers, RNnodes = TRUE) > > ~~~o~~o~~casnet~~o~~o~~~ > > Welcome to the multiplex... in layer similarity mode! > > > ~~~o~~o~~casnet~~o~~o~~~ >$MRN
> $MRN$window: 1 | start: 1 | stop: 495

>
>
> $interlayerMI >$interlayerMI$window: 1 | start: 1 | stop: 495 > g1 g2 g3 > g1 NA 0.2609506 0.3133915 > g2 NA NA 1.4864990 > g3 NA NA NA > > >$edgeOverlap
> $edgeOverlap$window: 1 | start: 1 | stop: 495
>    g1         g2         g3
> g1 NA 0.02568631 0.06124579
> g2 NA         NA 0.02625123
> g3 NA         NA         NA
>
>
> $meanValues > NULL mrn_plot(layers = layers, RNnodes = TRUE,MRNweightedBy = "EdgeOverlap", showEdgeColourLegend =TRUE) > > ~~~o~~o~~casnet~~o~~o~~~ > > Welcome to the multiplex... in layer similarity mode! > > > ~~~o~~o~~casnet~~o~~o~~~ >$MRN
> $MRN$window: 1 | start: 1 | stop: 495

>
>
> $interlayerMI >$interlayerMI$window: 1 | start: 1 | stop: 495 > g1 g2 g3 > g1 NA 0.2609506 0.3133915 > g2 NA NA 1.4864990 > g3 NA NA NA > > >$edgeOverlap
> $edgeOverlap$window: 1 | start: 1 | stop: 495
>    g1         g2         g3
> g1 NA 0.02568631 0.06124579
> g2 NA         NA 0.02625123
> g3 NA         NA         NA
>
>
> $meanValues > NULL ## Time-varying Multiplex Recurrence Networks The MRN can be calculated in a sliding window, by setting the arguments win and step.  # This will generate 26 windows MRN_win <- mrn(layers = layers, win = 250, step = 10) > > ~~~o~~o~~casnet~~o~~o~~~ > > Welcome to the multiplex... in layer similarity mode! > > > ~~~o~~o~~casnet~~o~~o~~~ # The MRN are returned as a list MRN_win$interlayerMI
> $window: 01 | start: 1 | stop: 250 > g1 g2 g3 > g1 NA 0.3377082 0.5216094 > g2 NA NA 1.6939426 > g3 NA NA NA > >$window: 02 | start: 11 | stop: 260
>    g1        g2        g3
> g1 NA 0.2687605 0.5649311
> g2 NA        NA 1.6881942
> g3 NA        NA        NA
>
> $window: 03 | start: 21 | stop: 270 > g1 g2 g3 > g1 NA 0.3708465 0.4491544 > g2 NA NA 1.6464434 > g3 NA NA NA > >$window: 04 | start: 31 | stop: 280
>    g1       g2        g3
> g1 NA 0.368284 0.4455847
> g2 NA       NA 1.7394935
> g3 NA       NA        NA
>
> $window: 05 | start: 41 | stop: 290 > g1 g2 g3 > g1 NA 0.3445478 0.5649311 > g2 NA NA 1.8526006 > g3 NA NA NA > >$window: 06 | start: 51 | stop: 300
>    g1        g2        g3
> g1 NA 0.2969798 0.5897381
> g2 NA        NA 1.8290879
> g3 NA        NA        NA
>
> $window: 07 | start: 61 | stop: 310 > g1 g2 g3 > g1 NA 0.475031 0.3404266 > g2 NA NA 1.6741508 > g3 NA NA NA > >$window: 08 | start: 71 | stop: 320
>    g1        g2        g3
> g1 NA 0.4140999 0.5532296
> g2 NA        NA 1.8609670
> g3 NA        NA        NA
>
> $window: 09 | start: 81 | stop: 330 > g1 g2 g3 > g1 NA 0.5824751 0.4637912 > g2 NA NA 1.9291099 > g3 NA NA NA > >$window: 10 | start: 91 | stop: 340
>    g1        g2        g3
> g1 NA 0.3726358 0.4256365
> g2 NA        NA 1.6247560
> g3 NA        NA        NA
>
> $window: 11 | start: 101 | stop: 350 > g1 g2 g3 > g1 NA 0.3427359 0.5163049 > g2 NA NA 2.0026265 > g3 NA NA NA > >$window: 12 | start: 111 | stop: 360
>    g1        g2        g3
> g1 NA 0.5663371 0.5236299
> g2 NA        NA 1.7002084
> g3 NA        NA        NA
>
> $window: 13 | start: 121 | stop: 370 > g1 g2 g3 > g1 NA 0.4577888 0.5714332 > g2 NA NA 1.7548832 > g3 NA NA NA > >$window: 14 | start: 131 | stop: 380
>    g1       g2        g3
> g1 NA 0.424447 0.5894762
> g2 NA       NA 1.7602788
> g3 NA       NA        NA
>
> $window: 15 | start: 141 | stop: 390 > g1 g2 g3 > g1 NA 0.673275 0.5958021 > g2 NA NA 1.9566326 > g3 NA NA NA > >$window: 16 | start: 151 | stop: 400
>    g1        g2        g3
> g1 NA 0.5742539 0.6372795
> g2 NA        NA 1.9088293
> g3 NA        NA        NA
>
> $window: 17 | start: 161 | stop: 410 > g1 g2 g3 > g1 NA 0.6137142 0.559716 > g2 NA NA 1.882073 > g3 NA NA NA > >$window: 18 | start: 171 | stop: 420
>    g1       g2        g3
> g1 NA 0.466476 0.5714332
> g2 NA       NA 2.0435263
> g3 NA       NA        NA
>
> $window: 19 | start: 181 | stop: 430 > g1 g2 g3 > g1 NA 0.506056 0.5119128 > g2 NA NA 2.0341960 > g3 NA NA NA > >$window: 20 | start: 191 | stop: 440
>    g1        g2        g3
> g1 NA 0.4166984 0.5894762
> g2 NA        NA 1.8045100
> g3 NA        NA        NA
>
> $window: 21 | start: 201 | stop: 450 > g1 g2 g3 > g1 NA 0.407035 0.4257224 > g2 NA NA 1.8711604 > g3 NA NA NA > >$window: 22 | start: 211 | stop: 460
>    g1        g2        g3
> g1 NA 0.5555665 0.4790611
> g2 NA        NA 1.5796812
> g3 NA        NA        NA
>
> $window: 23 | start: 221 | stop: 470 > g1 g2 g3 > g1 NA 0.5473118 0.3776754 > g2 NA NA 1.9126611 > g3 NA NA NA > >$window: 24 | start: 231 | stop: 480
>    g1        g2        g3
> g1 NA 0.3456384 0.3320625
> g2 NA        NA 1.7731634
> g3 NA        NA        NA
>
> $window: 25 | start: 241 | stop: 490 > g1 g2 g3 > g1 NA 0.5450731 0.3883475 > g2 NA NA 1.8579286 > g3 NA NA NA > >$window: 26 | start: 251 | stop: 495
>    g1        g2        g3
> g1 NA 0.4465384 0.4984407
> g2 NA        NA 1.6798012
> g3 NA        NA        NA

# It may be informative to create an animation [not shown here because it displays in the Viewer]
# MRN_ani <- mrn_plot(layers = layers, win = 250, step = 10, createAnimation = TRUE)

# The animation is stored in an output field, but is also saved as a .gif, see the man pages for more options.
# MRN_ani\$MRNanimationGG