Time-lagged cross-correlation of 168 observer/target dyads

Design. 28 observers (14 Test, 14 Control) x 6 target trials = 168 dyadic analyses. Continuous tracking sampled at 1000 ms.

Lag window. -5 s to +10 s evaluated at every 1000 ms increment. Both endpoints are included, so each trial yields exactly 16 coefficients.

Sign convention. r(tau) = pearson(x[t], y[t - tau]). A positive peak lag means the observer lags the target, which is the quantity the study calls temporal cognitive lag.

The data in this notebook is synthetic, generated with a known per-observer lag so that every step can be checked against ground truth. No study data is used.

import numpy as np, pandas as pd
from scipy.stats import pearsonr
import statsmodels.formula.api as smf

from tlcc.core import lag_grid, tlcc, peak_lag, fit_random_intercept_lmm
from tlcc.simulate import build_dataset

ds = build_dataset()
lags = lag_grid(-5, 10, step_ms=1000)
print("lags:", list(lags))
print("n coefficients per trial:", len(lags))
print("n dyads:", len(ds["dyads"]))

lags: [np.int64(-5), np.int64(-4), np.int64(-3), np.int64(-2), np.int64(-1), np.int64(0), np.int64(1), np.int64(2), np.int64(3), np.int64(4), np.int64(5), np.int64(6), np.int64(7), np.int64(8), np.int64(9), np.int64(10)]
n coefficients per trial: 16
n dyads: 168

1. The lag grid is 16 points, not 15 or 17

np.arange(-5, 10) silently drops the +10 endpoint and gives 15. lag_grid is endpoint-inclusive by construction, and the assertion below is in the test suite.

assert len(lags) == 16 and lags[0] == -5 and lags[-1] == 10
print("OK: 16 coefficients spanning -5 s .. +10 s inclusive")

OK: 16 coefficients spanning -5 s .. +10 s inclusive

2. At zero lag, the TLCC coefficient is the Pearson r

A cheap but load-bearing check: if this fails, the shift indexing is wrong.

d0 = ds["dyads"][0]
ours = tlcc(d0["x"], d0["y"], np.array([0]))[0]
scipy_r = pearsonr(d0["x"], d0["y"])[0]
print(f"tlcc(lag=0) = {ours:.15f}")
print(f"scipy       = {scipy_r:.15f}")
print("max abs diff:", abs(ours - scipy_r))

tlcc(lag=0) = 0.824175278245446
scipy       = 0.824175278245446
max abs diff: 0.0

3. One dyad, all 16 lags

Correlations at non-zero lag use the overlapping segment only. Zero padding would bias r toward 0 as |tau| grows; circular wrap would invent structure.

coeffs = tlcc(d0["x"], d0["y"], lags)
lag, r = peak_lag(coeffs, lags)
pd.DataFrame({"lag_s": lags, "r": np.round(coeffs, 4)}).T

	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
lag_s	-5.0000	-4.0000	-3.0000	-2.0000	-1.0000	0.0000	1.0000	2.0000	3.0000	4.0000	5.0000	6.0000	7.0000	8.0000	9.000	10.0000
r	0.5428	0.5953	0.6509	0.7085	0.7668	0.8242	0.8769	0.9224	0.9537	0.9641	0.9514	0.9175	0.8678	0.8089	0.746	0.6836

print(f"peak lag = {lag} s   r = {r:.4f}   true injected lag = {d0['true_lag']} s")

peak lag = 4 s   r = 0.9641   true injected lag = 4 s

4. The curve, not just the argmax

The target is a smoothed random walk, so it is strongly autocorrelated and the correlation peak is broad. Reporting the argmax alone hides that. The shaded band is where r is within 0.01 of the peak: several lags are nearly indistinguishable.

import matplotlib.pyplot as plt
plt.rcParams.update({"figure.dpi": 110, "font.size": 9})
fig, ax = plt.subplots(figsize=(7, 2.8))
near = coeffs >= (np.nanmax(coeffs) - 0.01)
ax.fill_between(lags, 0, 1, where=near, transform=ax.get_xaxis_transform(),
                color="#C4462F", alpha=0.10, lw=0)
ax.axvline(0, color="#B9B2A6", lw=0.8)
ax.plot(lags, coeffs, "-o", color="#1B1917", ms=3.5, lw=1.2)
ax.plot([lag], [r], "o", color="#C4462F", ms=7)
ax.set_xlabel("lag (s)  -  positive = observer lags target")
ax.set_ylabel("Pearson r"); ax.set_xticks(lags)
ax.set_title(f"Observer {d0['observer']} / target {d0['target']}: peak at {lag} s")
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout(); plt.show()

Matplotlib is building the font cache; this may take a moment.

5. All 168 dyads, and how often the peak recovers the truth

rows = []
for d in ds["dyads"]:
    c = tlcc(d["x"], d["y"], lags)
    pl, pr = peak_lag(c, lags)
    rows.append({"observer": d["observer"], "target": d["target"], "group": d["group"],
                 "group_code": d["group_code"], "true_lag": d["true_lag"],
                 "peak_lag": pl, "peak_r": pr})
df = pd.DataFrame(rows)
err = (df.peak_lag - df.true_lag).abs()
print(f"dyads analysed      : {len(df)}")
print(f"exact peak recovery : {(err == 0).mean():.4f}")
print(f"within +/- 1 sample : {(err <= 1).mean():.4f}")
print(f"worst error         : {err.max()} sample(s)")
df.head(8)

dyads analysed      : 168
exact peak recovery : 0.9464
within +/- 1 sample : 1.0000
worst error         : 1 sample(s)

	observer	target	group	group_code	true_lag	peak_lag	peak_r
0	0	0	Test	1	4	4	0.964150
1	0	1	Test	1	4	4	0.965023
2	0	2	Test	1	4	5	0.967165
3	0	3	Test	1	4	4	0.964453
4	0	4	Test	1	4	4	0.971515
5	0	5	Test	1	4	4	0.974981
6	1	0	Test	1	3	3	0.971505
7	1	1	Test	1	3	3	0.963569

The argmax misses by one sample on ~5% of dyads and never by more. That is a property of peak-picking on autocorrelated signals, not of the implementation. It is also the reason a study like this should model the peak lag per dyad and carry the uncertainty into the group model, rather than treat each argmax as exact.

6. Group comparison: trials nested within observers

Each observer contributes 6 trials, so the trials are not independent. A random intercept per observer absorbs the between-observer spread:

peak_lag_ij = b0 + b1 * Test_i + u_i + e_ij,     u_i ~ N(0, sigma_re^2)

fit_random_intercept_lmm fits this by maximum likelihood, profiling out the variance ratio. statsmodels fits the same model independently. They must agree.

ours = fit_random_intercept_lmm(df.peak_lag.astype(float), df.observer, df.group_code)
ref = smf.mixedlm("peak_lag ~ group_code", df, groups=df.observer).fit(reml=False)

cmp = pd.DataFrame({
    "ours":        [ours.beta[0], ours.beta[1], ours.se[1], ours.s2_re, ours.s2_resid, ours.loglike],
    "statsmodels": [ref.params["Intercept"], ref.params["group_code"], ref.bse["group_code"],
                    float(ref.cov_re.iloc[0, 0]), ref.scale, ref.llf],
}, index=["b0 (Control mean)", "b1 (Test - Control)", "SE(b1)",
          "sigma_re^2", "sigma_resid^2", "logLik"])
cmp["abs diff"] = (cmp.ours - cmp.statsmodels).abs()
cmp.round(6)

	ours	statsmodels	abs diff
b0 (Control mean)	0.833333	0.833333	0.000000
b1 (Test - Control)	2.035714	2.035714	0.000000
SE(b1)	0.332132	0.332126	0.000005
sigma_re^2	0.763251	0.763226	0.000025
sigma_resid^2	0.053571	0.053572	0.000000
logLik	-54.975016	-54.975016	0.000000

7. APA-style results line

Ground truth for the simulation: Control lag 1.0 s, Test lag 3.0 s, so the true group difference is 2.0 s.

b1, se1, z1, p1 = ours.beta[1], ours.se[1], ours.z[1], ours.p[1]
lo, hi = b1 - 1.96 * se1, b1 + 1.96 * se1
p_txt = "< .001" if p1 < .001 else f"= {p1:.3f}".replace("0.", ".")
print(f"Observers in the Test group showed a longer temporal lag than Control")
print(f"observers, b = {b1:.2f} s, SE = {se1:.2f}, 95% CI [{lo:.2f}, {hi:.2f}],")
print(f"z = {z1:.2f}, p {p_txt}. A random intercept for observer was included to")
print(f"account for the {ours.n_obs // ours.n_groups} trials nested within each of the")
print(f"{ours.n_groups} observers (sigma_re^2 = {ours.s2_re:.3f}, sigma_resid^2 = {ours.s2_resid:.3f}).")
print()
print(f"[simulation ground truth: true difference = 2.00 s]")

Observers in the Test group showed a longer temporal lag than Control
observers, b = 2.04 s, SE = 0.33, 95% CI [1.38, 2.69],
z = 6.13, p < .001. A random intercept for observer was included to
account for the 6 trials nested within each of the
28 observers (sigma_re^2 = 0.763, sigma_resid^2 = 0.054).

[simulation ground truth: true difference = 2.00 s]

8. Mirror export

The brief asks for the notebook and an automatically exported .py mirror. jupyter nbconvert --to script is what produces it, run from build_notebook.py so the two can never drift.