Windowed Cross-Correlation
Recovering a lead that moves
Overview
A single correlation between two partners’ physiological series asks one question: do they move together at the same instant, in the same way, for the whole recording? Three assumptions are buried in that sentence, and all three are usually wrong. Coupling can be delayed. The delay can change. And the coupling itself can come and go.
Windowed cross-correlation (WCC) relaxes all three. Rather than one coefficient, it computes a correlation for every combination of where in the recording you look and how far you shift one partner against the other. The result is a surface, not a number.
This tutorial starts from two equally spaced time series — the point at which beats have already been detected, screened, and placed on a common grid — and ends at that surface. It then asks what the surface is worth, and what it costs.
All data are simulated, and nothing is read from disk. The delay is therefore known exactly — but it is not shown until Section 5, so that the earlier sections can be read the way a real recording would be. The generating code is in the appendix.
Setup
PERSON <- c(A = "#1f6f8b", B = "#c05746")
# Lead direction. Purple-green, deliberately not the red-blue used for
# correlation values in the surfaces later on.
LEAD_COLS <- c(low = "#4d9221", mid = "grey98", high = "#762a83")
# Correlation values. Used for every r in the document, so that a colour means
# the same number in Figures 3b, 3c and 4.
R_COLS <- c(low = "#2c5f7c", mid = "grey96", high = "#a8322d")
scale_r_fill <- function(...) {
scale_fill_gradient2(low = R_COLS[["low"]], mid = R_COLS[["mid"]],
high = R_COLS[["high"]], midpoint = 0,
limits = c(-1, 1), ...)
}
scale_r_colour <- function(...) {
scale_colour_gradient2(low = R_COLS[["low"]], mid = R_COLS[["mid"]],
high = R_COLS[["high"]], midpoint = 0,
limits = c(-1, 1), ...)
}
# One fixed colour off the lead scale, for shading a single known window.
lead_colour <- function(lag_s) {
scales::div_gradient_pal(LEAD_COLS[["low"]], LEAD_COLS[["mid"]],
LEAD_COLS[["high"]], space = "Lab")(
scales::rescale(lag_s, from = c(-9, 9)))
}
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.5),
strip.text = element_text(face = "bold", size = 9.5),
legend.position = "top"
)
)1. Two dyads
Each partner’s series is 300 observations on a 2-second grid — ten minutes of recording, the resolution a beat-level signal comfortably supports.
These data are simulated, so there is a correct answer, but it is withheld until Section 5. The code that generated them is in the appendix. Until then the two dyads should be treated as any pair of recordings would be: something has to be recovered from them, and nothing about it is known in advance.
A note on the word “leads”
Throughout, A leads B by \(\delta\) is shorthand for a purely temporal statement: A’s trajectory reappears in B’s \(\delta\) seconds later. It is not a claim that A caused anything. Section 5.2 returns to the difference, which matters more than the shorthand suggests.
N <- 300 # observations
EPOCH_S <- 2 # seconds per observation
WINDOW <- 30 # observations per window = 60 s
STEP <- 5 # window increment = 10 s
MAX_LAG <- 10 # lags searched, in observations = +/- 20 sWith \(T = 300\) observations per person, here is everything we have.
dyads |>
pivot_longer(c(A, B), names_to = "person", values_to = "value") |>
ggplot(aes(t, value, colour = person)) +
geom_line(linewidth = 0.45) +
scale_colour_manual(values = PERSON) +
facet_wrap(~ dyad, ncol = 1) +
labs(title = "1. The two dyads",
subtitle = "Ten minutes each, on a 2-second grid. Two labels, and no other information.",
x = "time (s)", y = "value (SD units)", colour = NULL)
The panel titles are the only hint, and they are deliberately vague. Nothing in either picture shows a delay. That is the ordinary situation rather than a property of this simulation: lags in physiological data are a few seconds against recordings of many minutes, and no amount of staring at the raw series will reveal them.
2. What a single correlation reports
The usual analysis computes one coefficient over the whole recording, pairing each observation with the partner’s observation at the same instant:
\[r_0 = \frac{\sum_{t=1}^{T} (A_t - \bar{A})(B_t - \bar{B})} {\sqrt{\sum_{t=1}^{T} (A_t - \bar{A})^2 \; \sum_{t=1}^{T} (B_t - \bar{B})^2}}\]
## # A tibble: 2 × 2
## dyad r_0
## <chr> <dbl>
## 1 Dyad 1 0.22
## 2 Dyad 2 0.449Both dyads are strongly coupled by construction. A single contemporaneous correlation reports 0.22 for the first and 0.45 for the second.
The loss has two sources, and they compound. Because the coupling is delayed, comparing the partners at the same instant compares the wrong pair of moments. And if the delay is not constant, errors made early in the recording need not be in the same direction as errors made late, so averaging across the whole thing cancels part of what survives.
3. Windows and lags
Two indices replace the single coefficient of Section 2. Let \(w = 1, \dots, W\) index windows and \(\ell = -L, \dots, L\) index lags.
A window of length \(m\) advancing in steps of \(\Delta\) starts at observation \(s_w = 1 + (w - 1)\Delta\), so there are \(W = \lfloor (T - m)/\Delta \rfloor + 1\) of them. With \(T = 300\), \(m = 30\) and \(\Delta = 5\) that gives \(W = 55\) windows, and with \(L = 10\) there are \(2L + 1 = 21\) lags. The quantity computed at each combination is
\[r_{w\ell} = \operatorname{cor}\big( A_{\,s_w \,:\, s_w + m - 1},\; B_{\,s_w + \ell \,:\, s_w + m - 1 + \ell}\big)\]
— partner A’s window held fixed, partner B’s shifted by \(\ell\). Note that \(r_{w0}\) is not \(r_0\): the former is computed inside one window, the latter across the whole recording.
The windows overlap heavily. Each advances by \(\Delta = 5\) observations (10 s) while spanning \(m = 30\) (60 s), so consecutive windows share 83% of their data.
cross_correlate <- function(a, b, start, window = WINDOW, max_lag = MAX_LAG) {
tibble(lag = -max_lag:max_lag) |>
mutate(r = map_dbl(lag, \(l) {
ia <- start:(start + window - 1)
ib <- ia + l
ok <- ib >= 1 & ib <= length(b)
if (sum(ok) < 10) NA_real_ else cor(a[ia[ok]], b[ib[ok]])
}))
}
d1 <- filter(dyads, str_starts(dyad, "Dyad 1"))
W_TOTAL <- (N - WINDOW) %/% STEP + 1
six <- tibble(w = c(1:3, (W_TOTAL - 2):W_TOTAL)) |>
mutate(
first_obs = 1 + (w - 1) * STEP,
last_obs = first_obs + WINDOW - 1,
t_start = first_obs * EPOCH_S,
t_end = last_obs * EPOCH_S,
half = if_else(w <= 3, "first three windows", "last three windows"),
slot = factor(rep(c("1st", "2nd", "3rd"), 2),
levels = c("1st", "2nd", "3rd")),
cc = map(first_obs, \(s) cross_correlate(d1$A, d1$B, start = s)),
r_w0 = map_dbl(cc, \(x) x$r[x$lag == 0]),
peak_lag = map_dbl(cc, \(x) x |> drop_na(r) |> slice_max(r, n = 1) |> pull(lag)),
peak_r = map_dbl(cc, \(x) max(x$r, na.rm = TRUE))
)d1_long <- d1 |> pivot_longer(c(A, B), names_to = "person", values_to = "value")
panels <- six |>
mutate(lo = if_else(half == "first three windows", 1, N - 49),
hi = if_else(half == "first three windows", 50, N),
series = map2(lo, hi, \(a, b) filter(d1_long, i >= a, i <= b))) |>
unnest(series) |>
mutate(inside = i >= first_obs & i <= last_obs)
ggplot(panels, aes(t, value)) +
geom_rect(data = six, inherit.aes = FALSE,
aes(xmin = t_start, xmax = t_end, ymin = -Inf, ymax = Inf,
fill = peak_lag * EPOCH_S), alpha = 0.32) +
geom_line(aes(group = person), colour = "grey72", linewidth = 0.4) +
geom_line(data = \(d) filter(d, inside), aes(colour = person), linewidth = 0.6) +
geom_text(data = six, inherit.aes = FALSE, size = 2.9, colour = "grey10",
aes(x = (t_start + t_end) / 2, y = 4.3,
label = sprintf("w = %d r at lag 0 = %+.2f peak at %+g s",
w, r_w0, peak_lag * EPOCH_S))) +
scale_fill_gradient2(low = LEAD_COLS[["low"]], mid = LEAD_COLS[["mid"]],
high = LEAD_COLS[["high"]], midpoint = 0,
limits = c(-9, 9), breaks = c(-8, -4, 0, 4, 8),
labels = c("B precedes by 8", "4", "0", "4", "A precedes by 8")) +
scale_colour_manual(values = PERSON) +
coord_cartesian(ylim = c(-3.3, 4.8)) +
guides(fill = guide_colourbar(barwidth = 11, barheight = 0.5, order = 1),
colour = guide_legend(order = 2)) +
facet_grid(slot ~ half, scales = "free_x") +
labs(title = "3a. The first three windows and the last three",
subtitle = paste0("Dyad 1. Shading marks the window and is coloured by the lag at which its ",
"correlation peaks."),
x = "time (s)", y = "value (SD units)", colour = NULL, fill = NULL)
Two things in those annotations are worth pausing on. The within-window contemporaneous correlation \(r_{w0}\) ranges from -0.33 to +0.15 — near zero, and negative in four of the six. Yet the peak correlation in the same windows runs from 0.57 to 0.78. These partners are strongly coupled and a contemporaneous measure reports nothing, or worse than nothing.
And the shading flips colour between the columns. Every one of the first three windows peaks at +6 s and every one of the last three at -6 s. Whatever relationship holds at the start of this recording has reversed by the end.
Within any one window the building block is an ordinary cross-correlation: hold a stretch of partner A fixed, slide partner B across it, and correlate at each displacement.
The figure below takes one window and does that explicitly. It uses \(w = 14\) (observations 66–95, 132–190 s), chosen because it is far enough from the start of the recording that every lag has all 30 pairs available, and because Figure 3a suggests the relationship is stable across it.
Each row is one lag. The blue trace is partner A’s window and never moves. The grey dashed trace is partner B where it actually sits, also the same in every row. The red trace is B’s values drawn back at A’s time points — the pairing the correlation actually evaluates.
TIER_W <- 14L
TIER_LAGS <- c(-6L, -3L, 0L, 3L, 6L)
s_tier <- 1 + (TIER_W - 1) * STEP
idx <- s_tier:(s_tier + WINDOW - 1)
tier <- expand_grid(lag = TIER_LAGS, i = idx) |>
mutate(t = i * EPOCH_S, A = d1$A[i], B_rest = d1$B[i], B_shifted = d1$B[i + lag]) |>
mutate(r = cor(A, B_shifted), .by = lag) |>
mutate(panel = fct_reorder(sprintf("lag = %+g s", lag * EPOCH_S), lag))
tags <- distinct(tier, lag, panel, r)
# The lag-0 row is the reference: shaded as a window, as in Figure 3a, and
# boxed. Its shading colour is this window's own peak lag.
marked <- tags |>
filter(lag == 0L) |>
mutate(shade = lead_colour(tags$lag[which.max(tags$r)] * EPOCH_S))
ggplot(tier, aes(t)) +
geom_rect(data = marked, inherit.aes = FALSE,
aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf),
fill = marked$shade, alpha = 0.32) +
geom_rect(data = marked, inherit.aes = FALSE,
aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf),
fill = NA, colour = "grey10", linewidth = 1) +
geom_line(data = \(d) filter(d, lag != 0L), aes(y = B_rest),
colour = PERSON[["B"]], linewidth = 0.5, alpha = 0.38) +
geom_line(aes(y = A), colour = PERSON[["A"]], linewidth = 0.6) +
geom_line(data = \(d) filter(d, lag != 0L), aes(y = B_shifted),
colour = PERSON[["B"]], linewidth = 0.65, linetype = "31") +
geom_line(data = \(d) filter(d, lag == 0L), aes(y = B_shifted),
colour = PERSON[["B"]], linewidth = 0.65) +
geom_label(data = tags, inherit.aes = FALSE, size = 2.9, label.size = 0.25,
label.padding = unit(0.16, "lines"), colour = "grey10",
aes(x = mean(range(tier$t)), y = 3.9, fill = r,
label = sprintf("r = %+.3f", r))) +
scale_r_fill(guide = "none") +
coord_cartesian(ylim = c(-2.3, 4.4)) +
facet_wrap(~ panel, ncol = 1) +
labs(title = "3b. The same window, with B slid across it",
subtitle = paste0("w = ", TIER_W, ". Blue is A, and never moves.\n",
"The boxed row is the window itself, shaded\n",
"by its peak lag as in Figure 3a. Elsewhere\n",
"pale red is B at\n",
"rest, dashed red is B shifted by the lag.\n",
"Each r is boxed in the surface colours."),
x = "time (s)", y = "value (SD units)")
Reading down the tier, the red trace walks across the blue one. At \(-6\) s the two are close to mirror images — B’s peaks landing on A’s troughs — and the correlation is -0.41. At lag zero the red trace lies exactly on the grey dashed one, because nothing has been shifted, and it tracks A hardly at all: -0.10. At \(+6\) s the two lock together and the correlation reaches +0.85. By \(+12\) s they have begun to come apart again.
That single row is what a windowed cross-correlation is for. The same thirty pairs of observations yield anything from -0.41 to +0.85 depending only on how the two series are aligned before correlating.
Doing this at all 21 lags rather than five traces out the whole curve.
ccf_two <- bind_rows(
cross_correlate(d1$A, d1$B, start = s_tier) |>
mutate(label = paste0("w = ", TIER_W, " (A leads)")),
cross_correlate(d1$A, d1$B, start = 1 + (W_TOTAL - 1) * STEP) |>
mutate(label = paste0("w = ", W_TOTAL, " (B leads)"))
)
segments <- ccf_two |>
arrange(label, lag) |>
mutate(lag_to = lead(lag), r_to = lead(r), .by = label) |>
drop_na(lag_to) |>
mutate(r_mid = (r + r_to) / 2)
ggplot(ccf_two, aes(r, lag * EPOCH_S)) +
geom_hline(yintercept = 0, colour = "grey45", linewidth = 0.3, linetype = "dotted") +
geom_vline(xintercept = 0, colour = "grey85", linewidth = 0.3) +
geom_segment(data = segments, inherit.aes = FALSE, linewidth = 1.1,
aes(x = r, xend = r_to, y = lag * EPOCH_S, yend = lag_to * EPOCH_S,
colour = r_mid)) +
geom_point(aes(colour = r), size = 2.2) +
geom_point(data = \(d) slice_max(d, r, n = 1, by = label),
shape = 21, size = 4.6, stroke = 0.9, colour = "grey10", fill = NA) +
ggrepel::geom_text_repel(
data = \(d) filter(d, str_starts(label, paste0("w = ", TIER_W)),
lag %in% TIER_LAGS),
aes(label = sprintf("r_%+g = %+.3f", lag * EPOCH_S, r)),
size = 2.9, colour = "grey15", seed = 7,
box.padding = 0.8, point.padding = 0.4, min.segment.length = 0.2,
force = 6, max.overlaps = Inf, direction = "y", nudge_x = -0.55,
segment.colour = "grey60", segment.linewidth = 0.3
) +
scale_r_colour(guide = "none") +
scale_x_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1),
expand = expansion(mult = 0.02)) +
scale_y_continuous(breaks = seq(-20, 20, 10)) +
facet_wrap(~ label) +
labs(title = "3c. The whole cross-correlation function, in two windows",
subtitle = paste0("Axes swapped so that lag runs vertically, as it does in the surface. ",
"Colour is the correlation,\non the same scale. The ring marks the peak; ",
"the five labelled lags are those of Figure 3b."),
x = "correlation", y = "lag (s)")
The left panel is Figure 3b completed — the five labelled points are its five rows, and the other sixteen lags fill in the rest of the curve. The right panel is the last window of the recording, where the precedence runs the other way and the whole curve is mirrored.
Reading across either panel at \(\ell = 0\) gives the value a contemporaneous analysis would report, and in both it sits far from the curve’s extreme rather than at it.
The axes are swapped deliberately. With lag running vertically and correlation carried by colour, each panel is one vertical slice through the surface of the next section — a single window’s column, with the red bulge marking where that window’s coupling actually lives.
4. All windows, all lags
Repeating that across overlapping windows gives the surface.
wcc_surface <- function(a, b, window = WINDOW, step = STEP, max_lag = MAX_LAG) {
seq(1, length(a) - window + 1, by = step) |>
set_names() |>
map(\(s) cross_correlate(a, b, s, window, max_lag) |> mutate(start = s)) |>
list_rbind() |>
mutate(centre = (start + window / 2) * EPOCH_S)
}
surfaces <- dyads |>
group_split(dyad) |>
map(\(d) wcc_surface(d$A, d$B) |> mutate(dyad = d$dyad[1])) |>
list_rbind()
surfaces |>
summarise(windows = n_distinct(start), lags = n_distinct(lag),
cells = n(), .by = dyad)## # A tibble: 2 × 4
## dyad windows lags cells
## <chr> <int> <int> <int>
## 1 Dyad 1 55 21 1155
## 2 Dyad 2 55 21 1155surfaces |>
ggplot(aes(centre, lag * EPOCH_S, fill = r)) +
geom_raster(interpolate = TRUE) +
geom_hline(yintercept = 0, colour = "grey30", linewidth = 0.3, linetype = "dotted") +
scale_r_fill() +
facet_wrap(~ dyad, ncol = 1) +
labs(title = "4. The windowed cross-correlation surface",
subtitle = "Red is positive correlation. The dotted line marks lag 0.",
x = "window centre (s)", y = "lag (s)", fill = "r")
The surface replaces one number with 1155 of them per dyad, and the structure a single coefficient reported as 0.22 is now plainly visible. In dyad 1 a band of high correlation sits above zero for the first half of the recording and below it for the second. In dyad 2 the band descends steadily from positive to negative.
Neither pattern was visible in Figure 1, and neither was recoverable from a single coefficient.
5. Reading the lead off the surface
Taking the highest correlation in each window gives an estimate for each moment in the recording:
\[\hat{\ell}_w = \operatorname*{arg\,max}_{\ell \in \{-L, \dots, L\}} r_{w\ell}, \qquad \hat{r}_w = \max_{\ell} r_{w\ell}\]
peaks <- surfaces |>
drop_na(r) |>
slice_max(r, n = 1, by = c(dyad, start)) |>
mutate(mid = floor(start + WINDOW / 2)) |>
left_join(select(dyads, dyad, i, true_lead), by = c("dyad", "mid" = "i"))5.1 Peak lag, and its variability
Two numbers summarise the sequence of \(\hat{\ell}_w\), and both are needed.
peaks |>
summarise(
`mean peak lag (s)` = round(mean(lag) * EPOCH_S, 2),
`SD of peak lag (s)` = round(sd(lag) * EPOCH_S, 2),
`mean peak r` = round(mean(r), 3),
.by = dyad
)## # A tibble: 2 × 4
## dyad `mean peak lag (s)` `SD of peak lag (s)` `mean peak r`
## <chr> <dbl> <dbl> <dbl>
## 1 Dyad 1 -0.69 6.57 0.69
## 2 Dyad 2 -0.07 4.29 0.735Reported alone, the mean would say there is nothing here. It is -0.69 s for dyad 1 and -0.07 s for dyad 2 — in both cases indistinguishable from no lead at all, because each dyad spends as much of the recording with one partner in front as with the other, and the two halves cancel.
The standard deviation is where the finding is. At 6.57 s and 4.29 s it says the precedence is large and moves. A dyad in which one partner is consistently ahead would show the opposite pattern — a mean far from zero and a small SD — and the two cases are not distinguishable from the mean alone.
This is why a windowed analysis should report the peak lag and its variability. Section 7.1 shows the variability doing real work: across five dyads it is the single best predictor of a dyad-level covariate, better than any correlation-based summary.
5.2 Now the answer
Only at this point is it worth looking at what actually generated the data.
peaks |>
ggplot(aes(centre)) +
geom_hline(yintercept = 0, colour = "grey85", linewidth = 0.3) +
geom_step(aes(y = true_lead * EPOCH_S, colour = "truth"), linewidth = 0.6) +
geom_point(aes(y = lag * EPOCH_S, colour = "recovered"), size = 1.8) +
scale_colour_manual(values = c("truth" = "grey45", "recovered" = PERSON[["A"]])) +
facet_wrap(~ dyad) +
labs(title = "5. Recovered lead against the truth",
subtitle = "One point per window. The grey step is the simulation's own parameter.",
x = "window centre (s)", y = "A precedes B by (s)", colour = NULL)
peaks |>
summarise(
`% lag exactly right` = round(100 * mean(lag == true_lead), 1),
`% within one step` = round(100 * mean(abs(lag - true_lead) <= 1), 1),
`true SD of lead (s)` = round(sd(true_lead) * EPOCH_S, 2),
`recovered SD (s)` = round(sd(lag) * EPOCH_S, 2),
.by = dyad
)## # A tibble: 2 × 5
## dyad `% lag exactly right` `% within one step` `true SD of lead (s)`
## <chr> <dbl> <dbl> <dbl>
## 1 Dyad 1 92.7 94.5 6.05
## 2 Dyad 2 89.1 100 4.29
## # ℹ 1 more variable: `recovered SD (s)` <dbl>Dyad 1 was built with A ahead by 6 s for the first five minutes and B ahead by 6 s for the second; dyad 2’s lead slides smoothly from \(+8\) s to \(-8\) s. Both are recovered, and the variability reported in Section 5.1 before the answer was known lands close to the truth in each case.
Where the estimate is worst is also informative: the recovered lead scatters most where the true lead is changing fastest, because a window spanning a change contains two different relationships and the peak must pick one.
5.3 What “A precedes B” does not establish
The data were generated from a single shared process reaching B later than A:
\[A_t = x_t + \varepsilon^A_t, \qquad B_t = x_{t - \delta_t} + \varepsilon^B_t\]
with \(\varepsilon^A\) and \(\varepsilon^B\) independent of each other and of \(x\). There is no term carrying A into B. Both series are effects of a common cause, sampled at different delays, and the arrows run \(x \to A\) and \(x \to B\) — never \(A \to B\).
That matters, because the surface looks exactly the same as it would under genuine transmission. At least three structures produce a peak at \(+\delta\):
| Structure | What is happening |
|---|---|
| Common cause, differential delay | both respond to the same thing; B’s response is slower. This is the simulation. |
| Transmission | \(B_t = \alpha A_{t-\delta} + \varepsilon\). A really does drive B. |
| Differential measurement latency | the same event, recorded through pathways or sensors with different delays. No interpersonal process at all. |
A cross-correlation cannot separate these, and neither can Granger causality or transfer entropy in any practical sense. There is a signature in principle — under transmission, A’s idiosyncratic noise reappears in B but B’s never reappears in A — but detecting it needs far more than 30 observations per window.
So the warranted claim from a peak at \(+\delta\) is temporal precedence, and only that. Separating precedence from influence is a matter of design rather than estimator: a condition in which the partners cannot perceive each other, a perturbation with recovery, or a model that represents the shared driver explicitly.
6. Choosing the window and the lag range
Two parameters dominate, and they pull against each other.
The hard constraint. A window of W
observations compared at a lag of L has only
W - L observations in common. When L
approaches W the overlap vanishes and the correlation is
undefined.
The soft constraint. Well before that point the estimates at extreme lags rest on far fewer observations than those near zero, so the edges of the surface are noisier than the middle — and the peak-picking step searches the noisy edges on equal terms with the reliable centre.
The opposing constraint. Widening the window buys stability but costs time resolution. A window long enough to be stable may be long enough to span a change in the lead, in which case it averages two relationships and reports neither.
recovery <- function(d, window, max_lag) {
if (window - max_lag < 10) {
return(tibble(`% within one step` = NA_real_, `mean peak r` = NA_real_))
}
pk <- wcc_surface(d$A, d$B, window = window, step = STEP, max_lag = max_lag) |>
drop_na(r) |>
slice_max(r, n = 1, by = start) |>
mutate(mid = floor(start + window / 2)) |>
left_join(select(d, i, true_lead), by = c("mid" = "i")) |>
drop_na(true_lead)
tibble(`% within one step` = 100 * mean(abs(pk$lag - pk$true_lead) <= 1),
`mean peak r` = mean(pk$r))
}
expand_grid(window = c(20, 30, 60, 120), max_lag = c(5, 10, 20)) |>
mutate(`overlap at max lag` = window - max_lag) |>
mutate(map2(window, max_lag, \(w, l) recovery(d1, w, l)) |> list_rbind()) |>
mutate(`% within one step` = round(`% within one step`, 1),
`mean peak r` = round(`mean peak r`, 3))## # A tibble: 12 × 5
## window max_lag `overlap at max lag` `% within one step` `mean peak r`
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 20 5 15 93 0.673
## 2 20 10 10 87.7 0.679
## 3 20 20 0 NA NA
## 4 30 5 25 96.4 0.69
## 5 30 10 20 94.5 0.69
## 6 30 20 10 85.5 0.704
## 7 60 5 55 91.8 0.697
## 8 60 10 50 91.8 0.697
## 9 60 20 40 91.8 0.7
## 10 120 5 115 94.6 0.649
## 11 120 10 110 94.6 0.649
## 12 120 20 100 94.6 0.649All three constraints appear in that table. A 20-observation window searched to a lag of 20 has no overlap at all and returns nothing. Holding the window at 30, widening the lag range from 5 to 20 shrinks the overlap from 25 to 10 and drops recovery from 96.4% to 85.5% — the extra lags searched are estimated from too little data. And the 120-observation window, which is the most stable of the four, has the lowest mean peak correlation, because at four minutes it straddles the midpoint switch and dilutes both regimes.
A workable rule: set the maximum lag from the delay the theory predicts, not from the data, and then make the window several times longer than that maximum — but no longer than the interval over which the coupling can be assumed constant. If those two requirements conflict, the design cannot answer the question, and that is worth discovering before the analysis rather than after.
Applied to the surfaces above: autonomic coupling between interacting partners is expected to operate over a few seconds rather than a few minutes, so the lag range was fixed at ±20 s — generous cover for a delay of that order — before any surface was computed, rather than widened until something appeared. The 60-second window is three times that maximum, so even at the extreme lags 20 of its 30 observations still overlap, and it remains well short of the five-minute stretches over which each dyad’s lead is constant.
7. Reducing the surface to a number
Many published WCC analyses do not report a surface. This may be because the driving research question is about between-dyad differences.
A model of the form synchrony ~ relationship quality needs one row per dyad. A surface of 1155 cells is not a row. Reduction is therefore not laziness; it is what makes the second-stage analysis possible at all. The question is which number, because the rules disagree.
Each summary rule corresponds to a different feature of the surface. Taking dyad 1, the three most common are marked below.
s1 <- surfaces |> filter(str_starts(dyad, "Dyad 1")) |> drop_na(r)
p1 <- s1 |> slice_max(r, n = 1, by = start) # each window's peak
g1 <- s1 |> slice_max(r, n = 1) # the single largest cell
ggplot(s1, aes(centre, lag * EPOCH_S, fill = r)) +
geom_raster(interpolate = TRUE) +
# the lag-0 row: what a contemporaneous correlation looks at
annotate("rect", xmin = min(s1$centre) - 5, xmax = max(s1$centre) + 5,
ymin = -1, ymax = 1, fill = NA, colour = "grey10", linewidth = 0.6) +
# each window's peak
geom_point(data = p1, aes(centre, lag * EPOCH_S), inherit.aes = FALSE,
colour = "grey10", size = 1.1) +
# the global maximum
geom_point(data = g1, aes(centre, lag * EPOCH_S), inherit.aes = FALSE,
shape = 21, size = 6, stroke = 1.1, colour = "grey10", fill = NA) +
annotate("segment", x = 152, xend = 152, y = 16, yend = 10,
colour = "grey10", linewidth = 0.4,
arrow = arrow(length = unit(0.1, "cm"))) +
annotate("label", x = 152, y = 19, label = "maximum r anywhere",
size = 3, label.size = 0, fill = "white") +
annotate("segment", x = 600, xend = 578, y = 0, yend = 0,
colour = "grey10", linewidth = 0.4,
arrow = arrow(length = unit(0.1, "cm"))) +
annotate("text", x = 608, y = 0, label = "the lag-0 row", hjust = 0, size = 3) +
annotate("segment", x = 600, xend = 572, y = -11, yend = -6.5,
colour = "grey10", linewidth = 0.4,
arrow = arrow(length = unit(0.1, "cm"))) +
annotate("text", x = 608, y = -12, label = "each window's peak",
hjust = 0, size = 3) +
scale_r_fill() +
scale_x_continuous(limits = c(0, 790), breaks = seq(0, 600, 200)) +
labs(title = "7. The same surface, with the summarised features marked",
subtitle = "Dyad 1. Each summary rule below reads a different part of this picture.",
x = "window centre (s)", y = "lag (s)", fill = "r")
Using the indices from Section 3 — \(w = 1, \dots, W\) over windows and \(\ell = -L, \dots, L\) over lags, with \(W = 55\) and \(L = 10\) — the rules differ only in which cells of \(r_{w\ell}\) they touch and how they combine them.
tibble(
`summary rule` = c("r at lag 0, whole series",
"mean r along the lag-0 row",
"mean r across the whole surface",
"mean |r| across the whole surface",
"mean of each window's peak r",
"maximum r anywhere on the surface"),
# \lvert and \rvert rather than a literal "|": kable escapes pipes inside
# table cells, which would break the maths. Bars here mean absolute value
# and nothing else -- window and lag counts are the scalars W and L.
`how it is computed` = c(
"$r_0$",
"$\\frac{1}{W}\\sum_{w=1}^{W} r_{w0}$",
"$\\frac{1}{W(2L+1)}\\sum_{w=1}^{W}\\sum_{\\ell=-L}^{L} r_{w\\ell}$",
"$\\frac{1}{W(2L+1)}\\sum_{w=1}^{W}\\sum_{\\ell=-L}^{L} \\lvert r_{w\\ell}\\rvert$",
"$\\frac{1}{W}\\sum_{w=1}^{W} \\max_{\\ell} r_{w\\ell}$",
"$\\max_{w,\\,\\ell}\\ r_{w\\ell}$"
),
value = round(c(cor(d1$A, d1$B), mean(s1$r[s1$lag == 0]), mean(s1$r),
mean(abs(s1$r)), mean(p1$r), max(s1$r)), 3)
) |>
knitr::kable(escape = FALSE, align = "llr")| summary rule | how it is computed | value |
|---|---|---|
| r at lag 0, whole series | \(r_0\) | 0.220 |
| mean r along the lag-0 row | \(\frac{1}{W}\sum_{w=1}^{W} r_{w0}\) | 0.028 |
| mean r across the whole surface | \(\frac{1}{W(2L+1)}\sum_{w=1}^{W}\sum_{\ell=-L}^{L} r_{w\ell}\) | 0.074 |
| mean |r| across the whole surface | \(\frac{1}{W(2L+1)}\sum_{w=1}^{W}\sum_{\ell=-L}^{L} \lvert r_{w\ell}\rvert\) | 0.262 |
| mean of each window’s peak r | \(\frac{1}{W}\sum_{w=1}^{W} \max_{\ell} r_{w\ell}\) | 0.690 |
| maximum r anywhere on the surface | \(\max_{w,\,\ell}\ r_{w\ell}\) | 0.858 |
The first row is the odd one out, and the column makes it visible: it is the only expression with no \(w\) or \(\ell\) in it. It is not a summary of the surface at all — it is \(r_0\) from Section 2, computed on the whole series before any window, step or lag range had been chosen, and it appears here only as the baseline the other five are compared against.
One dyad, one surface, five defensible summary rules, and a spread from 0.03 to 0.86.
The summary rules are not interchangeable, because they read different parts of the picture. The lag-0 row — the black band in the figure — is what a contemporaneous correlation attends to, and almost all of the structure in this surface lies outside it. Note also that the row averages to 0.028 while the whole-series correlation is 0.22: correlations computed within windows and then averaged are not the same quantity as one correlation computed across the whole recording, and the two can differ by a good deal.
The mean over the surface is near zero only because the surface holds roughly equal positive and negative regions, an artifact of this dyad’s reversal rather than evidence of weak coupling. The mean absolute value fixes that and destroys the direction. The maximum — the circled cell — is the largest of 1155 estimates and is therefore an extreme-value statistic, biased upward by the number of cells searched.
The mean of the window peaks — the dotted trace along the ridge — is usually the most defensible, and it still requires reporting the window width, the step, and the lag range that produced it.
Note that all six rows summarise correlations. The surface also supports the two lag summaries of Section 5.1 — mean peak lag -0.69 s and its SD 6.57 s for this dyad — which have no counterpart in a single coefficient at all. A complete report of one window’s worth of analysis therefore needs one number from each family.
7.1 Why the reduction earns its keep
Five dyads, differing on a dyad-level covariate — call it closeness. Closer dyads are simulated to be both more strongly coupled (less private noise) and steadier in precedence (a smaller swing in the lead). Nothing else varies.
set.seed(2026)
five <- 1:5 |>
map(\(closeness) {
simulate_lagged(
N,
lead = as.integer(round((6 - closeness) * sin(2 * pi * seq_len(N) / N))),
noise = 0.85 - 0.13 * (closeness - 1)
) |>
mutate(dyad = paste("Dyad", closeness), closeness = closeness)
}) |>
list_rbind()
surfaces_5 <- five |>
group_split(dyad) |>
map(\(d) wcc_surface(d$A, d$B) |> mutate(dyad = d$dyad[1], closeness = d$closeness[1])) |>
list_rbind()
metrics_5 <- surfaces_5 |>
drop_na(r) |>
slice_max(r, n = 1, by = c(dyad, start)) |>
summarise(peak_r = mean(r), mean_lag_s = mean(lag) * EPOCH_S,
lag_sd_s = sd(lag) * EPOCH_S, .by = c(dyad, closeness)) |>
left_join(summarise(five, r_lag0 = cor(A, B), .by = dyad), by = "dyad")left <- surfaces_5 |>
ggplot(aes(centre, lag * EPOCH_S, fill = r)) +
geom_raster(interpolate = TRUE) +
geom_hline(yintercept = 0, colour = "grey30", linewidth = 0.25, linetype = "dotted") +
scale_r_fill() +
facet_wrap(~ dyad, ncol = 1) +
labs(title = "7.1 Five dyads, ordered by closeness",
x = "window centre (s)", y = "lag (s)", fill = "r") +
theme(legend.position = "bottom")
right <- metrics_5 |>
pivot_longer(c(closeness, r_lag0, peak_r, mean_lag_s, lag_sd_s)) |>
mutate(
dyad = fct_rev(factor(dyad)),
name = factor(name,
levels = c("closeness", "r_lag0", "peak_r", "mean_lag_s", "lag_sd_s"),
labels = c("closeness", "r at\nlag 0", "mean\npeak r",
"mean\nlag (s)", "SD of\nlag (s)")),
label = if_else(name == "closeness", sprintf("%.0f", value), sprintf("%.2f", value))
) |>
ggplot(aes(name, dyad, label = label)) +
geom_text(size = 3.3, colour = "grey15") +
scale_x_discrete(position = "top") +
labs(title = "Summary metrics", x = NULL, y = NULL) +
theme(panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(face = "bold", size = 8.5))
left + right + plot_layout(widths = c(2, 1.15))
metrics_5 |>
summarise(across(c(r_lag0, peak_r, mean_lag_s, lag_sd_s),
\(x) round(cor(x, closeness), 3))) |>
pivot_longer(everything(), names_to = "metric", values_to = "r with closeness")## # A tibble: 4 × 2
## metric `r with closeness`
## <chr> <dbl>
## 1 r_lag0 0.844
## 2 peak_r 0.979
## 3 mean_lag_s 0.489
## 4 lag_sd_s -0.984Three things are visible in that comparison, and each matters for the second-stage model.
The surface-derived metrics track the covariate far better than the contemporaneous correlation does. Mean peak r correlates 0.98 with closeness and the SD of the peak lag correlates -0.98, against 0.84 for r at lag zero — which is also non-monotone here, ranking dyad 4 below dyad 3. The single coefficient is not merely a weaker measure of the same thing; it reorders the dyads.
The surface yields a metric with no single-correlation counterpart. The variability of the peak lag indexes how stable the precedence is, and in this simulation it is the best predictor of all. A design that reduced each dyad to one correlation could not have produced this column.
The right summary still has to match the question. Mean lag correlates only 0.49 with closeness, correctly — the simulation varied the swing of the lead, not its average, so the mean is near zero in every dyad. Reporting mean lag here would return a null for an effect that is plainly present.
broom::tidy(lm(peak_r ~ closeness, data = metrics_5)) |>
mutate(across(where(is.numeric), \(x) round(x, 4)))## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.398 0.0355 11.2 0.0015
## 2 closeness 0.09 0.0107 8.41 0.0035With five dyads this regression is an illustration of structure rather than a result; the standard errors are not to be taken seriously. The point is the shape of the data it consumes — one row per dyad, with the within-dyad temporal structure already compressed into columns — and that the compression chosen determines whether the effect survives into it.
A note on the software
Everything above is computed from base cor() in two
short functions — cross_correlate() over lags,
wcc_surface() over windows — because the point was to show
the mechanism rather than call into a routine that hides it. For
analysis rather than teaching, use a maintained implementation.
The method is Boker, Rotondo, Xu & King (2002), Windowed cross-correlation and peak picking for the analysis of variability in the association between behavioral time series, Psychological Methods 7(3), 338 — including the peak-picking step of Section 5, which is theirs.
| Package | Notes |
|---|---|
wcc |
The reference implementation, maintained by Boker.
wccCalc() and wccCalcBatch() build the
surface, wccPeakPick() does Section 5,
wccAggregate() does Section 7,
wccFindDyadParam() searches the parameter grid of Section
6, and wccSurrogateDyads() /
wccTestSurrogateDyads() give the surrogate test this
tutorial does not cover. |
rMEA |
Kleinbub & Ramseyer (2020). Built for motion-energy data but general in its windowed cross-correlation and surrogate routines. |
⚠️ Name collision. The function defined here is
called wcc_surface() rather than wcc()
precisely because wcc is the package name. Calling a local
function wcc() masks the package and produces errors that
are hard to read.
What the packages add beyond this tutorial is mostly inference: surrogate construction and testing, which Section 7 gestures at and never performs.
Checklist
Each item was demonstrated above; the section is given so the claim can be checked rather than taken on trust.
| Do this | Shown in | The evidence | |
|---|---|---|---|
| 1 | Plot the series first, and do not read a flat-looking picture as absence of coupling | §1, §4 | Figure 1 shows nothing; the same data in Figure 4 show a band that reverses |
| 2 | Set the maximum lag from theory, not by widening until something appears | §6 | ±20 s fixed in advance from the few-seconds timescale of autonomic coupling |
| 3 | Make the window several times the maximum lag | §6 | holding \(m = 30\) and widening \(L\) from 5 to 20 shrinks overlap 25→10 and drops recovery 96.4% → 85.5% |
| 4 | But keep the window shorter than the interval over which coupling is constant | §6 | the \(m = 120\) window is the most stable and has the lowest peak \(r\), because it straddles the reversal |
| 5 | Report window, step and lag range with every coefficient | §3, §6 | none is recoverable from the number; the §6 table shows four of twelve combinations returning nothing at all |
| 6 | State the summary rule | §7 | one surface, six rules, 0.03 to 0.86 |
| 7 | Keep the sign | §7 | mean \(r\) over the surface is 0.07; mean \(\lvert r\rvert\) is 0.26. The first is small because the dyad reverses, which is the finding; the second hides it |
| 8 | Report the peak lag and its variability, not only the peak correlation | §5.1, §7.1 | both dyads have a mean peak lag near zero and an SD above 4 s; across five dyads the SD is the best predictor of a dyad-level covariate |
| 9 | Claim precedence, not influence | §5.3 | the data were generated by a common cause with differential delay, and the surface is indistinguishable from genuine transmission |
Appendix: how the data were made
Withheld until now so that Sections 1–4 could be read the way a real recording would be.
ar_series <- function(n, phi = 0.8) {
as.numeric(scale(stats::filter(rnorm(n), phi, method = "recursive")))
}
# `lead` is in observations. lead > 0 means B reproduces the shared process
# `lead` steps AFTER A does -- A's trajectory reappears in B's, later.
simulate_lagged <- function(n, lead, phi = 0.8, noise = 0.55) {
pad <- max(abs(lead)) + 1
x <- ar_series(n + 2 * pad, phi)
idx <- seq_len(n) + pad
tibble(
i = seq_len(n),
t = i * EPOCH_S,
true_lead = lead,
A = x[idx] + noise * rnorm(n),
B = x[idx - lead] + noise * rnorm(n)
)
}
set.seed(814)
dyads <- bind_rows(
simulate_lagged(N, lead = rep(c(3L, -3L), each = N / 2)) |>
mutate(dyad = "Dyad 1"),
simulate_lagged(N, lead = as.integer(round(seq(4, -4, length.out = N)))) |>
mutate(dyad = "Dyad 2")
)set.seed(2026)
five <- 1:5 |>
map(\(closeness) {
simulate_lagged(
N,
lead = as.integer(round((6 - closeness) * sin(2 * pi * seq_len(N) / N))),
noise = 0.85 - 0.13 * (closeness - 1)
) |>
mutate(dyad = paste("Dyad", closeness), closeness = closeness)
}) |>
list_rbind()
surfaces_5 <- five |>
group_split(dyad) |>
map(\(d) wcc_surface(d$A, d$B) |> mutate(dyad = d$dyad[1], closeness = d$closeness[1])) |>
list_rbind()
metrics_5 <- surfaces_5 |>
drop_na(r) |>
slice_max(r, n = 1, by = c(dyad, start)) |>
summarise(peak_r = mean(r), mean_lag_s = mean(lag) * EPOCH_S,
lag_sd_s = sd(lag) * EPOCH_S, .by = c(dyad, closeness)) |>
left_join(summarise(five, r_lag0 = cor(A, B), .by = dyad), by = "dyad")Session information
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Golden Gate 27.0.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.3.2 lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0
## [5] dplyr_1.2.1 purrr_1.2.2 readr_2.2.0 tidyr_1.3.2
## [9] tibble_3.3.1 ggplot2_4.0.3 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 utf8_1.2.6 generics_0.1.4 stringi_1.8.9
## [5] hms_1.1.4 digest_0.6.39 magrittr_2.0.5 evaluate_1.0.5
## [9] grid_4.6.1 timechange_0.4.0 RColorBrewer_1.1-3 bookdown_0.47
## [13] fastmap_1.2.0 jsonlite_2.0.0 ggrepel_0.9.8 backports_1.5.1
## [17] scales_1.4.0 codetools_0.2-20 jquerylib_0.1.4 cli_3.6.6
## [21] rlang_1.3.0 withr_3.0.3 cachem_1.1.0 yaml_2.3.12
## [25] otel_0.2.0 tools_4.6.1 tzdb_0.5.0 broom_1.0.13
## [29] vctrs_0.7.3 R6_2.6.1 lifecycle_1.0.5 pkgconfig_2.0.3
## [33] pillar_1.11.1 bslib_0.11.0 gtable_0.3.6 glue_1.8.1
## [37] rmdformats_1.0.4 Rcpp_1.1.2 xfun_0.60 tidyselect_1.2.1
## [41] rstudioapi_0.19.0 knitr_1.51 farver_2.1.2 htmltools_0.5.9
## [45] rmarkdown_2.31 labeling_0.4.3 compiler_4.6.1 S7_0.2.2