Wavelet Coherence: An Illustration

One therapy dyad, sixteen sessions, and 199 coded moments

library(tidyverse)
library(patchwork)
# biwavelet and signal are called with :: throughout, so it is always visible
# where a computation leaves the tidyverse.

# The wavelet transforms take about two minutes in total. They are written to
# disk the first time and read back afterwards, so re-rendering is immediate.
# Delete the directory to force a recomputation.
DERIVED <- here::here("tutorials", "_derived", "09i")
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")
BANDS  <- tibble(band = c("phasic (2–20 s)", "tonic (20–128 s)"),
                 lo = c(2, 20), hi = c(20, 128))
BAND_FILL <- set_names(c("#c05746", "#7fb3d5"), BANDS$band)

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, which discards both which rhythm the two shared and when they shared it. This illustration recovers both axes: the coupling is estimated separately at each timescale and at each moment, and the result is then 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 and are read here rather than rebuilt. Two decisions were taken in that script 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 the content above the new Nyquist frequency \(f_s/2\) — 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 with a 30-tap FIR kernel. The FIR path is specified rather than the default IIR because the IIR design carries passband ripple, and because both are zero-phase — a filter that shifted the series in time would move every episode alignment with it.

The episode times were read from the per-session marker files, not from the pooled file supplied beside them. The two disagree by one episode, and only the per-session files carry the session identifier unambiguously. One episode of a fifth type occurs exactly once and was dropped.

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 durations govern Step 4 and are worth reading now. The median episode lasts 17 seconds and only 28 of the 199 reach a minute. The types also differ sharply in length — Reconceptualization episodes have 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

Before any estimator, plot the thing. Two sessions are shown: the first, which is sparsely coded, and session 5, the most densely coded of the sixteen.

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 are visible and each returns later. The two signals have different character: the patient’s 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 and dominated by slow movement. Both carry a strong upward drift across the session, which a whole-recording correlation would read as agreement. And the episodes are short and scattered against a 45-minute recording.

---

Step 2 — Choose the frequency band

A frequency-domain method requires a band, and the band is a claim about which rhythms could plausibly carry the coupling of interest. The usual way to choose one is to locate a peak in the spectrum and take the frequencies around it. That is not available here.

spectrum_of <- function(x) {
  s <- spectrum(x - mean(x), spans = c(51, 51), taper = 0.1,
                detrend = TRUE, plot = FALSE)
  tibble(freq = s$freq * FS, power = s$spec)
}

spectra <- cached("spectra", {
  signals |>
    pivot_longer(c(patient, therapist), names_to = "who", values_to = "x") |>
    reframe(spectrum_of(x), .by = c(session, who))
})

shading <- tibble(band = BANDS$band,
                  xmin = 1 / BANDS$hi, xmax = 1 / BANDS$lo) |>
  mutate(band = factor(band, BANDS$band))

spectra |>
  filter(freq >= .002, freq <= 1) |>
  mutate(relative = power / max(power), .by = c(session, who)) |>
  ggplot() +
  geom_rect(data = shading,
            aes(xmin = xmin, xmax = xmax, ymin = 1e-8, ymax = 1.5, fill = band),
            alpha = .3) +
  geom_line(aes(freq, relative, colour = who), linewidth = .25) +
  facet_wrap(~ session, ncol = 8) +
  scale_fill_manual(values = BAND_FILL, name = NULL) +
  scale_colour_manual(values = PERSON, name = NULL) +
  scale_x_log10(breaks = c(.01, .1, 1), labels = c(".01", ".1", "1")) +
  scale_y_log10(breaks = 10^c(-6, -3, 0), labels = c("1e-6", "1e-3", "1"),
                limits = c(1e-8, 1.5)) +
  labs(title = "Spectra of all 16 sessions, each scaled to its own maximum",
       subtitle = paste("Log-log. There is no peak anywhere: power falls",
                        "smoothly from the slowest frequencies upward."),
       x = "frequency (Hz)", y = "relative power")

Every panel shows the same shape. Power declines steadily with frequency and there is no bump at any scale, in either person, in any session. This is not a defect of the recording: skin conductance has no respiratory or baroreflex rhythm to generate a peak, unlike heart period, where the high- and low-frequency bands of heart rate variability are defined by exactly such features. Nor does the power sit anywhere useful.

spectra |>
  mutate(band = cut(freq, c(0, 1/128, 1/20, 1/2, Inf),
                    labels = c("below 0.008 Hz", "tonic, 0.008–0.05 Hz",
                               "phasic, 0.05–0.5 Hz", "above 0.5 Hz"))) |>
  summarise(power = sum(power), .by = c(who, band)) |>
  mutate(`share (%)` = round(100 * power / sum(power), 1), .by = who) |>
  select(who, band, `share (%)`) |>
  pivot_wider(names_from = who, values_from = `share (%)`)

band	patient	therapist
below 0.008 Hz	56.2	64.6
tonic, 0.008–0.05 Hz	37.4	33.7
phasic, 0.05–0.5 Hz	6.4	1.6
above 0.5 Hz	0.0	0.0

More than half of each person’s power lies below 0.008 Hz, at periods longer than two minutes. That is the session-long drift seen in Step 1, and it cannot distinguish one 17-second moment from the next.

So the band must come from the physiology of the signal rather than from its spectrum. Electrodermal activity has an established decomposition into two components, and the division between them is a frequency:

• a phasic component — discrete skin conductance responses, rising in one to three seconds and recovering over several more. The standard high-pass used to isolate it is 0.05 Hz, the default in Biopac’s AcqKnowledge, in NeuroKit2 and in the PhysioData Toolbox. The signal is band-limited above at roughly 0.5 Hz; the original authors of this dataset low-pass filtered individual sessions at 0.45–0.6 Hz for exactly that reason.
• a tonic component — the slowly varying background level, below 0.05 Hz.

This gives a phasic band of 0.05–0.5 Hz, or periods of 2–20 s, and a tonic band which we take as periods of 20–128 s. The lower bound of 128 s is a choice: a 45-minute session contains only about 21 cycles of a 128-second rhythm, and fewer cycles than that would not support an estimate.

Two things follow that should be stated before any result appears. The phasic band holds only 6.4% of the patient’s power and 1.6% of the therapist’s, so if coupling lives there it will be coupling in a small part of the signal. And the phasic band is the one whose timescale is comparable to a 17-second episode; the tonic band is not, and is carried forward partly as a comparison.

---

Step 3 — Estimate the coherence surface

The transform is the Morlet wavelet with \(\omega_0 = 6\). It is complex-valued, so it returns a phase as well as an amplitude at every point in time and scale, which is what a coherence requires — a real-valued wavelet would give amplitude only. The value \(\omega_0 = 6\) is the conventional compromise: narrow enough in frequency to separate neighbouring scales, short enough in time to localise an event to within a few of its own cycles.

Writing \(W_A(a,\tau)\) and \(W_B(a,\tau)\) for the two transforms at scale \(a\) and time \(\tau\), the squared wavelet coherence is

\[R^{2}(a,\tau) \;=\; \frac{\bigl\lvert S\!\left(a^{-1} W_A \overline{W_B}\right)\bigr\rvert^{2}} {S\!\left(a^{-1} \lvert W_A \rvert^{2}\right)\, S\!\left(a^{-1} \lvert W_B \rvert^{2}\right)}\]

where the bar denotes complex conjugation and \(S(\cdot)\) is a smoothing operator applied in both time and scale.

Why the smoothing is not optional

Without \(S(\cdot)\) the expression is identically 1. This is algebra, not approximation: for any two complex numbers, \(\lvert W_A \overline{W_B}\rvert = \lvert W_A\rvert\,\lvert W_B\rvert\), so numerator and denominator are equal at every point. A coherence is a correlation, and a correlation computed from a single observation is always perfect, for the same reason a straight line always passes through one point.

The smoothing operator supplies the sample size. Once each term is a sum over a neighbourhood, the numerator sums complex quantities, which can cancel, while the denominator sums magnitudes, which cannot. Coherence is high when the phase difference between the two people is consistent across that neighbourhood, so that the products align and add; it is low when the lead–lag relationship wanders and they partially cancel.

MIN_PERIOD <- 2; MAX_PERIOD <- 128
MORLET_K   <- 1.033   # Morlet-6: period = 1.033 * scale

smoothing_demo <- cached("smoothing_demo", {
  d  <- filter(signals, session == 5)
  wa <- biwavelet::wt(cbind(d$t, d$patient),
                      s0 = MIN_PERIOD / MORLET_K,
                      max.scale = MAX_PERIOD / MORLET_K, do.sig = FALSE)
  wb <- biwavelet::wt(cbind(d$t, d$therapist),
                      s0 = MIN_PERIOD / MORLET_K,
                      max.scale = MAX_PERIOD / MORLET_K, do.sig = FALSE)
  unsmoothed <- Mod(wa$wave * Conj(wb$wave))^2 /
                (Mod(wa$wave)^2 * Mod(wb$wave)^2)
  smoothed <- biwavelet::wtc(cbind(d$t, d$patient), cbind(d$t, d$therapist),
                             s0 = MIN_PERIOD / MORLET_K,
                             max.scale = MAX_PERIOD / MORLET_K,
                             nrands = 0, quiet = TRUE)$rsq
  tibble(smoothing = c("none", "Torrence–Webster (the standard)"),
         `mean R2` = round(c(mean(unsmoothed), mean(smoothed)), 4),
         `min R2`  = round(c(min(unsmoothed), min(smoothed)), 4),
         `% above 0.9` = round(100 * c(mean(unsmoothed > .9),
                                       mean(smoothed > .9)), 1))
})

smoothing_demo

smoothing	mean R2	min R2	% above 0.9
none	1.0000	1	100.0
Torrence–Webster (the standard)	0.3357	0	1.5

The unsmoothed surface is exactly 1 at every one of its points, on real data.

Which smoothing

The operator is not a free parameter. The standard for wavelet coherence is Torrence and Webster’s, as set out by Grinsted and colleagues, and it is fixed entirely by the properties of the Morlet wavelet: a Gaussian in time whose width is the scale itself, and a boxcar in scale of 0.6 octaves, the decorrelation length of the wavelet. biwavelet::wtc() implements it directly.

The scale-dependence matters. A smoothing window of fixed duration — say 100 seconds, applied uniformly — would average over fifty cycles of a 2-second rhythm while covering less than one cycle of a 128-second rhythm. It would over-smooth one end of the surface and barely smooth the other, leaving slow scales close to the trivial value of 1.

wtc_session <- function(s) {
  d <- filter(signals, session == s)
  biwavelet::wtc(cbind(d$t, d$patient), cbind(d$t, d$therapist),
                 s0 = MIN_PERIOD / MORLET_K,
                 max.scale = MAX_PERIOD / MORLET_K,
                 nrands = 0, quiet = TRUE)
}

derived <- cached("derived", {
  W <- map(1:16, wtc_session)

  # Band-mean series at full resolution. biwavelet returns `coi`, the longest
  # period still trustworthy at each time point, so a band is unusable wherever
  # the cone of influence has risen above it.
  bands <- map2_dfr(W, 1:16, \(w, s) {
    map_dfr(seq_len(nrow(BANDS)), \(k) {
      in_band <- w$period >= BANDS$lo[k] & w$period < BANDS$hi[k]
      tibble(session = s, band = BANDS$band[k], t = w$t,
             coh  = colMeans(w$rsq[in_band, , drop = FALSE]),
             edge = w$coi < BANDS$hi[k])
    })
  })

  # Thinned surfaces for plotting. The full set is ~450 MB; at 5 s resolution
  # it is finer than the pixels of the figure below.
  surfaces <- map2_dfr(W, 1:16, \(w, s) {
    keep <- seq(1, length(w$t), by = 50)
    expand_grid(i = keep, period = w$period) |>
      mutate(session = s, t = w$t[i], coh = as.numeric(w$rsq[, keep]),
             edge = period > w$coi[i]) |>
      select(session, t, period, coh, edge)
  })

  list(bands = bands, surfaces = surfaces)
})

bands    <- derived$bands
surfaces <- derived$surfaces

A final feature of the surface has to be marked before it can be read. Near either end of a session a slow wavelet extends past the available data, and the estimate there is partly made of nothing. That region is the cone of influence. It widens with the period, so it costs the tonic band far more than the phasic one.

bands |>
  summarise(`% of session inside the cone` = round(100 * mean(edge), 1),
            .by = band)

band	% of session inside the cone
phasic (2–20 s)	1.9
tonic (20–128 s)	12.1

ggplot(surfaces, aes(t / 60, period, fill = coh)) +
  geom_raster(interpolate = TRUE) +
  geom_tile(data = filter(surfaces, edge), fill = "white", alpha = .55) +
  geom_hline(yintercept = 20, colour = "grey10", linetype = 2, linewidth = .35) +
  facet_wrap(~ session, ncol = 2, labeller = labeller(session = \(s) paste("Session", s))) +
  scale_y_continuous(trans = "log2", breaks = c(2, 8, 32, 128), expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_fill_viridis_c(option = "magma", limits = c(0, 1), name = expression(R^2)) +
  labs(title = "Wavelet coherence surfaces, all 16 sessions",
       subtitle = paste("Pale overlay is the cone of influence. The dashed line",
                        "at 20 s divides the phasic band below from the tonic band above."),
       x = "minutes into session", y = "period (s)")

The surfaces are structured, and structured differently at the two ends. Above the dashed line the pattern is broad: patches of high coherence several minutes across, with the cone of influence removing the upper corners where a 128-second wavelet runs off the end of the recording. Below the line the pattern is fine — narrow vertical filaments, bright and dark alternating on a scale of tens of seconds. That is the only part of the surface whose grain is comparable to a 17-second episode.

The two bands nonetheless average to almost the same value.

bands |>
  filter(!edge) |>
  summarise(`mean R2` = round(mean(coh), 3),
            `sd across the session` = round(sd(coh), 3), .by = band)

band	mean R2	sd across the session
phasic (2–20 s)	0.315	0.117
tonic (20–128 s)	0.316	0.113

---

Step 4 — Read the surface at the coded moments

A 17-second episode is too short to support an estimate in a band whose slowest rhythm is 20 seconds long, and far too short for one whose slowest is 128. The resolution limit is exact: over a stretch of length \(T\) seconds, frequencies closer together than \(\Delta f = 1/T\) cannot be distinguished at all.

The way through is not to estimate anything inside an episode. The surface is computed over the whole 45-minute session, where there is ample data, and then sampled at the episode’s times. Let episode \(e\) span the \(n_e\) time points \(\tau_1, \dots, \tau_{n_e}\), and let the band contain the \(m\) scales \(a_1, \dots, a_m\). Then

\[\bar{R}^{2}_{e} \;=\; \frac{1}{n_e \, m} \sum_{j=1}^{n_e} \sum_{k=1}^{m} R^{2}(a_k, \tau_j)\]

with every cell inside the cone of influence excluded from the sum. What this yields is not an estimate made from the episode but the value of a session-scale estimate at the episode’s location, and it is smoothed across a neighbourhood wider than the episode itself.

S <- 5
bands |>
  filter(session == S) |>
  ggplot(aes(t / 60, coh)) +
  geom_rect(data = filter(episodes, session == S), inherit.aes = FALSE,
            aes(xmin = start / 60, xmax = end / 60,
                ymin = -Inf, ymax = Inf, fill = type), alpha = .6) +
  geom_line(aes(alpha = !edge), linewidth = .3, colour = "grey15") +
  facet_wrap(~ band, ncol = 1) +
  scale_alpha_manual(values = c(`FALSE` = .18, `TRUE` = 1), guide = "none") +
  scale_fill_manual(values = TYPE, limits = TYPES, name = NULL, drop = FALSE) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(title = sprintf("Band-mean coherence across session %d", S),
       subtitle = "Faded where the cone of influence makes the estimate unreliable.",
       x = "minutes into session", y = expression(mean~R^2))

\(R^2\) is bounded between 0 and 1 and reads directly as the proportion of variance shared at that scale and moment. The question is not what it equals during a coded moment but whether that differs from an uncoded one, which requires a duration-matched control window: same session, same length, drawn at random from a stretch no coder marked. Five are drawn per episode and averaged.

window_mean <- function(d, s, from, to) {
  v <- d |> filter(session == s, t >= from, t < to, !edge) |> pull(coh)
  if (length(v) < 2) NA_real_ else mean(v)
}

matched <- cached("matched", {
  set.seed(11)
  map_dfr(BANDS$band, \(b) {
    d <- filter(bands, band == b)
    map_dfr(seq_len(nrow(episodes)), \(i) {
      e <- episodes[i, ]
      occupied <- filter(episodes, session == e$session)
      t_max <- max(d$t[d$session == e$session])
      draws <- map_dbl(1:5, \(.) {
        for (try in 1:40) {
          from  <- runif(1, 0, t_max - e$duration)
          clash <- any(from < occupied$end &
                       (from + e$duration) > occupied$start)
          if (!clash) {
            v <- window_mean(d, e$session, from, from + e$duration)
            if (!is.na(v)) return(v)
          }
        }
        NA_real_
      })
      tibble(id = e$id, band = b, type = e$type, duration = e$duration,
             episode = window_mean(d, e$session, e$start, e$end),
             control = mean(draws, na.rm = TRUE))
    })
  }) |>
    filter(!is.na(episode), !is.na(control)) |>
    mutate(gap = episode - control)
})

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

band	n	episode	control	difference	t	p
phasic (2–20 s)	189	0.314	0.318	-0.0042	-0.54	0.590
tonic (20–128 s)	176	0.318	0.311	0.0063	0.72	0.475

Coded episodes do not differ from control windows. In the phasic band the episode mean is \(0.314\) against a control of \(0.318\), a difference of \(-0.004\), \(t = -0.54\), \(p = .59\). In the tonic band it is \(0.318\) against \(0.311\), \(+0.006\), \(p = .48\). Both differences are smaller than a hundredth of a unit on a scale that runs from 0 to 1.

A mean and an interval conceal what that comparison looks like, so the two distributions are given in full.

dists <- matched |>
  select(band, episode, control) |>
  pivot_longer(c(episode, control), names_to = "window", values_to = "coh") |>
  mutate(window = factor(window, c("episode", "control")))

panel_a <- ggplot(dists, aes(coh, fill = window)) +
  geom_histogram(bins = 40, position = "identity", alpha = .55, colour = NA) +
  facet_wrap(~ band, ncol = 2, scales = "free_x") +
  scale_fill_manual(values = c(episode = "#c05746", control = "grey55"),
                    name = NULL) +
  labs(title = "A. Where the two distributions sit",
       subtitle = "Coded episodes and their duration-matched control windows.",
       x = expression(mean~R^2), y = "episodes")

annot <- matched |>
  summarise(mean = mean(gap), sd = sd(gap), .by = band) |>
  mutate(txt = sprintf("mean %+.4f\nSD %.3f\n(%.0f× the mean)",
                       mean, sd, sd / abs(mean)))

panel_b <- 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(~ band, ncol = 2) +
  labs(title = "B. The per-episode difference, episode minus control",
       subtitle = "Black line is zero; red line is the mean difference.",
       x = expression(Delta~mean~R^2), y = "episodes")

panel_a / panel_b

The distributions are nearly coincident, and panel B gives the reason the mean is uninformative on its own: the episode-to-episode standard deviation of the difference is \(0.107\) in the phasic band and \(0.117\) in the tonic, roughly twenty-five and nineteen times the mean differences being tested. Whatever drives coherence up or down from one moment to the next, it is an order of magnitude larger than anything the coding distinguishes.

One asymmetry in panel A is an artefact of the construction rather than a result: the control averages five draws while each episode is a single window, so the control distribution is the narrower of the two by design. This is the right choice for comparing means and the wrong one for comparing spreads.

Does the episode difference need a multilevel model?

It does, because the \(t\)-tests above treat 199 episodes as independent and they are not — they arrive in clusters of seven to twenty-three from fifteen sessions. A random intercept for session is the minimum correction.

library(lme4)
library(lmerTest)

mlm <- matched |> left_join(select(episodes, id, session), by = "id")

map_dfr(BANDS$band, \(b) {
  fit <- lmer(gap ~ 1 + (1 | session), data = filter(mlm, band == b))
  co  <- summary(fit)$coefficients
  vc  <- as.data.frame(VarCorr(fit))
  tibble(band = b,
         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),
         `session SD` = round(sqrt(vc$vcov[1]), 4),
         ICC = round(vc$vcov[1] / sum(vc$vcov), 3))
})

band	estimate	SE	df	t	p	session SD	ICC
phasic (2–20 s)	-0.0043	0.0080	14.6	-0.54	0.594	0.0062	0.003
tonic (20–128 s)	0.0057	0.0102	17.2	0.56	0.585	0.0195	0.027

The correction changes nothing, and the intraclass correlations say why: \(0.003\) in the phasic band and \(0.027\) in the tonic. Almost none of the variance in the episode-minus-control difference lies between sessions — it is nearly all within. The degrees of freedom fall from 188 to about 15, which is the honest accounting, and the estimates are unmoved at \(-0.004\) and \(+0.006\).

This is worth stating plainly because the conclusion could have gone the other way: the concern about non-independence was legitimate, and testing it showed the clustering to be negligible. What a session random effect cannot repair is the level above it. All sixteen sessions come from one dyad, and there is no between-dyad variance to model because there is no second dyad.

Adding episode type to the model leaves the picture unchanged. Action episodes sit \(+0.038\) above the others in the phasic band (\(p = .048\)), but the omnibus test for type does not reach significance (\(F(3, 185) = 2.00\), \(p = .116\)), and the model is a singular fit with the session variance driven to zero. On \(n = 44\) Action episodes from a single dyad, inspected after the fact, this is a hypothesis rather than a result.

A final feature of the construction: because the cone of influence widens with period, it removes more of the tonic band than the phasic one, costing 23 episodes against 10.

---

Step 5 — Test against a null

The question Step 4 cannot answer on its own is whether these two people are coupled at all. A mean \(R^2\) of \(0.315\) has no interpretation in isolation: two unrelated skin conductance series will also share some variance at some scales, because both are slow, autocorrelated, drifting signals.

A surrogate must destroy the coupling while preserving everything else. With a single dyad there is no second dyad to shuffle against, but this pair was recorded sixteen times, which supplies the needed comparison: pair the patient’s session \(i\) with the therapist’s session \(j \neq i\). Each person keeps their own spectrum, their own autocorrelation and their own session-long drift. What is destroyed is only the shared timing. Three mismatched partners are drawn for each of the sixteen sessions.

null_a <- cached("null_a", {
  mismatched <- function(s_patient, s_therapist) {
    a <- filter(signals, session == s_patient)
    b <- filter(signals, session == s_therapist)
    n <- min(nrow(a), nrow(b))
    w <- biwavelet::wtc(cbind(a$t[1:n], a$patient[1:n]),
                        cbind(a$t[1:n], b$therapist[1:n]),
                        s0 = MIN_PERIOD / MORLET_K,
                        max.scale = MAX_PERIOD / MORLET_K,
                        nrands = 0, quiet = TRUE)
    map_dfr(seq_len(nrow(BANDS)), \(k) {
      in_band <- w$period >= BANDS$lo[k] & w$period < BANDS$hi[k]
      usable  <- w$coi >= BANDS$hi[k]
      tibble(band = BANDS$band[k],
             coh = mean(colMeans(w$rsq[in_band, , drop = FALSE])[usable]))
    })
  }
  set.seed(2)
  pairs <- map_dfr(1:16, \(s) tibble(p = s, t = sample(setdiff(1:16, s), 3)))
  map2_dfr(pairs$p, pairs$t, mismatched)
})

real_a <- bands |> filter(!edge) |> summarise(coh = mean(coh), .by = c(session, band))

map_dfr(BANDS$band, \(b) {
  r <- filter(real_a, band == b)$coh
  m <- filter(null_a, band == b)$coh
  test <- t.test(r, m)
  tibble(band = b,
         `real (n=16)` = round(mean(r), 3), `sd` = round(sd(r), 3),
         `mismatched (n=48)` = round(mean(m), 3), `sd ` = round(sd(m), 3),
         difference = round(mean(r) - mean(m), 4),
         t = round(test$statistic, 2), p = round(test$p.value, 4))
})

band	real (n=16)	sd	mismatched (n=48)	sd	difference	t	p
phasic (2–20 s)	0.315	0.010	0.303	0.011	0.0125	4.16	0.0003
tonic (20–128 s)	0.316	0.028	0.301	0.027	0.0152	1.91	0.0685

bind_rows(mutate(real_a, kind = "real dyad, same session"),
          mutate(null_a, kind = "different sessions")) |>
  ggplot(aes(coh, kind, fill = kind)) +
  geom_boxplot(width = .5, outlier.size = .7, linewidth = .3) +
  facet_wrap(~ band, ncol = 2, scales = "free_x") +
  scale_fill_manual(values = c("real dyad, same session" = "#c05746",
                               "different sessions" = "grey78"), guide = "none") +
  labs(title = "Whole-session mean coherence, real pairings against mismatched ones",
       x = expression(mean~R^2), y = NULL)

In the phasic band the real pairings exceed the mismatched ones, \(0.315\) against \(0.303\), a difference of \(+0.013\) with \(t = 4.16\) and \(p = .0003\). These two people are coupled at the timescale of skin conductance responses, and the coupling is specific to being in the room together on the same day.

In the tonic band the difference is of similar size, \(+0.015\), but \(t = 1.91\) and \(p = .069\). The reason is visible in the spread: the between-session standard deviation is \(0.028\) in the tonic band against \(0.010\) in the phasic. A 45-minute session contains hundreds of cycles of a 10-second rhythm and about twenty of a 128-second one, so the tonic estimate rests on far fewer independent observations and is correspondingly noisier. The tonic result is suggestive and unresolved; the phasic result is not in doubt.

A second surrogate, not run here

The mismatch test answers whether the dyad is coupled. It does not test the Step 4 question — whether coupling concentrates at the coded moments — and the \(t\)-tests reported there rest on an independence assumption that the data violate.

The appropriate surrogate for that question is a circular shift of the markers. Hold both physiological series exactly as they are and move all of a session’s episode windows by a common random offset, wrapping at the end of the session. The physiology is untouched, the number of episodes is untouched, and their durations and spacing relative to one another are untouched; the only thing broken is the alignment between the codes and the signal. Repeating this a thousand times gives a null distribution for the mean episode coherence that requires no assumption of independence, because each draw preserves whatever dependence the real episodes have. Comparing the observed mean against that distribution yields a p-value directly, as the proportion of shifts reaching the observed value.

It would be the natural next step, and it is inexpensive: the surfaces are already computed, so each draw is a resampling of windows rather than a new transform.

---

What the analysis supports

1. The dyad is coupled in the phasic band. Real pairings exceed mismatched pairings at periods of 2–20 s, \(p = .0003\). This is the band in which skin conductance responses occur.
2. The tonic band is unresolved, \(p = .069\), with three times the between-session variance of the phasic band.
3. Coupling is not elevated at the coded moments in either band (\(p = .59\), \(p = .48\)), subject to the independence caveat above.

The third 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 tightly these two autonomic signals share their rhythms. A measure sensitive to amplitude rather than to coupling, or to one person rather than to the pair, could still separate them.

Would a different choice have changed the answer?

The band would, and it is the decision that matters most. Had the band been chosen where the power is, more than half of it lies below 0.008 Hz at periods longer than two minutes, and the analysis would have measured shared drift. Had it been chosen from the heart rate literature, the high- and low-frequency bands of heart rate variability sit at 0.15–0.4 Hz and 0.04–0.15 Hz, which straddle the phasic–tonic boundary of a signal that has neither a respiratory nor a baroreflex rhythm.

The smoothing would not, because it is not a choice: the Torrence–Webster operator is fixed by the wavelet. Had it been a free parameter, it would have mattered a great deal — the same surface runs from a mean of \(1.000\) unsmoothed to \(0.336\) under the standard, and a fixed-duration window would have inflated the tonic band specifically.

The remaining choices are minor. The control-matching rule averages five draws where one would be noisier but unbiased; the decimation to 10 Hz discards content above 5 Hz, where these signals carry no measurable power.

What no choice repairs is the sample. One dyad is one dyad.

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
Wavelet and parameter	Morlet, \(\omega_0 = 6\)
Period range estimated	2–128 s
Scale resolution	dj = 1/12, the package default
Smoothing	Torrence–Webster; Gaussian in time with width equal to the scale, boxcar of 0.6 octaves in scale
Bands	Phasic 2–20 s and tonic 20–128 s, from the EDA decomposition
Cone of influence	Excluded, using the wavelet’s own cone
Episode aggregation	Mean of the surface within the episode window
Comparison	Five duration-matched control windows per episode
Episode model	Random intercept for session, fitted with lme4::lmer
Surrogate	48 cross-session pairings
Software	R version 4.6.1 (2026-06-24); biwavelet 0.20.22; 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 here 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 raise power at high frequencies, which is where the phasic band sits and where the Step 5 effect is, so that effect should be re-checked with the log applied.

The episode comparison has no between-dyad level. A session random intercept corrects for clustering within the sixteen recordings and changes nothing (ICC \(\le 0.027\)), but the episodes ultimately come from one relationship, and no model can supply variance that was never sampled. The marker-shift surrogate described in Step 5 would additionally remove the assumption that episodes are placed independently of the signal.

Skin conductance indexes sympathetic activity only. It says nothing about parasympathetic or cardiac coupling, which may behave differently in the same recordings.

---

References

Torrence, C., & Compo, G. P. (1998). A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 79, 61–78.

Grinsted, A., Moore, J. C., & Jevrejeva, S. (2004). Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Processes in Geophysics, 11, 561–566.

Boucsein, W., et al. (2012). Publication recommendations for electrodermal measurements. Psychophysiology, 49, 1017–1034.