Unidimensional and Multidimensional Methods for Recurrence Quantification Analysis
The crqa R package allows users to conduct a wide range of recurrence-based analyses on single (e.g., auto-recurrence) and multivariate time series (e.g., multidimensional cross-recurrence quantification), examine coupling properties underlying leader-follower relationships (i.e., diagonal-profile methods), as well as track the evolution of recurrence rate over the time course (i.e., windowed methods).
# You can install the latest version of crqa on CRAN by running:
install.packages("crqa")
# Or for the development version from GitHub:
# install.packages("devtools")
# devtools::install_github("morenococo/crqa")
install.packages("crqa", repos = c('https://morenococo.r-universe.dev', 'https://cloud.r-project.org'))
5–15× faster, lower memory. A fused C++ (Rcpp) kernel replaces the old R pipeline for the three most common distance metrics (Euclidean, maximum, Manhattan). It runs a single pass over the phase-space matrix and never allocates the full N×N distance or recurrence matrix — memory scales as O(N + nnz) instead of O(N²). The practical limit grew from ~10 000 to ≥ 20 000 data points within a 30-second budget on a laptop.
OpenMP parallelism. The fused kernel runs in parallel across CPU cores out of the box. Sparse recurrence regimes (the typical CRQA case) scale near-linearly with the number of cores. Thread count is controlled via the standard OMP_NUM_THREADS environment variable.
Approximative RQA for very long series. method = "aRQA" in crqa() dispatches to a phase-space histogram algorithm (Schultz et al. 2015) that computes RR, DET and L in O(N) memory. Designed for N ≥ 10 000 where exact computation is expensive.
Theiler-aware recurrence rate. All functions now accept rr_denom = "valid" (default "full" for backward compatibility), which excludes Theiler-blanked and side-masked cells from the RR denominator — making RR internally consistent with DET and ENTR.
Global normalisation fix. wincrqa(), windowdrp() and piecewiseRQA() now apply normalize once to the full series before windowing, so all windows share the same scale reference. Previously normalisation was applied independently per window, shifting the effective threshold across windows.
New functions. aRQA(), line_stats(), theiler_exclusion(), rosslerattractor().
crqa comes with bundled datasets that can be used to explore the different functions.
data(crqa) # loads eyemovement, handmovement, and text datasets
Run auto-recurrence on a nursery rhyme (“The wheels on the bus”, 120 words):
delay = 1; embed = 1; rescale = 0; radius = 0.0001
normalize = 0; mindiagline = 2; minvertline = 2
tw = 1; whiteline = FALSE; recpt = FALSE
side = “both”; method = “rqa”; metric = “euclidean”
datatype = “categorical”
ans = crqa(ts1 = text, ts2 = text,
delay = delay, embed = embed, rescale = rescale,
radius = radius, normalize = normalize,
mindiagline = mindiagline, minvertline = minvertline,
tw = tw, whiteline = whiteline, recpt = recpt,
side = side, method = method, metric = metric,
datatype = datatype)
print(ans[!names(ans) %in% “RP”])
Eye-tracking data from a joint task: 2,000 observations of six screen locations looked at by a narrator and a listener.
narrator = eyemovement$narrator
listener = eyemovement$listener
delay = 1; embed = 1; rescale = 0; radius = 0.001
normalize = 0; mindiagline = 2; minvertline = 2
tw = 0; whiteline = FALSE; recpt = FALSE; side = “both”
method = “crqa”; metric = “euclidean”; datatype = “categorical”
ans = crqa(ts1 = narrator, ts2 = listener,
delay = delay, embed = embed, rescale = rescale,
radius = radius, normalize = normalize,
mindiagline = mindiagline, minvertline = minvertline,
tw = tw, whiteline = whiteline, recpt = recpt,
side = side, method = method, metric = metric,
datatype = datatype)
print(ans[!names(ans) %in% “RP”])
Extract the diagonal cross-recurrence profile (DCRP) to quantify leader–follower dynamics across a range of lags:
timecourse = round(seq(-3300, 3300, 33) / 1000, digits = 2)
res = drpfromts(ts1 = narrator, ts2 = listener,
windowsize = 100,
delay = 1, embed = 1, rescale = 0,
radius = 0.001, normalize = 0,
mindiagline = 2, minvertline = 2,
tw = 0, whiteline = FALSE, recpt = FALSE,
side = “both”, method = “crqa”,
metric = “euclidean”, datatype = “categorical”)
profile = res$profile * 100 # recurrence rate in %
plot(timecourse, profile, type = “l”, lwd = 2.5,
xlab = “Lag (seconds)”, ylab = “Recurrence Rate %”)
Hand-movement data from a LEGO joint construction task (5,799 observations, two position channels per person):
P1 = cbind(handmovement$P1_TT_d, handmovement$P1_TT_n)
P2 = cbind(handmovement$P2_TT_d, handmovement$P2_TT_n)
delay = 5; embed = 2; rescale = 0; radius = 0.1
normalize = 0; mindiagline = 10; minvertline = 10
tw = 0; whiteline = FALSE; recpt = FALSE; side = “both”
method = “mdcrqa”; metric = “euclidean”; datatype = “continuous”
ans = crqa(ts1 = P1, ts2 = P2,
delay = delay, embed = embed, rescale = rescale,
radius = radius, normalize = normalize,
mindiagline = mindiagline, minvertline = minvertline,
tw = tw, whiteline = whiteline, recpt = recpt,
side = side, method = method, metric = metric,
datatype = datatype)
print(ans[!names(ans) %in% “RP”])