Dyadic Time-Lagged Cross-Correlation Pipeline¶
28 Observers x 6 Targets = 168 dyadic analyses, 1000 ms sampling¶
This notebook is a working preview of the Milestone 1 + Milestone 2 deliverable structure for the HCI continuous-tracking study:
- Data processing (here: synthetic dyadic tracking data, since the real human-subject data is not shared pre-contract).
- Baseline zero-lag Pearson correlations per dyad.
- Time-lagged cross-correlation, lag window -5 s to +10 s, stepped at every 1000 ms increment. That grid is exactly 16 lags, asserted numerically below.
- Peak coefficient + peak lag extraction per dyad.
- Group-level inferential model: linear mixed-effects (statsmodels MixedLM), trials nested within observers, Test vs Control.
- APA-style reporting of the model output.
Every number shown in the execution readouts is computed live in this notebook.
The mirror production script (crosscorr_pipeline_export.py) is auto-exported
from this notebook via jupyter nbconvert --to script.
In [1]:
import json
from pathlib import Path
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf
RNG = np.random.default_rng(42)
# Study design constants (from the project description)
N_OBSERVERS = 28 # split Test / Control
N_TARGETS = 6 # unique target trials per observer
SAMPLE_MS = 1000 # continuous sampling interval
TRIAL_SECONDS = 120 # trial length used for the synthetic preview
LAG_MIN_S, LAG_MAX_S = -5, 10 # uniform lag window, in seconds
The lag grid: -5 s to +10 s at 1000 ms steps = exactly 16 coefficients¶
The window is inclusive on both ends, so the count is (10 - (-5)) + 1 = 16. We assert it rather than trust it.
In [2]:
LAGS_S = np.arange(LAG_MIN_S, LAG_MAX_S + 1) # seconds; step = 1 sample at 1000 ms
assert len(LAGS_S) == 16, f"Expected exactly 16 lags, got {len(LAGS_S)}"
print(f"Lag grid ({len(LAGS_S)} lags, seconds): {LAGS_S.tolist()}")
Lag grid (16 lags, seconds): [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
Synthetic dyadic tracking data¶
28 observers (14 Test, 14 Control) each track 6 unique continuous target trajectories. Targets are smooth band-limited signals. Each observer tracks with a group-dependent response lag and coupling strength, plus motor noise: Test observers are simulated as faster, tighter trackers; Control observers as slower, looser ones. Positive lag means the observer trails the target.
In [3]:
N_SAMPLES = TRIAL_SECONDS # one sample per second at 1000 ms
def make_target(rng: np.random.Generator, n: int) -> np.ndarray:
"""Smooth pseudo-random target trajectory: sum of low-frequency sinusoids."""
t = np.arange(n)
sig = np.zeros(n)
for _ in range(4):
freq = rng.uniform(0.01, 0.06) # cycles per sample
amp = rng.uniform(0.5, 1.5)
phase = rng.uniform(0, 2 * np.pi)
sig += amp * np.sin(2 * np.pi * freq * t + phase)
return sig
def make_observer_trace(rng, target, lag_samples: int, coupling: float,
noise_sd: float) -> np.ndarray:
"""Observer output = lagged copy of the target, scaled + noise."""
lagged = np.roll(target, lag_samples)
# roll wraps; overwrite the wrapped head with the target's first value + noise
lagged[:lag_samples] = target[0]
return coupling * lagged + rng.normal(0, noise_sd, len(target))
GROUPS = {f"O{i+1:02d}": ("Test" if i < 14 else "Control") for i in range(N_OBSERVERS)}
TARGETS = {f"T{j+1}": make_target(RNG, N_SAMPLES) for j in range(N_TARGETS)}
rows, series = [], {}
for obs_id, group in GROUPS.items():
# Group-dependent tracking parameters (per observer, jittered per trial)
base_lag = RNG.normal(2.0, 0.6) if group == "Test" else RNG.normal(4.0, 0.9)
base_coupling = RNG.normal(0.85, 0.05) if group == "Test" else RNG.normal(0.60, 0.08)
for tgt_id, target in TARGETS.items():
lag_s = int(np.clip(round(base_lag + RNG.normal(0, 0.5)), 1, 7))
coupling = float(np.clip(base_coupling + RNG.normal(0, 0.05), 0.2, 0.98))
observer = make_observer_trace(RNG, target, lag_s, coupling, noise_sd=0.55)
rows.append({"observer": obs_id, "group": group, "target": tgt_id,
"gen_lag_s": lag_s})
series[f"{obs_id}|{tgt_id}"] = (target, observer)
dyads = pd.DataFrame(rows)
n_dyads = len(dyads)
assert n_dyads == N_OBSERVERS * N_TARGETS == 168, f"Expected 168 dyads, got {n_dyads}"
print(f"Dyads: {n_dyads} (28 observers x 6 targets); "
f"Test n={sum(g == 'Test' for g in GROUPS.values())} observers, "
f"Control n={sum(g == 'Control' for g in GROUPS.values())} observers")
Dyads: 168 (28 observers x 6 targets); Test n=14 observers, Control n=14 observers
Baseline zero-lag Pearson + time-lagged cross-correlation per dyad¶
For each dyad and each lag k in the 16-lag grid, we shift the observer series by k samples relative to the target and compute Pearson's r on the overlap: positive k tests "observer trails the target by k seconds". Zero-lag r is the baseline. Peak r and its lag are extracted from the 16-coefficient function.
In [4]:
def lagged_pearson(x: np.ndarray, y: np.ndarray, lag: int) -> float:
"""Pearson r between x[t] and y[t + lag], computed on the overlapping window."""
if lag > 0:
a, b = x[:-lag], y[lag:]
elif lag < 0:
a, b = x[-lag:], y[:lag]
else:
a, b = x, y
return float(stats.pearsonr(a, b)[0])
records = []
for row in dyads.itertuples(index=False):
target, observer = series[f"{row.observer}|{row.target}"]
ccf = np.array([lagged_pearson(target, observer, int(k)) for k in LAGS_S])
assert len(ccf) == 16, f"Dyad {row.observer}|{row.target}: {len(ccf)} coefficients"
peak_idx = int(np.argmax(ccf))
records.append({
"observer": row.observer, "group": row.group, "target": row.target,
"gen_lag_s": row.gen_lag_s,
"r_zero_lag": ccf[LAGS_S.tolist().index(0)],
"r_peak": float(ccf[peak_idx]),
"peak_lag_s": int(LAGS_S[peak_idx]),
"ccf": ccf.round(4).tolist(),
})
results = pd.DataFrame(records)
total_coeffs = int(results["ccf"].apply(len).sum())
print(f"Computed {total_coeffs} lagged coefficients "
f"({n_dyads} dyads x {len(LAGS_S)} lags); "
f"all dyads have exactly 16: {bool((results['ccf'].apply(len) == 16).all())}")
# Sanity check: recovered peak lag should match the generating lag
lag_recovery = (results["peak_lag_s"] == results["gen_lag_s"]).mean()
print(f"Peak-lag recovery vs. generating lag: {lag_recovery:.1%} of 168 dyads")
Computed 2688 lagged coefficients (168 dyads x 16 lags); all dyads have exactly 16: True Peak-lag recovery vs. generating lag: 94.6% of 168 dyads
Descriptives: zero-lag baseline vs peak, by group¶
In [5]:
desc = results.groupby("group")[["r_zero_lag", "r_peak", "peak_lag_s"]].agg(["mean", "std"])
print(desc.round(3).to_string())
r_zero_lag r_peak peak_lag_s
mean std mean std mean std
group
Control 0.442 0.169 0.838 0.047 4.060 0.841
Test 0.798 0.100 0.908 0.028 1.821 0.809
Group-level inferential model: linear mixed-effects (MixedLM)¶
Trials are nested within observers (6 repeated dyads per observer), so we fit a linear mixed-effects model with a random intercept per observer. Peak r is Fisher z-transformed before modeling (correlations are not interval-scaled). A second model tests the group difference in peak lag (tracking latency).
In [6]:
results["z_peak"] = np.arctanh(results["r_peak"])
results["group"] = pd.Categorical(results["group"], categories=["Control", "Test"])
m_sync = smf.mixedlm("z_peak ~ group", results, groups=results["observer"]).fit(reml=True)
print(m_sync.summary())
m_lag = smf.mixedlm("peak_lag_s ~ group", results, groups=results["observer"]).fit(reml=True)
print(m_lag.summary())
Mixed Linear Model Regression Results
=======================================================
Model: MixedLM Dependent Variable: z_peak
No. Observations: 168 Method: REML
No. Groups: 28 Scale: 0.0233
Min. group size: 6 Log-Likelihood: 64.2171
Max. group size: 6 Converged: Yes
Mean group size: 6.0
-------------------------------------------------------
Coef. Std.Err. z P>|z| [0.025 0.975]
-------------------------------------------------------
Intercept 1.237 0.022 55.271 0.000 1.193 1.281
group[T.Test] 0.299 0.032 9.461 0.000 0.237 0.361
Group Var 0.003 0.014
=======================================================
Mixed Linear Model Regression Results
========================================================
Model: MixedLM Dependent Variable: peak_lag_s
No. Observations: 168 Method: REML
No. Groups: 28 Scale: 0.3095
Min. group size: 6 Log-Likelihood: -170.7018
Max. group size: 6 Converged: Yes
Mean group size: 6.0
--------------------------------------------------------
Coef. Std.Err. z P>|z| [0.025 0.975]
--------------------------------------------------------
Intercept 4.060 0.179 22.726 0.000 3.709 4.410
group[T.Test] -2.238 0.253 -8.859 0.000 -2.733 -1.743
Group Var 0.395 0.243
========================================================
/Users/mhnd/Documents/upwork/proposal-timeseries-crosscorr/.venv/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space. warnings.warn(msg, ConvergenceWarning)
APA-style results reporting (7th edition statistical style)¶
The write-up below is generated from the fitted models, in the shape a results section expects: descriptives as M and SD, fixed effects as b, SE, z, p, and 95% CI.
In [7]:
def apa_p(p: float) -> str:
return "p < .001" if p < 0.001 else f"p = {p:.3f}".replace("0.", ".")
def apa_stat(v: float, dp: int = 2) -> str:
return f"{v:.{dp}f}".replace("0.", ".", 1) if abs(v) < 1 else f"{v:.{dp}f}"
def mixedlm_apa(m, term: str) -> dict:
ci = m.conf_int().loc[term]
return {"b": float(m.params[term]), "se": float(m.bse[term]),
"z": float(m.tvalues[term]), "p": float(m.pvalues[term]),
"ci_lo": float(ci[0]), "ci_hi": float(ci[1])}
g = results.groupby("group", observed=True)
stats_by_group = {
grp: {"r_peak_m": float(d["r_peak"].mean()), "r_peak_sd": float(d["r_peak"].std()),
"r0_m": float(d["r_zero_lag"].mean()), "r0_sd": float(d["r_zero_lag"].std()),
"lag_m": float(d["peak_lag_s"].mean()), "lag_sd": float(d["peak_lag_s"].std())}
for grp, d in g
}
fx_sync = mixedlm_apa(m_sync, "group[T.Test]")
fx_lag = mixedlm_apa(m_lag, "group[T.Test]")
t, c = stats_by_group["Test"], stats_by_group["Control"]
apa_paragraph = (
f"Across the 168 dyadic trials, peak cross-correlation was higher in the Test group "
f"(M = {apa_stat(t['r_peak_m'])}, SD = {apa_stat(t['r_peak_sd'])}) than in the Control group "
f"(M = {apa_stat(c['r_peak_m'])}, SD = {apa_stat(c['r_peak_sd'])}). A linear mixed-effects model "
f"on Fisher z-transformed peak correlations, with random intercepts for observers, showed a "
f"significant group effect, b = {apa_stat(fx_sync['b'])}, SE = {apa_stat(fx_sync['se'])}, "
f"z = {fx_sync['z']:.2f}, {apa_p(fx_sync['p'])}, 95% CI [{apa_stat(fx_sync['ci_lo'])}, "
f"{apa_stat(fx_sync['ci_hi'])}]. Test observers also tracked with shorter latency "
f"(M = {t['lag_m']:.2f} s, SD = {apa_stat(t['lag_sd'])}) than Control observers "
f"(M = {c['lag_m']:.2f} s, SD = {apa_stat(c['lag_sd'])}), b = {fx_lag['b']:.2f}, "
f"SE = {apa_stat(fx_lag['se'])}, z = {fx_lag['z']:.2f}, {apa_p(fx_lag['p'])}, "
f"95% CI [{fx_lag['ci_lo']:.2f}, {fx_lag['ci_hi']:.2f}]."
)
print(apa_paragraph)
Across the 168 dyadic trials, peak cross-correlation was higher in the Test group (M = .91, SD = .03) than in the Control group (M = .84, SD = .05). A linear mixed-effects model on Fisher z-transformed peak correlations, with random intercepts for observers, showed a significant group effect, b = .30, SE = .03, z = 9.46, p < .001, 95% CI [.24, .36]. Test observers also tracked with shorter latency (M = 1.82 s, SD = .81) than Control observers (M = 4.06 s, SD = .84), b = -2.24, SE = .25, z = -8.86, p < .001, 95% CI [-2.73, -1.74].
Figures¶
Four views of the same computed objects: one raw dyad, its 16-lag cross-correlation function, the group-mean CCFs, and the peak-r distributions.
In [8]:
import matplotlib
try:
get_ipython() # noqa: F821 (defined inside Jupyter, where figures render inline)
except NameError:
matplotlib.use("Agg") # headless backend for plain-script runs
import matplotlib.pyplot as plt
mean_ccf_plot = {grp: np.vstack(d["ccf"].to_list()).mean(axis=0) for grp, d in g}
ex_key = "O03|T2" # example Test dyad
ex_row = results[(results.observer == "O03") & (results.target == "T2")].iloc[0]
tgt, obs = series[ex_key]
fig, axes = plt.subplots(2, 2, figsize=(11, 7))
ax = axes[0, 0]
ax.plot(tgt, lw=1.4, label="Target", color="#1a6fb0")
ax.plot(obs, lw=1.0, label="Observer O03", color="#d1495b", alpha=0.85)
ax.set_title(f"Dyad {ex_key}: raw traces (1000 ms sampling)")
ax.set_xlabel("Time (s)"); ax.legend(frameon=False)
ax = axes[0, 1]
ccf = np.array(ex_row["ccf"])
ax.plot(LAGS_S, ccf, "o-", color="#1a6fb0", lw=1.4)
ax.axvline(0, color="#888", ls=":", lw=1)
ax.plot(ex_row["peak_lag_s"], ex_row["r_peak"], "o", ms=11, mfc="none",
mec="#d1495b", mew=2)
ax.set_title(f"CCF, 16 lags; peak r = {ex_row['r_peak']:.3f} at +{ex_row['peak_lag_s']} s")
ax.set_xlabel("Lag (s), + = observer trails"); ax.set_ylabel("Pearson r")
ax = axes[1, 0]
for grp, color in [("Test", "#1a6fb0"), ("Control", "#d1495b")]:
ax.plot(LAGS_S, mean_ccf_plot[grp], "o-", lw=1.4, color=color, label=grp)
ax.axvline(0, color="#888", ls=":", lw=1)
ax.set_title("Group-mean CCF across 168 dyads")
ax.set_xlabel("Lag (s)"); ax.set_ylabel("Mean r"); ax.legend(frameon=False)
ax = axes[1, 1]
data = [results.loc[results.group == grp, "r_peak"] for grp in ("Test", "Control")]
ax.boxplot(data, tick_labels=["Test", "Control"], widths=0.5)
for i, d in enumerate(data, start=1):
ax.plot(np.full(len(d), i) + RNG.normal(0, 0.05, len(d)), d, ".",
color="#1a6fb0" if i == 1 else "#d1495b", alpha=0.4)
ax.set_title("Peak r by group (168 dyads)"); ax.set_ylabel("Peak Pearson r")
fig.tight_layout()
plt.show()
Export machine-readable results¶
Everything the demo page shows is written here, from this execution: the lag grid, all 168 cross-correlation functions, per-dyad peaks, group descriptives, both MixedLM fixed effects, and the APA paragraph. Nothing is hand-typed.
In [9]:
out_dir = Path(__file__).resolve().parent / "results" if "__file__" in globals() \
else Path.cwd() / "results"
out_dir.mkdir(exist_ok=True)
mean_ccf = {grp: np.vstack(d["ccf"].to_list()).mean(axis=0).round(4).tolist()
for grp, d in g}
payload = {
"design": {"observers": N_OBSERVERS, "targets": N_TARGETS, "dyads": n_dyads,
"sample_ms": SAMPLE_MS, "trial_seconds": TRIAL_SECONDS,
"lag_window_s": [LAG_MIN_S, LAG_MAX_S], "n_lags": int(len(LAGS_S)),
"total_coefficients": total_coeffs,
"lag_recovery_rate": round(float(lag_recovery), 4)},
"lags_s": LAGS_S.tolist(),
"group_stats": stats_by_group,
"mean_ccf_by_group": mean_ccf,
"mixedlm_sync": fx_sync,
"mixedlm_lag": fx_lag,
"apa_paragraph": apa_paragraph,
"dyads": results.drop(columns=["z_peak"]).to_dict(orient="records"),
}
(out_dir / "results.json").write_text(json.dumps(payload, indent=1))
series_payload = {key: {"target": np.round(tgt, 2).tolist(),
"observer": np.round(obs, 2).tolist()}
for key, (tgt, obs) in series.items()}
(out_dir / "series.json").write_text(json.dumps(series_payload))
print(f"Wrote {out_dir / 'results.json'} and series.json "
f"({n_dyads} dyads, {total_coeffs} coefficients)")
Wrote /Users/mhnd/Documents/upwork/proposal-timeseries-crosscorr/analysis/results/results.json and series.json (168 dyads, 2688 coefficients)