Windowed Cross-Correlation: An Illustration

One therapy dyad, sixteen sessions, and 199 coded moments

library(tidyverse)
library(patchwork)
library(ggrepel)
# wcc, lme4 and lmerTest are called with :: where used, so it is always visible
# where a computation leaves the tidyverse.

# The parameter sweep in Step 2 takes about fifteen minutes. Results are
# written to disk the first time and read back afterwards. Delete the
# directory to force a recomputation.
DERIVED <- here::here("tutorials", "_derived", "08i")
dir.create(DERIVED, recursive = TRUE, showWarnings = FALSE)

cached <- function(name, expr) {
  path <- file.path(DERIVED, paste0(name, ".rds"))
  if (file.exists(path)) return(readRDS(path))
  value <- force(expr)
  saveRDS(value, path)
  value
}

PERSON <- c(patient = "#1f6f8b", therapist = "#c05746")
TYPES  <- c("Reflection", "Action", "Protest", "Reconceptualization")
TYPE   <- c(Reflection = "#7fb3d5", Action = "#f0b27a",
            Protest = "#c39bd3", Reconceptualization = "#76d7c4")

theme_set(
  theme_minimal(base_size = 11) +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(face = "bold", size = 11),
          plot.subtitle = element_text(colour = "grey35", size = 9),
          plot.caption = element_text(colour = "grey35", size = 8.5, hjust = 0),
          strip.text = element_text(face = "bold", size = 9.5, hjust = 0),
          legend.position = "top")
)

# --- EDIT THIS ONE LINE ------------------------------------------------------
CASE_DIR <- file.path(
  "/Users/nilamram/Library/CloudStorage/GoogleDrive-nilamram@stanford.edu",
  "Shared drives/ChangeLab/Studies_DataSets/Psychophysiological Synchrony",
  "Workshop/data_check/kleinbub_case"
)

---

The question

A correlation between two people’s physiological series answers one question — do they rise and fall together over the recording as a whole — and answers it with a single number. Two things are lost. The first is that the association changes over the course of an interaction, so a single number describes no particular moment. The second is that two people rarely respond simultaneously; one leads and the other follows, by an interval that itself shifts as the roles change. This illustration recovers both: it estimates the association in short overlapping windows, and within each window allows the two series to slide against each other to find the alignment at which they match best. The result is read at 199 moments that a clinical coder, working from the transcript and blind to the physiology, had already marked as meaningful.

The data

One patient–therapist dyad in psychotherapy, recorded over 16 sessions of roughly 45 minutes each. Skin conductance was sampled from both people simultaneously at 100 Hz. Fifteen of the sessions additionally carry coded Innovative Moments — stretches of talk in which the patient departs from the problematic pattern that brought them to therapy. Each is marked with a start time, an end time and a type.

FS <- 10  # Hz, as analysed

signals  <- readRDS(file.path(CASE_DIR, "kleinbub_case_10hz.rds"))
episodes <- readRDS(file.path(CASE_DIR, "kleinbub_episodes.rds"))

c(sessions = n_distinct(signals$session),
  rows     = nrow(signals),
  episodes = nrow(episodes))

## sessions     rows episodes 
##       16   465276      199

Those two objects were built once by a preparation script. Two decisions were taken there and both belong in the open.

The series were decimated from 100 Hz to 10 Hz. Reducing a sampling rate is two operations whose order cannot be swapped: low-pass filter first, then keep every \(q\)-th sample. Omitting the filter does not remove content above the new Nyquist frequency — it folds that content down onto lower frequencies, where it becomes indistinguishable from signal. The script calls signal::decimate(x, q = 10, ftype = "fir", n = 30), which performs both steps.

The episode times were read from the per-session marker files, not the pooled file supplied beside them: the two disagree by one episode, and only the per-session files carry the session identifier unambiguously.

episodes |>
  summarise(n = n(), shortest = min(duration), median = median(duration),
            longest = max(duration), `n >= 60 s` = sum(duration >= 60),
            .by = type) |>
  arrange(desc(n))

type	n	shortest	median	longest	n >= 60 s
Reflection	116	1	16.0	233	14
Action	45	2	16.0	109	5
Protest	30	1	15.0	127	2
Reconceptualization	8	52	224.5	478	7

The median episode lasts 17 seconds. The types differ sharply in length — Reconceptualization has a median of 224 s against roughly 16 s for the other three — so any comparison between types is partly a comparison between durations.

---

Step 1 — Look at the signals with the episodes on them

label_for <- function(s) {
  sprintf("Session %d — %d coded episodes", s, sum(episodes$session == s))
}
SHOW <- c(1, 5)

series_show <- signals |>
  filter(session %in% SHOW) |>
  pivot_longer(c(patient, therapist), names_to = "who", values_to = "sc") |>
  mutate(panel = map_chr(session, label_for))

episodes_show <- episodes |>
  filter(session %in% SHOW) |>
  mutate(panel = map_chr(session, label_for))

ggplot() +
  geom_rect(data = episodes_show,
            aes(xmin = start / 60, xmax = end / 60,
                ymin = -Inf, ymax = Inf, fill = type), alpha = .6) +
  geom_line(data = series_show,
            aes(t / 60, sc, colour = who), linewidth = .25) +
  facet_wrap(~ panel, ncol = 1, scales = "free") +
  scale_fill_manual(values = TYPE, limits = TYPES, name = NULL, drop = FALSE) +
  scale_colour_manual(values = PERSON, name = NULL) +
  labs(x = "minutes into session", y = "skin conductance (µS)",
       caption = paste("Episode windows are drawn at true width; the median",
                       "episode is 17 s, so most appear as slivers.")) +
  guides(colour = guide_legend(order = 1, override.aes = list(linewidth = 1.2)),
         fill = guide_legend(order = 2, override.aes = list(alpha = .85)))

Three features matter later. The patient’s signal carries frequent sharp rises — discrete skin conductance responses, each climbing over one to three seconds and recovering over several more — while the therapist’s is smoother. Both drift upward across the session, which a whole-recording correlation would read as agreement. And the episodes are short against a 45-minute recording.

---

Step 2 — Choose the window and the maximum lag

Windowed cross-correlation has two parameters that determine what it can see. The window sets how much of the recording contributes to each estimate, and therefore the slowest fluctuation the method can register: a rhythm longer than the window appears within it only as a trend. The maximum lag sets how far apart in time two responses may be and still count as related.

Two published sources disagree about what these should be for skin conductance, and the disagreement is instructive because it is a disagreement about method.

Behrens and colleagues (2020) set them from the physiology. Reviewing studies of the time course of individual physiological responses, they found that “most studies used response durations between 5 and 10 seconds for both heart rate and skin conductance level,” and chose a window of 8 s with a maximum lag of 4 s — a lag at half the window.

Behrens and colleagues (2026), largely the same authors, set them statistically instead, sweeping a grid of configurations and scoring each by how well it separated real dyads from surrogate pairings. For skin conductance they recommend “a small window size and a large maximum lag that is around four times the window size,” their best cell being a 6 s window with a 24 s lag.

These are close to opposites. Rather than choose between them we can run the sweep here, because the criterion is reproducible: compute the synchrony measure for every configuration, for real pairings and for surrogate pairings, and ask which configurations separate them.

WIN_GRID <- c(5, 6, 8, 10, 12, 15, 20, 25)
LAG_GRID <- c(4, 6, 8, 12, 16, 20, 24)

wcc_pair <- function(a, b, winSec, lagSec) {
  wMax <- winSec * FS; tMax <- lagSec * FS
  wInc <- max(1, round(wMax / 10))     # window increment, after Behrens
  n <- min(length(a), length(b))
  M <- wcc::wccCalc(a[1:n], b[1:n], wMax = wMax, tMax = tMax,
                    wInc = wInc, tInc = 1, method = "cumc")
  P <- wcc::wccPeakPick(M, Lsize = 8, pspan = 0.25, type = "Max")
  tibble(mean_r    = mean(P$maxValue, na.rm = TRUE),
         mean_absr = mean(abs(M), na.rm = TRUE),
         pct_na    = 100 * mean(is.na(P$maxIndex)))
}

patient   <- map(1:16, \(s) filter(signals, session == s)$patient)
therapist <- map(1:16, \(s) filter(signals, session == s)$therapist)

grid_score <- cached("grid_score", {
  set.seed(2)
  mismatches <- map_dfr(1:16, \(s) tibble(p = s, t = sample(setdiff(1:16, s), 3)))
  expand_grid(winSec = WIN_GRID, lagSec = LAG_GRID) |>
    pmap_dfr(\(winSec, lagSec) {
      real <- map_dfr(1:16, \(s)
        wcc_pair(patient[[s]], therapist[[s]], winSec, lagSec)) |> mutate(kind = "real")
      surr <- map2_dfr(mismatches$p, mismatches$t, \(p, t)
        wcc_pair(patient[[p]], therapist[[t]], winSec, lagSec)) |> mutate(kind = "surrogate")
      bind_rows(real, surr) |>
        summarise(real = mean(mean_r[kind == "real"]),
                  surrogate = mean(mean_r[kind == "surrogate"]),
                  t = t.test(mean_r[kind == "real"],
                             mean_r[kind == "surrogate"])$statistic,
                  pct_na = mean(pct_na)) |>
        mutate(winSec = winSec, lagSec = lagSec)
    })
})

marks <- tribble(~winSec, ~lagSec, ~label,
                 10,  4, "Kleinbub",
                  8,  4, "Behrens 2020",
                  6, 24, "Behrens 2026\nrecommended")

panel_a <- ggplot(grid_score, aes(factor(lagSec), factor(winSec), fill = t)) +
  geom_tile(colour = "white", linewidth = .6) +
  geom_text(aes(label = sprintf("%.1f", t)), size = 2.9, colour = "grey15") +
  geom_point(data = marks, aes(factor(lagSec), factor(winSec)), inherit.aes = FALSE,
             shape = 21, size = 11, stroke = 1.1, colour = "black", fill = NA) +
  scale_fill_gradient2(low = "#2c5f7c", mid = "grey96", high = "#a8322d",
                       midpoint = 0, name = "t") +
  labs(title = "A. Separation of real from surrogate pairings, 56 configurations",
       subtitle = paste("t for mean peak r, real dyads (n = 16) against",
                        "cross-session surrogates (n = 48)."),
       x = "maximum lag (s)", y = "window size (s)") +
  theme(legend.position = "right")

panel_b <- ggplot(grid_score, aes(pct_na, t)) +
  geom_smooth(method = "lm", se = FALSE, colour = "grey75", linewidth = .5) +
  geom_point(aes(colour = factor(lagSec)), size = 2) +
  geom_text_repel(data = left_join(marks, grid_score, by = c("winSec", "lagSec")),
                  aes(label = label), size = 3, colour = "grey20",
                  box.padding = .6, seed = 1) +
  scale_colour_viridis_d(option = "mako", end = .85, name = "lag (s)") +
  labs(title = "B. The ranking tracks how much data is discarded",
       subtitle = sprintf("r = %+.2f between separation and the proportion of windows with no peak.",
                          cor(grid_score$t, grid_score$pct_na)),
       x = "% of windows with no peak", y = "t") +
  theme(legend.position = "right")

panel_a / panel_b + plot_layout(heights = c(1.25, 1))

Every configuration separates the two groups in the right direction, but the gradient runs opposite to the 2026 recommendation. Short lags do better than long ones — mean \(t\) falls from 2.68 at a 4 s lag to 1.11 at 24 s — and larger windows do better than small, rising from 0.74 at 5 s to 2.38 at 25 s. The recommended 6 s / 24 s cell ranks 54th of 56.

Panel B shows why that ranking cannot be taken at face value.

grid_score |>
  arrange(desc(t)) |>
  mutate(rank = row_number()) |>
  filter((winSec == 10 & lagSec == 4) | (winSec == 8 & lagSec == 4) |
         (winSec == 6 & lagSec == 24) | rank <= 3) |>
  transmute(rank, `window (s)` = winSec, `lag (s)` = lagSec,
            t = round(t, 2), `% no peak` = round(pct_na, 1)) |>
  arrange(rank)

rank	window (s)	lag (s)	t	% no peak
1	25	4	3.76	40.8
2	20	4	3.43	37.1
3	12	4	3.29	26.3
6	10	4	2.88	22.1
16	8	4	2.10	17.1
54	6	24	0.46	0.0

The configurations that separate best are the configurations that throw the most data away. Peak-picking returns a missing value when no qualifying peak exists, and short lags produce many such windows; the best cell discards 41% of them. Separation and discard correlate at \(r = +0.83\) across the grid, so the sweep is ranking selectivity at least as much as sensitivity. A criterion that rewards discarding data is not a criterion one should follow to its maximum.

Taking both panels together, we use a 10-second window with a maximum lag of 4 seconds, stepped every 2 seconds. The window sits inside the 5–10 s range that the 2020 paper derived from the time course of a skin conductance response; the lag is half the window, as there; the configuration ranks 6th of 56 on separation while discarding less than the five above it; and it is the configuration under which these particular recordings have been analysed before.

---

Step 3 — Compute the surface and pick the peaks

Within each window the two series are slid against one another and a correlation computed at every offset, giving a matrix with one row per window and one column per lag. For window \(w\) and lag \(\ell\),

\[r_{w\ell} \;=\; \mathrm{cor}\!\left(A_{[w]},\; B_{[w + \ell]}\right)\]

and the peak for that window is the correlation closest to zero lag that qualifies as a local maximum,

\[\hat{\ell}_w \;=\; \arg\max_{\ell} \; r_{w\ell}, \qquad \hat{r}_w \;=\; r_{w\hat{\ell}_w}\]

WIN_S <- 10; LAG_S <- 4; STEP_S <- 2
WMAX <- WIN_S * FS; TMAX <- LAG_S * FS; WINC <- STEP_S * FS

peaks_for <- function(a, b) {
  n <- min(length(a), length(b))
  M <- wcc::wccCalc(a[1:n], b[1:n], wMax = WMAX, tMax = TMAX,
                    wInc = WINC, tInc = 1, method = "cumc")
  P <- wcc::wccPeakPick(M, Lsize = 8, pspan = 0.25, type = "Max")
  hi <- (WMAX + TMAX) + (seq_len(nrow(M)) - 1) * WINC   # last sample of window
  tibble(t         = (hi - WMAX / 2 - 1) / FS,          # window centre
         peak_r    = P$maxValue,
         peak_lag  = P$maxIndex / 2 / FS,               # maxIndex is HALF-samples
         mean_absr = rowMeans(abs(M), na.rm = TRUE))
}

peaks <- cached("peaks", {
  map_dfr(1:16, \(s) {
    d <- filter(signals, session == s)
    peaks_for(d$patient, d$therapist) |> mutate(session = s)
  })
})

peaks |>
  summarise(windows = n(),
            `mean peak r` = round(mean(peak_r, na.rm = TRUE), 3),
            `mean peak lag (s)` = round(mean(peak_lag, na.rm = TRUE), 3),
            `SD of peak lag (s)` = round(sd(peak_lag, na.rm = TRUE), 3),
            `% with no peak` = round(100 * mean(is.na(peak_r)), 1))

windows	mean peak r	mean peak lag (s)	SD of peak lag (s)	% with no peak
23143	0.53	0.033	1.74	21.1

Two features of that output need explaining before anything is built on it.

The lag is returned in half-samples. wccPeakPick smooths each row across the lag axis onto a grid of twice the original resolution before searching, so its index must be divided by two as well as by the sampling rate. Read without that correction every lag would be reported at twice its true value.

One window in five yields no peak at all. This is deliberate in the algorithm, not a failure of it: a peak must have eight descending values on each side and must not sit at the edge of the searched range, and when nothing qualifies a missing value is returned. The figure below shows what that looks like.

S <- 5; T0 <- 600; T1 <- 960
d5 <- filter(signals, session == S)
M5 <- wcc::wccCalc(d5$patient, d5$therapist, wMax = WMAX, tMax = TMAX,
                   wInc = WINC, tInc = 1, method = "cumc")
P5 <- wcc::wccPeakPick(M5, Lsize = 8, pspan = 0.25, type = "Max")
tc <- ((WMAX + TMAX) + (seq_len(nrow(M5)) - 1) * WINC - WMAX / 2 - 1) / FS
pk <- tibble(t = tc, lag = P5$maxIndex / 2 / FS, r = P5$maxValue,
             found = !is.na(P5$maxIndex))

sel  <- which(tc >= T0 & tc <= T1)
surf <- expand_grid(i = sel, lag = seq(-TMAX, TMAX) / FS) |>
  mutate(t = tc[i], r = as.numeric(t(M5[sel, , drop = FALSE])))
pks  <- filter(pk, t >= T0, t <= T1)

p_surface <- ggplot(surf, aes(t / 60, lag, fill = r)) +
  geom_raster(interpolate = TRUE) +
  geom_rect(data = filter(pks, !found), inherit.aes = FALSE,
            aes(xmin = (t - STEP_S / 2) / 60, xmax = (t + STEP_S / 2) / 60,
                ymin = -LAG_S, ymax = LAG_S), fill = "white", alpha = .75) +
  geom_point(data = filter(pks, found), inherit.aes = FALSE,
             aes(t / 60, lag), colour = "#1b9e77", size = .5) +
  geom_hline(yintercept = 0, colour = "grey20", linewidth = .3, linetype = 2) +
  scale_fill_gradient2(low = "#2c5f7c", mid = "grey97", high = "#a8322d",
                       midpoint = 0, limits = c(-1, 1), name = "r") +
  scale_y_continuous(breaks = seq(-4, 4, 2), expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0)) +
  labs(title = "A. The correlation surface, six minutes of session 5",
       subtitle = "Green dots are the picked peak. White bands are windows where no peak qualified.",
       x = NULL, y = "lag (s)   ↑ therapist leads")

p_track <- ggplot(pks, aes(t / 60, r)) +
  geom_hline(yintercept = 0, colour = "grey70", linewidth = .3) +
  geom_point(colour = "#1b9e77", size = .6) +
  coord_cartesian(ylim = c(-1, 1)) + scale_x_continuous(expand = c(0, 0)) +
  labs(subtitle = "B. The peak correlation that survives. Gaps are the white bands above.",
       x = "minutes into session", y = "peak r")

SMOOTH <- wcc::wccSmoothMatrix(ncol(M5), 0.25)
gridx  <- (seq_len(2 * ncol(M5) - 1) - ncol(M5)) / 2 / FS
idx    <- c(head(sel[pk$found[sel]], 2), head(sel[!pk$found[sel]], 2))
labs_v <- make.unique(sprintf("t = %.2f min — %s", tc[idx] / 60,
                              if_else(pk$found[idx], "peak found", "no peak")), sep = "  ")
examples <- map2_dfr(idx, labs_v, \(i, lb)
  tibble(label = lb, lag = gridx, r = as.numeric(M5[i, ] %*% t(SMOOTH)),
         pk_lag = pk$lag[i], pk_r = pk$r[i], found = pk$found[i])) |>
  mutate(label = factor(label, levels = labs_v))

p_slices <- ggplot(examples, aes(lag, r)) +
  geom_hline(yintercept = 0, colour = "grey85", linewidth = .3) +
  geom_line(linewidth = .45, colour = "grey25") +
  geom_point(data = distinct(filter(examples, found), label, pk_lag, pk_r),
             aes(pk_lag, pk_r), colour = "#1b9e77", size = 2.4) +
  facet_wrap(~ label, nrow = 1) +
  labs(title = "C. Why a window does or does not yield a peak",
       subtitle = "The smoothed correlation across lags within a single window.",
       x = "lag (s)", y = "r")

(p_surface / p_track / p_slices) + plot_layout(heights = c(1.5, .7, 1.1))

How to read the surface

Panel A is the central object of the method, and it repays a slow look.

Each vertical stripe is one window — one ten-second stretch of the session, centred at the time given on the horizontal axis. Reading left to right is reading forward through the session.

Within a stripe, height is the lag. At the centre line the two series are compared as recorded, simultaneously. Higher up, the therapist’s series has been shifted so that it is compared against the patient’s later responses — that is, the therapist is treated as leading. Lower down, the patient leads. The vertical extent of the plot is the ±4 s we allowed.

Colour is the correlation at that lag in that window. Red is positive, the two rising and falling together; blue is negative, one rising while the other falls; pale is no relationship. A single stripe is therefore the whole cross-correlation function for that moment, drawn as a column of colour — exactly the curves plotted in panel C, turned on their side.

The green dot is the chosen peak: the reddest point in that stripe, subject to being a genuine interior maximum. Where the dots trace a smooth path, the lead–follow relationship is holding steady. Where they jump, it has reorganised. Where the stripe is white, no peak qualified.

Read that way, two things are visible. The surface is organised into broad diagonal bands several seconds wide that drift slowly up or down across the lag axis. A band drifting upward means the therapist is progressively leading by more; drifting downward, by less. These bands take tens of seconds to cross the ±4 s range, which is longer than the range itself.

That is the explanation for the white stripes, and panel C confirms it. The two windows that yield a peak have an interior hump descending on both sides. The two that do not are monotone across the whole lag range, rising in one case and falling in the other, so the best alignment lies at or beyond ±4 s. The algorithm is not reporting an absence of coupling; it is reporting that the coupling peaked further away in time than we allowed it to look. The missing windows cluster exactly where a band is running off the top or bottom edge.

Panel B is panel A with everything but the peak discarded — the single number per window that the rest of the analysis is built on.

The same surface for all sixteen sessions

WINDOW_MIN <- 6   # minutes shown per session, from the start of usable data

all_surf <- cached("all_surf", {
  map_dfr(1:16, \(s) {
    d <- filter(signals, session == s)
    M <- wcc::wccCalc(d$patient, d$therapist, wMax = WMAX, tMax = TMAX,
                      wInc = WINC, tInc = 1, method = "cumc")
    P <- wcc::wccPeakPick(M, Lsize = 8, pspan = 0.25, type = "Max")
    tt <- ((WMAX + TMAX) + (seq_len(nrow(M)) - 1) * WINC - WMAX / 2 - 1) / FS
    keep <- which(tt <= min(tt) + WINDOW_MIN * 60)
    expand_grid(i = keep, lag = seq(-TMAX, TMAX) / FS) |>
      mutate(session = s, t = tt[i], r = as.numeric(t(M[keep, , drop = FALSE])),
             peak = rep(P$maxIndex[keep] / 2 / FS, each = 2 * TMAX + 1)) |>
      select(session, t, lag, r, peak)
  })
})

peak_track <- distinct(all_surf, session, t, peak) |> filter(!is.na(peak))

ggplot(all_surf, aes(t / 60, lag, fill = r)) +
  geom_raster(interpolate = TRUE) +
  geom_point(data = peak_track, inherit.aes = FALSE, aes(t / 60, peak),
             colour = "#1b9e77", size = .35) +
  geom_hline(yintercept = 0, colour = "grey20", linewidth = .25, linetype = 2) +
  facet_wrap(~ session, ncol = 2, scales = "free_x",
             labeller = labeller(session = \(s) paste("Session", s))) +
  scale_fill_gradient2(low = "#2c5f7c", mid = "grey97", high = "#a8322d",
                       midpoint = 0, limits = c(-1, 1), name = "r") +
  scale_y_continuous(breaks = seq(-4, 4, 2), expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0)) +
  labs(title = sprintf("The first %g minutes of every session", WINDOW_MIN),
       subtitle = paste("Green dots are picked peaks; blank stripes are windows with none.",
                        "Same colour scale throughout."),
       x = "minutes into session", y = "lag (s)   ↑ therapist leads")

The banded structure is present in all sixteen, with the same characteristic width, which is reassurance that it is a property of the signals rather than of one recording. What differs is how much of each surface is strongly coloured and how continuous the peak track is — sessions 2 and 11 carry long unbroken runs, while 4 and 16 are more fragmented.

Two quantities come out of this step, and both are of interest. The peak correlation says how closely the two series matched at their best alignment. The peak lag says what that alignment was, and its variability over a session says how stable the lead–follow relationship was.

---

Step 4 — Read the peaks at the coded moments

Each episode spans a set of windows, and the measures are averaged over them. Let episode \(e\) cover the \(n_e\) windows \(w_1, \dots, w_{n_e}\) whose centres fall inside it; then

\[\bar{r}_e \;=\; \frac{1}{n_e} \sum_{j=1}^{n_e} \hat{r}_{w_j}\]

and likewise for the lag. Windows with no peak contribute nothing, so an episode’s mean rests on however many of its windows yielded one.

Each episode is compared against a duration-matched control window: same session, same length, drawn at random from a stretch no coder marked, five draws averaged.

MEASURES <- c("peak_r", "peak_lag", "mean_absr")

window_mean <- function(s, from, to, v) {
  x <- peaks[[v]][peaks$session == s & peaks$t >= from & peaks$t < to]
  if (sum(!is.na(x)) < 2) NA_real_ else mean(x, na.rm = TRUE)
}

matched <- cached("matched", {
  set.seed(11)
  map_dfr(seq_len(nrow(episodes)), \(i) {
    e <- episodes[i, ]
    occupied <- filter(episodes, session == e$session)
    t_rng <- range(peaks$t[peaks$session == e$session])
    draws <- map(1:5, \(.) {
      for (k in 1:60) {
        from <- runif(1, t_rng[1], t_rng[2] - e$duration)
        if (!any(from < occupied$end & (from + e$duration) > occupied$start))
          return(map_dbl(MEASURES, \(v) window_mean(e$session, from, from + e$duration, v)))
      }
      rep(NA_real_, length(MEASURES))
    })
    tibble(id = e$id, session = e$session, type = e$type, duration = e$duration,
           measure = MEASURES,
           episode = map_dbl(MEASURES, \(v) window_mean(e$session, e$start, e$end, v)),
           control = colMeans(do.call(rbind, draws), na.rm = TRUE))
  }) |>
    filter(!is.na(episode), !is.na(control)) |>
    mutate(gap = episode - control)
})

matched |>
  summarise(n = n(), episode = round(mean(episode), 4),
            control = round(mean(control), 4),
            difference = round(mean(gap), 4),
            t = round(t.test(gap)$statistic, 2),
            p = round(t.test(gap)$p.value, 3),
            .by = measure)

measure	n	episode	control	difference	t	p
peak_r	158	0.4934	0.5321	-0.0387	-1.40	0.163
peak_lag	158	0.0287	0.0001	0.0285	0.33	0.742
mean_absr	172	0.5161	0.5036	0.0125	1.54	0.125

dists <- matched |>
  filter(measure == "peak_r") |>
  select(episode, control) |>
  pivot_longer(everything(), names_to = "window", values_to = "r") |>
  mutate(window = factor(window, c("episode", "control")))

panel_dist <- ggplot(dists, aes(r, fill = window)) +
  geom_histogram(bins = 40, position = "identity", alpha = .55, colour = NA) +
  scale_fill_manual(values = c(episode = "#c05746", control = "grey55"), name = NULL) +
  labs(title = "A. Peak correlation at coded episodes and at matched controls",
       x = expression(mean~peak~r), y = "episodes")

annot <- matched |>
  summarise(mean = mean(gap), sd = sd(gap), .by = measure) |>
  mutate(txt = sprintf("mean %+.4f\nSD %.3f", mean, sd))

panel_gap <- ggplot(matched, aes(gap)) +
  geom_histogram(bins = 45, fill = "grey78", colour = "white", linewidth = .15) +
  geom_vline(xintercept = 0, colour = "grey20", linewidth = .4) +
  geom_vline(data = annot, aes(xintercept = mean), colour = "#c05746", linewidth = .7) +
  geom_text(data = annot, aes(x = Inf, y = Inf, label = txt), hjust = 1.05,
            vjust = 1.25, size = 3, colour = "grey30", lineheight = .95) +
  facet_wrap(~ measure, ncol = 3, scales = "free") +
  labs(title = "B. Episode minus control, one value per episode",
       subtitle = "Black line is zero; red line is the mean difference.",
       x = "difference", y = "episodes")

panel_dist / panel_gap

Nothing distinguishes the coded moments. Peak correlation is \(0.493\) at episodes against \(0.532\) at controls, \(t = -1.40\), \(p = .16\) — if anything slightly lower during coded episodes. Peak lag is unchanged, \(p = .74\). The distributions in panel A are almost coincident, and panel B shows the difference centred near zero with a spread an order of magnitude larger.

The \(t\)-tests treat 199 episodes as independent when they arrive in clusters of seven to twenty-three from fifteen sessions, so a session random intercept is the minimum correction.

library(lme4)
library(lmerTest)

map_dfr(MEASURES, \(v) {
  fit <- lmer(gap ~ 1 + (1 | session), data = filter(matched, measure == v))
  co  <- summary(fit)$coefficients
  vc  <- as.data.frame(VarCorr(fit))
  tibble(measure = v,
         estimate = round(co[1, "Estimate"], 4),
         SE = round(co[1, "Std. Error"], 4),
         df = round(co[1, "df"], 1),
         t = round(co[1, "t value"], 2),
         p = round(co[1, "Pr(>|t|)"], 3),
         ICC = round(vc$vcov[1] / sum(vc$vcov), 3))
})

measure	estimate	SE	df	t	p	ICC
peak_r	-0.0387	0.0276	157.0	-1.40	0.163	0.000
peak_lag	0.0298	0.1131	15.4	0.26	0.796	0.069
mean_absr	0.0125	0.0081	171.0	1.54	0.125	0.000

It changes nothing. The intraclass correlations are at or near zero, so almost none of the variance in the episode-minus-control difference lies between sessions. The concern was legitimate and the correction was worth making; the answer is that the clustering is negligible. What no random effect can supply is the level above it, since all sixteen sessions come from one dyad.

---

Step 5 — Test against a null

Step 4 asks whether coupling is elevated at particular moments. It does not ask whether these two people are coupled at all, and a mean peak correlation of \(0.53\) cannot answer that on its own — two unrelated skin conductance series will also share variance at their best alignment, and the peak is the maximum of 81 candidates, which raises it mechanically.

A surrogate must destroy the coupling and preserve everything else. With one dyad there is no second dyad to shuffle against, but this pair was recorded sixteen times: pair the patient’s session \(i\) with the therapist’s session \(j \neq i\). Each person keeps their own responsiveness, their own autocorrelation and their own drift. Only the shared timing is destroyed. Note how stringent this is — the surrogate contains the same two people, so anything stable about them as individuals survives into the null, and only coupling specific to being in the room together on that day can show up.

null_a <- cached("null_a", {
  summarise_pair <- function(P) {
    tibble(peak_r = mean(P$peak_r, na.rm = TRUE),
           sd_lag = sd(P$peak_lag, na.rm = TRUE),
           mean_absr = mean(P$mean_absr, na.rm = TRUE))
  }
  set.seed(2)
  mismatches <- map_dfr(1:16, \(s) tibble(p = s, t = sample(setdiff(1:16, s), 3)))
  bind_rows(
    map_dfr(1:16, \(s) summarise_pair(peaks_for(patient[[s]], therapist[[s]]))) |>
      mutate(kind = "real"),
    map2_dfr(mismatches$p, mismatches$t, \(p, t)
      summarise_pair(peaks_for(patient[[p]], therapist[[t]]))) |>
      mutate(kind = "surrogate"))
})

map_dfr(c("peak_r", "sd_lag", "mean_absr"), \(v) {
  r <- null_a[[v]][null_a$kind == "real"]
  m <- null_a[[v]][null_a$kind == "surrogate"]
  test <- t.test(r, m)
  tibble(measure = v,
         `real (n=16)` = round(mean(r), 4), sd = round(sd(r), 4),
         `surrogate (n=48)` = round(mean(m), 4),
         difference = round(mean(r) - mean(m), 4),
         t = round(test$statistic, 2), p = round(test$p.value, 4))
})

measure	real (n=16)	sd	surrogate (n=48)	difference	t	p
peak_r	0.5297	0.0363	0.4983	0.0315	2.89	0.0072
sd_lag	1.7362	0.0535	1.7878	-0.0517	-3.42	0.0022
mean_absr	0.5135	0.0505	0.5104	0.0031	0.23	0.8172

null_a |>
  pivot_longer(c(peak_r, sd_lag, mean_absr), names_to = "measure", values_to = "value") |>
  mutate(kind = factor(kind, c("real", "surrogate"),
                       c("real dyad, same session", "different sessions"))) |>
  ggplot(aes(value, kind, fill = kind)) +
  geom_boxplot(width = .5, outlier.size = .7, linewidth = .3) +
  facet_wrap(~ measure, ncol = 3, scales = "free_x") +
  scale_fill_manual(values = c("real dyad, same session" = "#c05746",
                               "different sessions" = "grey78"), guide = "none") +
  labs(title = "Whole-session measures, real pairings against mismatched ones",
       x = NULL, y = NULL)

Two of the three measures separate the real pairings from the mismatched ones.

Peak correlation is \(0.530\) for real sessions against \(0.498\) for mismatched, \(t = 2.89\), \(p = .007\).

The variability of the peak lag separates them more strongly, and in the direction that matters: \(1.736\) s for real sessions against \(1.788\) s for mismatched, \(t = -3.42\), \(p = .002\). When these two people are actually in the room together, the interval by which one leads the other is more stable than when their sessions are artificially paired. This is the quantity windowed cross-correlation exists to provide — not how strong the association is, but how consistently timed it is — and here it is the better discriminator.

The third measure does not separate them at all. Averaging the absolute correlation over every lag and every window, rather than picking peaks, gives \(0.5135\) against \(0.5104\), \(t = 0.23\), \(p = .82\). This is a widely used alternative aggregation, and on these data it finds nothing where peak-picking finds something. The information is in which alignment is best and how steady that alignment is, and averaging across all alignments destroys it.

A second surrogate, not run here

This test establishes that the dyad is coupled. It does not test the Step 4 question, and the appropriate surrogate for that one is a circular shift of the markers: hold both signals fixed and move all of a session’s episode windows by a common random offset, wrapping at the end. The physiology, the number of episodes and their durations and relative spacing are all preserved; only the alignment between codes and signal is broken. A thousand such shifts give a null distribution for the mean episode value that requires no assumption of independence, because each draw preserves whatever dependence the real episodes have. It is inexpensive, since the peaks are already computed.

---

What the analysis supports

1. The dyad is coupled, by both peak correlation (\(p = .007\)) and the stability of the peak lag (\(p = .002\)), against a surrogate that preserves both individuals and destroys only their shared timing.
2. Coupling is not elevated at the coded moments — peak correlation \(p = .16\), peak lag \(p = .74\), unchanged under a session random intercept.
3. The result depends on peak-picking. The same windows aggregated as a mean absolute correlation show nothing (\(p = .82\)).

The second point is a null and a specific one. It does not say the coded moments are physiologically unremarkable; it says that whatever distinguishes them is not a change in how strongly or how steadily these two autonomic signals track one another.

Would a different choice have changed the answer?

The parameters would, and Step 2 shows by how much: across 56 configurations the separation statistic ranges from 0.3 to 3.8, and a published recommendation for this signal lands near the bottom. That range is wider than most of the effects reported here, which is the strongest argument in the document for running the sweep rather than adopting a number.

The aggregation rule would too, and more sharply: peak-picking and absolute-mean disagree completely on the same matrices. These are not variants of one method but different claims about what synchrony is.

The missing-peak rule is the remaining degree of freedom and it has no settled convention. Here the peaks are simply absent from the averages. At these parameters that discards 21% of windows, and because those windows are mostly ones where the best alignment lay outside \(\pm 4\) s, the surviving estimate is conditioned on the coupling having been reasonably prompt.

Reporting checklist

	This analysis
Signal and sampling rate	Skin conductance, recorded at 100 Hz
Downsampling	signal::decimate, FIR, 30 taps, to 10 Hz
Artifact handling	None applied; see Limits
Window	10 s (100 samples)
Window increment	2 s (20 samples)
Maximum lag	±4 s (±40 samples)
Lag increment	1 sample
Peak picking	wcc::wccPeakPick, Lsize = 8, pspan = 0.25, type = "Max"
Missing peaks	21.1% of windows; excluded from averages
Parameter justification	56-cell sweep, reported in Step 2
Episode aggregation	Mean over windows centred inside the episode
Comparison	Five duration-matched control windows per episode
Episode model	Random intercept for session, lme4::lmer
Surrogate	48 cross-session pairings
Software	R version 4.6.1 (2026-06-24); wcc 0.4.2; signal 1.8.1; lme4 2.0.6

Limits

One dyad. Sixteen sessions are repeated measurements of one relationship, not a sample of relationships. Nothing generalises beyond this pair.

Artifacts were not removed. The deposit includes a log of artifact periods, but its timestamps mix mm:ss and hh:mm:ss formats and one entry reads “until the end”. Reconciling it was out of scope. Brief artifacts create sharp excursions that a 10-second window will register, so the Step 5 effect should be re-checked with the log applied.

A fifth of the windows contribute nothing, and not at random: they are mostly windows whose best alignment fell outside the lag range. A longer lag would recover them at the cost of a weaker separation, as Step 2 shows.

Skin conductance indexes sympathetic activity only.

---

References

Boker, S. M., Xu, M., Rotondo, J. L., & King, K. (2002). Windowed cross-correlation and peak picking for the analysis of variability in the association between behavioral time series. Psychological Methods, 7(3), 338–355.

Behrens, F., Snijdewint, J. A., Moulder, R. G., Prochazkova, E., Sjak-Shie, E. E., Boker, S. M., & Kret, M. E. (2020). Physiological synchrony is associated with cooperative success in real-life interactions. Scientific Reports, 10, 19609.

Behrens, F., Diana, F., Moulder, R. G., Boker, S. M., & Kret, M. E. (2026). Quantifying physiological synchrony using windowed cross-correlation analysis: Statistical and theoretical considerations. bioRxiv.