Tutorial: Quantifying Interpersonal Dynamics using CRQA


  1. a basic understanding of (C)RQA and its applicability to research in cognitive science, and apply it to interpersonal dynamics.
  2. a primer, hands-on, of (C)RQA to behavioral data using the crqa() package in R.


Interpersonal dynamics has become a growing topic in cognitive science because of its implications for our ‘social cognitive system.’ Advances in this research topic have been possible through the application of statistical methods providing a quantification for the dynamic structure of cognitive responses observed when individuals interact. Recurrence Quantification Analysis (RQA) is one such framework, as it makes possible to quantify how, and the extent to which, a signal is revisiting a similar state in time. In the context of the dyad, it can be used to reveal the phases during which interlocutor responses synchronize, and the strength of such synchronization.

Focus of this source code

  1. R basics: Brief primer on Markdown, recap of basic R functions, and pre-processing of raw data
  2. Categorical data: Auto and Cross Recurrence
  3. Continuous data: Auto and Cross Recurrence

R basics

Very brief primer on Markdown and Rknit

Paragraphing: continuous = sign makes the title, single # makes a section, a double ## makes a paragraph

This is a quote

To compile html make, library(knitr) must be loaded, then either press Knit HTML (R-Studio) or knit('tutorial_knitted.Rmd'), with command line.

Hyperlinks: see the borrowed R Markdown example from here.

italic and bold.

A code chunk opens with three backticks and an r in curly bracket, and closes using three backticks.

## [1] 2

inline code is enclosed as 3.1415927.

Code can also be nested into numbered (or not) lists:

  1. String split

    strsplit('hello indented world', ' ')[[1]]
    ## [1] "hello"    "indented" "world"
  2. Substitute strings

    a = 'hate'
    ## [1] "hate"
    sub('hate', 'love', a)
    ## [1] "love"

Equations are enclosed using the dollar sign: : \(\alpha+\beta=\gamma\)

Plots are simply embedded in the code chunk

hist(rnorm(10000), breaks = 500)

With these basics, you can write almost any R document in MarkDown.

Load libraries and data for the tutorial.

## Loading required package: Matrix
## Loading required package: tseriesChaos
## Loading required package: deSolve
## Loading required package: fields
## Loading required package: spam
## Loading required package: grid
## Spam version 1.0-1 (2014-09-09) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##     backsolve, forwardsolve
## Loading required package: maps
## Attaching package: 'fields'
## The following object is masked from 'package:plyr':
##     ozone
## Loading required package: plot3D
## Loading required package: pracma
## Attaching package: 'pracma'
## The following object is masked from 'package:deSolve':
##     rk4
## The following objects are masked from 'package:Matrix':
##     expm, lu, tril, triu
# projects/papers/RickMorenoRQA/tutorial_CRQA_cogsci2015/')  # change to your own path (after download)  
source('lorenzattractor.R') # unless you have the last version of crqa() 1.0.6
load('alldata.Rdata') # load the Rdata object
# take out of the Rdata object all datasets to be used

emotData = alldata$emotData # emotional states (parent, adolescent)
jepsen =  alldata$jepsen # lyrics from carly jepsen (raw text)
obama = alldata$obama # speech from barak obama (raw text)
eyemovement = alldata$eyemovement # from Richardson & Dale, 2005 (eye movements, speaker, listener)
bodymovement =  alldata$bodymovement # from Paxton & Dale, 2013 (bodily motion of conversants)

Categorical events over a time-course

Many transcription and video-coding systems (e.g., ELAN, etc.) utilize an event-interval coding scheme that specifies a given event and its start / end times.

Let’s consider, for example, data about emotional dynamics between parent and child during an interaction (Dr. Alexandra Main’s data, UC Merced); refer to emotData.

Initially, the emotional time line of the parent looked like this:

00:00:44:23,00:00:50:09,P Neutral,,,

The timestamps from start and end of the emotional state, the identifier for (P)arent and (A)dolescent, and the actual state. The data do not start at time 0s, but at 44s, and “parent humor” is the emotional event taking place in the interation.

The first thing one must do to start recurrence is to obtain a new data format which encodes events with regular time sampling.

Put it differently, recurrence assumes that you are inputting time series which represent regularly sampled measurements (either scalars or nominal identifiers).

This is just about the sole “pre-analysis” assumption that recurrence requires before stepping into its magical world.

(Minor note: In a sense, recurrence also requires the assumption of “exclusive” codes. In any given time series, we can only have only one measure or event at a given time; which means that all time points have a discrete code. However, in special cases, such as non-events, recurrence might not require it.)

So, we reshaped the data into a two-column second-by-second “time series”.

##   s    parent adolescent
## 1 0 P Neutral    A Humor
## 2 1 P Neutral    A Humor
## 3 2 P Neutral    A Humor
## 4 3 P Neutral    A Humor
## 5 4 P Neutral    A Humor
## 6 5 P Neutral    A Humor

Now the left column shows time in seconds, starting from 0s. The file is tab-delimited and the second and third columns reflect the emotional events taking place in each person.

Pre-processing is not over. In fact, we need to code the data seeking for the “recurring states.”

Recode the data into event states

In order to determine state matches, we recode the time series by stripping the “A” and “P” prefixes so that we can play with the time series and leverage matching strings.

parent = as.character(emotData[,2]) # making it characters rather than factor list
adoles = as.character(emotData[,3]) 
# let's get emotion state out (string), stripping prefixes
for (i in 1:length(parent)) { 
  parent[i] = unlist(strsplit(parent[i]," "))[2] # space delimiter
  adoles[i] = unlist(strsplit(adoles[i]," "))[2]
## [1] "Neutral" "Neutral" "Neutral" "Neutral"
## [1] "Humor" "Humor" "Humor" "Humor"

Now let’s get the set of unique codes across parent and adolescent.

uniques = unique(c(parent,adoles))
##  [1] "Neutral"        "Tension"        "Domineering"    "Anger"         
##  [5] "Defensiveness"  "Elaboration"    "Backchannels"   "Sadness"       
##  [9] "Open-ended"     "Mirroring"      "Humor"          "Positive"      
## [13] "Direct"         "Identification" "Whining"        "Belligerence"  
## [17] "Sentence"

There are 17 unique codes. We can use the ordering of elements in this list to consistently renumber each of the time series we’re analyzing (parent / adolescent). Let’s use a lame little loop to do so.

pix = (1:length(parent))*0 # pix = parent's index in uniques
aix = (1:length(parent))*0 # default to 0's; aix = adolescent's
for (i in 1:length(parent)) {
  # which returns the index in uniques[j] of the parent[i] string
  pix[i] = which(uniques==parent[i]) 
  aix[i] = which(uniques==adoles[i]) 
} # below is a more efficient example of this with "apply"
##  [1] 1 1 1 1 1 1 2 2 2 2
##  [1] 11 11 11 11 11 11 11 11 11 11

Another solution would be to use the function checkts(), which takes as input two time-series of categorical nature, check whether their length is equal, and accept or reject according to a user specified threshold. Series can be made equal by padding the shortest, or chopping the longest. Finally, categorical states are recoded into numeric identifiers.

ans = checkts(parent, adoles, datatype = 'categorical', thrshd = 5, pad = FALSE)
recoded = ans[[1]]

pix = recoded[,1]
aix = recoded[,2]
# was the difference between the series smaller than the threshold?
## [1] TRUE

We can plot these changes now, if we want to eyeball the tendency for emotions to change during their interaction.

First we create an x-axis to unfold the states (realtime):

realtime = emotData$s
indnorm = seq(1, max(realtime), 42)
indnorm = c(indnorm, max(realtime))

Generate as many colors as there are unique states (look at the function rainbow())

cols = rainbow(length(uniques)) 

If we wish to plot the states with different colors plot(), we can implement a loop that goes through all points of the two series, identify their state, and color code accordingly on the plot as points().

par(mar = c(4, 5, 2, 2), font.lab = 2, cex.lab = 2.5, font.axis = 2, cex.axis = 2.2)

plot(realtime, rep(1, length(realtime)),
     lwd = 3, cex = .8,
     col = "gray", xaxt = "n", yaxt = "n",
     ylim = c(0, 3.5), xlab = "Time-Course (seconds)",
     ylab = "Interlocutor",
     type = "n", pch = 19

mtext("Parent", at = 1, cex = 2, font = 2, line = 1, side = 2)
mtext("Adoles", at = 2.2, cex = 2, font = 2, line = 1, side = 2)

r = 2
for (r in 1:length(parent)){
    curP = as.character(parent[r])
    curS = as.character(adoles[r])

    indcP = which(uniques == curP)
    indcS = which(uniques == curS)

    points(realtime[r], 1, col = cols[indcP], cex = 2, pch = 16)
    points(realtime[r], 2, col = cols[indcS], cex = 2, pch = 16)


mtext(seq(1, max(realtime), 1)[indnorm], at = seq(1, max(realtime), 1)[indnorm], cex = 2, font = 2, line = 1, side = 1)

legend("top", uniques, col = cols, ncol = 3, cex = 1, bty = "n",
       pch = 16)