H.1 Assignment: Cross Recurrence Quantification Analysis

| Jump to assignment |

  1. Create two variables for CRQA analysis:
y1 <- sin(1:900*2*pi/67)
y2 <- sin(.01*(1:900*2*pi/67)^2)
  1. You have just created two sine waves. We’ll examine if and how they are coupled in a shared phase space. As a first step plot them.

  2. Find an embedding delay (using mutual information) and an embedding dimension (if you calculate an embedding dimension using package fractal for each signal seperately, as a rule of thumb use the highest embedding dimension you find in further analyses).

# General settings for `crqa()`
par0 <- list(rescale = 0,
             normalize = 0,
             mindiagline = 2,
             minvertline = 2,
             tw = 0,
             whiteline = FALSE,
             recpt = FALSE,
             side = "both",
             checkl = list(do = FALSE, thrshd = 3, datatype = "categorical",pad = TRUE)
             )

# Settings for `optimizeParam()`
par <- list(lgM =  20, steps = seq(1, 6, 1),
           radiusspan = 100, radiussample = 40,
           normalize = par0$normalize, 
           rescale = par0$rescale, 
           mindiagline = par0$mindiagline, minvertline = par0$minvertline,
           tw = par0$tw, 
           whiteline = par0$whiteline, 
           recpt = par0$recpt, 
           fnnpercent = 10, typeami = "mindip")
  1. We can now create a cross recurrence matrix. Fill in the values you decided on. You can choose a radius automatically, look in the crqa manual.

Note: There is no rescaling of data, the sines were created in the same range. You can plot a matrix using image(). You could also check package nonlinearTseries. If you sourced the nlRtsa library, use plotRP.crqa()

  • Get the optimal parameters using a radius which will give us 2%-5% recurrent points.
(ans <- optimizeParam(ts1 = y1, ts2 = y2, par = par, min.rec = 2, max.rec = 5))
$radius
[1] 0.2128611

$emddim
[1] 2

$delay
[1] 15
  1. Run the CRQA.
crqaOutput <- crqa(ts1 = y1, ts2 = y2,  
                  delay = ans$delay, 
                  embed = ans$emddim, 
                  radius = ans$radius, 
                  normalize = par0$normalize,
                  rescale = par0$rescale, 
                  mindiagline = par0$mindiagline, minvertline = par0$minvertline,
                  tw = par0$tw, 
                  whiteline = par0$whiteline, 
                  recpt = par0$recpt, 
                  side = par0$side, checkl = par0$checkl
                  )
  1. Can you understand what is going on? Explain the the lack of recurrent points at the beginning of the time series.
plotRP.crqa(crqaOutput)

  1. Examine the synchronisation under the diagonal LOS. Look in the manual of crqa or Coco & Dale (2014).

  2. Make a plot of the diagonal profie. How far is the peak in RR removes from 0 (Line of Synchronisation)?

This is based on the timeseries (without embedding)

win <- 50
(res <-  drpdfromts(y1, y2, ws = win, datatype = 'continuous', radius = ans$radius))[2:3]
$maxrec
[1] 0.3080357

$maxlag
[1] 55
plot(-win:win,res$profile,type = "l", lwd = 5, xlab = "Delays", ylab = "Recurrence")
abline(v=0)

To get the diagonal profile from the recurrence plot, use spdiags().

dprofile <- spdiags(crqaOutput$RP)
plot(-win:win,colMeans(dprofile$B[, between(dprofile$d, -win,win)]), type="l", 
     lwd = 5, xlab = "Delays", ylab = "Recurrence")
abline(v=0)

  1. Perform the same steps with a shuffled version (or use surrogate analysis!) of the data of timeseries \(y1\). You can use the embedding parameters you found earlier.