# Setup: dependencies, the shared helpers, non-interactive plotting.
import importlib.util
import subprocess
import sys
import time
import warnings
from pathlib import Path

# 1. Dependencies are pinned in notebooks/requirements.txt.  statsmodels is the one this notebook
#    needs beyond the core stack; if it is missing the notebook falls back to a documented
#    two-stage estimator and says so in every cell that uses it.
_needed = ("mne", "scipy", "matplotlib", "pandas", "pooch")
_missing = [p for p in _needed if importlib.util.find_spec(p) is None]
if _missing:
    _req = next((d / "requirements.txt" for d in (Path.cwd(), *Path.cwd().parents)
                 if (d / "requirements.txt").exists()), None)
    _cmd = [sys.executable, "-m", "pip", "install", "-q"]
    _cmd += ["-r", str(_req)] if _req else ["mne==1.10.2", "pooch>=1.8"]
    subprocess.check_call(_cmd)

HAVE_STATSMODELS = importlib.util.find_spec("statsmodels") is not None

# 2. Shared helpers, located relative to the working directory -- never through an absolute path.
_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
                if (d / "_shared" / "helpers_l6.py").exists()), None)
if _shared is None:
    raise FileNotFoundError("start the kernel in notebooks/L6/ (or notebooks/) so that _shared/helpers_l6.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l3 as L3
import helpers_l6 as L6

# 3. Plotting: static PNGs through the inline backend; every figure cell ends with plt.show().
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import mne
from scipy import stats

mne.set_log_level("WARNING")
plt.rcParams["figure.dpi"] = 72

SEED = L6.SEED
ALPHA = L6.ALPHA
FULL_RUN = False          # this notebook has no permutation test; FULL_RUN only widens the cohort
SUBJECTS = list(L6.SUBSET_DEFAULT) if not FULL_RUN else list(range(1, 41))

if HAVE_STATSMODELS:
    import statsmodels
    import statsmodels.formula.api as smf
    print(f"PATH TAKEN: statsmodels {statsmodels.__version__} is present -- every model below is a real "
          f"MixedLM fit.")
else:
    print("PATH TAKEN: statsmodels is NOT present.  Every mixed model below falls back to a documented "
          "two-stage estimator (fit the trial-level regression inside each subject, then test the "
          "coefficients across subjects).  For a balanced design the two-stage fixed effect is very close "
          "to the mixed model's, but the variance components, the crossed-item model and the likelihood "
          "comparisons are NOT available and those cells say so.")
print(f"MNE {mne.__version__}; helpers_l6 imported from notebooks/_shared")
print(f"ERP CORE cache: {L3.erpcore_root().name}/ (resolved relative to the working directory, or "
      f"$EEG_COURSE_ERPCORE); each subject's EEGLAB pair is deleted as soon as its measurements exist")

before = L6.disk_report("before any download")
t0 = time.time()
cohort = L6.load_cohort(SUBJECTS)
print(f"{len(cohort)} subjects in {time.time() - t0:.0f} s")
L6.disk_report("after the download loop")

i_prespec = L6.path_index(L6.PRESPECIFIED)
MIN_TRIALS = int(L6.FIXED_CHOICES["min_trials_per_condition"])
rows = []
for s in sorted(cohort):
    p = cohort[s]
    keep = p["keep"][i_prespec]
    n_t = int((keep & p["is_target"]).sum())
    n_s = int((keep & ~p["is_target"]).sum())
    if min(n_t, n_s) < MIN_TRIALS:
        print(f"EXCLUDED {s}: {n_t} target / {n_s} standard trials survive, below the stated minimum of "
              f"{MIN_TRIALS}")
        continue
    idx = np.where(keep)[0]
    rows.append(pd.DataFrame({
        "subject": s,
        "trial": p["trial"][idx],                              # position in the session, 0-based
        "condition": np.where(p["is_target"][idx], "target", "standard"),
        "letter": [chr(ord("A") + c % 10 - 1) for c in p["stim_code"][idx]],    # the letter shown
        "block_target": [chr(ord("A") + c // 10 - 1) for c in p["stim_code"][idx]],
        "amp_uv": p["amp"][i_prespec][idx].astype(float),
        "rt_ms": p["rt_ms"][idx],
        "accuracy": p["accuracy"][idx],
    }))
df = pd.concat(rows, ignore_index=True)
df["condition"] = pd.Categorical(df["condition"], categories=["standard", "target"])   # standard = reference
df["is_target"] = (df["condition"] == "target").astype(float)
df["trial_z"] = (df["trial"] - df["trial"].mean()) / df["trial"].std()
print(f"\n{len(df):,} single trials from {df.subject.nunique()} subjects, "
      f"{df.letter.nunique()} stimulus letters")
print(df.groupby("condition", observed=True)["amp_uv"].agg(["count", "mean", "std"]).round(3))
print(f"\ntrials per subject (target / standard):")
per = df.pivot_table(index="subject", columns="condition", values="amp_uv", aggfunc="count", observed=True)
print(f"   target   min {int(per['target'].min())}, median {per['target'].median():.0f}, "
      f"max {int(per['target'].max())}")
print(f"   standard min {int(per['standard'].min())}, median {per['standard'].median():.0f}, "
      f"max {int(per['standard'].max())}")
print(df.head(4).to_string(index=False))

means = df.pivot_table(index="subject", columns="condition", values="amp_uv", aggfunc="mean", observed=True)
counts = df.pivot_table(index="subject", columns="condition", values="amp_uv", aggfunc="count", observed=True)
sds = df.pivot_table(index="subject", columns="condition", values="amp_uv", aggfunc="std", observed=True)
diff = (means["target"] - means["standard"]).to_numpy()
n_sub = len(diff)
t_avg, p_avg = stats.ttest_1samp(diff, 0)
ci_avg = stats.t.ppf(1 - ALPHA / 2, n_sub - 1) * diff.std(ddof=1) / np.sqrt(n_sub)
print(f"AVERAGE-THEN-TEST (the Level 3 analysis), {n_sub} subjects:")
print(f"   effect {diff.mean():+.4f} uV, 95% CI [{diff.mean() - ci_avg:+.4f}, {diff.mean() + ci_avg:+.4f}], "
      f"t({n_sub - 1}) = {t_avg:.3f}, p = {p_avg:.3g}, dz = {diff.mean() / diff.std(ddof=1):.3f}")
print()
# the standardized measurement error of each subject's difference score (L3.5)
sme = np.sqrt(sds["target"] ** 2 / counts["target"] + sds["standard"] ** 2 / counts["standard"])
print(f"{'subject':9s} {'nT':>4s} {'nS':>4s} {'difference':>11s} {'SME':>7s} "
      f"{'implied weight':>15s}")
w = 1 / sme ** 2
for s in means.index:
    print(f"{s:9s} {int(counts.loc[s, 'target']):4d} {int(counts.loc[s, 'standard']):4d} "
          f"{means.loc[s, 'target'] - means.loc[s, 'standard']:+11.3f} {sme[s]:7.3f} "
          f"{w[s] / w.sum():15.3f}")
print(f"\nThe paired t-test gives every subject the weight 1/{n_sub} = {1 / n_sub:.3f}.  Weighted by "
      f"measurement precision they would range from {(w / w.sum()).min():.3f} to {(w / w.sum()).max():.3f} -- "
      f"a factor of {(w.max() / w.min()):.1f} between the best- and worst-measured subject.")
print(f"Ratio of the noisiest subject's SME to the quietest: {sme.max() / sme.min():.2f}")

fig, axes = plt.subplots(1, 3, figsize=(12.5, 4.0))
axes[0].hist(df.loc[df.condition == "standard", "amp_uv"], bins=60, alpha=0.6, label="standard",
             color="tab:blue", density=True)
axes[0].hist(df.loc[df.condition == "target", "amp_uv"], bins=60, alpha=0.6, label="target",
             color="tab:orange", density=True)
axes[0].axvline(0, color="gray", lw=0.7)
axes[0].set(xlabel="Single-trial P3 amplitude at Pz, 300–600 ms (µV)", ylabel="Density",
            title=f"{len(df):,} single trials, both conditions")
axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3)

for s in means.index:
    axes[1].plot([0, 1], [means.loc[s, "standard"], means.loc[s, "target"]], "o-", color="0.6", lw=0.8, ms=3)
axes[1].plot([0, 1], [means["standard"].mean(), means["target"].mean()], "o-", color="tab:red", lw=2.4, ms=7,
             label="group mean")
axes[1].set_xticks([0, 1], ["standard", "target"])
axes[1].set(ylabel="Mean amplitude (µV)", title=f"Subject means ({n_sub} subjects)")
axes[1].legend(fontsize=8)
axes[1].grid(alpha=0.3)

axes[2].scatter(counts["target"], sme, s=28, color="tab:blue")
for s in means.index:
    axes[2].annotate(s.replace("sub-", ""), (counts.loc[s, "target"], sme[s]), fontsize=6,
                     textcoords="offset points", xytext=(3, 2))
axes[2].set(xlabel="Surviving target trials (count)", ylabel="SME of the difference score (µV)",
            title="Precision differs between subjects,\nand the paired test ignores it")
axes[2].grid(alpha=0.3)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

def two_stage(formula_terms, data, response="amp_uv"):
    """Fallback estimator: fit the trial-level regression inside each subject, then test the
    coefficients across subjects (Holmes & Friston's summary-statistic approach).  Returns the
    fixed-effect table.  It cannot report variance components or a likelihood."""
    import itertools
    coefs = {}
    for s, g in data.groupby("subject"):
        Xc = np.column_stack([np.ones(len(g))] + [g[t].to_numpy(float) for t in formula_terms])
        y = g[response].to_numpy(float)
        beta, *_ = np.linalg.lstsq(Xc, y, rcond=None)
        coefs[s] = beta
    B = np.array(list(coefs.values()))
    out = []
    for j, name in enumerate(["Intercept"] + list(formula_terms)):
        t, p = stats.ttest_1samp(B[:, j], 0)
        out.append({"name": name, "coef": float(B[:, j].mean()),
                    "se": float(B[:, j].std(ddof=1) / np.sqrt(len(B))), "z": float(t), "p": float(p)})
    return pd.DataFrame(out)


def fixed_effect(fit, term="condition[T.target]"):
    """(estimate, standard error, z, p) of one fixed effect, from either path."""
    if HAVE_STATSMODELS and hasattr(fit, "params"):
        return (float(fit.params[term]), float(fit.bse[term]), float(fit.tvalues[term]),
                float(fit.pvalues[term]))
    row = fit.loc[fit["name"] == term].iloc[0]
    return float(row["coef"]), float(row["se"]), float(row["z"]), float(row["p"])


t0 = time.time()
if HAVE_STATSMODELS:
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        m_ri = smf.mixedlm("amp_uv ~ condition", df, groups=df["subject"]).fit(reml=True)
    print(m_ri.summary())
    print(f"\nconverged: {m_ri.converged}; method {m_ri.method}; fitted in {time.time() - t0:.1f} s")
    ri = fixed_effect(m_ri)
else:
    m_ri = two_stage(["is_target"], df)
    print(m_ri.to_string(index=False))
    print("(two-stage fallback: no variance components, no likelihood)")
    ri = fixed_effect(m_ri, "is_target")
print(f"\nRANDOM INTERCEPT -- fixed effect of condition (target minus standard): "
      f"{ri[0]:+.4f} uV, SE {ri[1]:.4f}, z = {ri[2]:.3f}, p = {ri[3]:.3g}")
print(f"the average-then-test estimate on the same trials was {diff.mean():+.4f} uV "
      f"(difference {ri[0] - diff.mean():+.4f} uV)")

t0 = time.time()
if HAVE_STATSMODELS:
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        m_rs = smf.mixedlm("amp_uv ~ condition", df, groups=df["subject"],
                           re_formula="~condition").fit(reml=True)
    print(m_rs.summary())
    print(f"\nconverged: {m_rs.converged}; fitted in {time.time() - t0:.1f} s")
    rs = fixed_effect(m_rs)
    vc = m_rs.cov_re
    print(f"\nrandom-effect covariance (uV^2):")
    print(vc.round(4).to_string())
    sd_int = float(np.sqrt(max(vc.iloc[0, 0], 0)))
    sd_slope = float(np.sqrt(max(vc.iloc[1, 1], 0)))
    corr = float(vc.iloc[0, 1] / (sd_int * sd_slope)) if sd_int > 0 and sd_slope > 0 else float("nan")
    print(f"   SD of the subject intercepts : {sd_int:.3f} uV")
    print(f"   SD of the subject slopes     : {sd_slope:.3f} uV   <- how much the condition effect itself "
          f"varies between people")
    print(f"   their correlation            : {corr:+.3f}")
    print(f"   residual (trial-level) SD    : {np.sqrt(m_rs.scale):.3f} uV")
else:
    m_rs = None
    rs = ri
    sd_int = sd_slope = corr = float("nan")
    print("statsmodels is absent, so there is no random-slope model.  The two-stage fallback already "
          "estimates a separate condition effect inside every subject, which is the same structure fitted "
          "without pooling; its across-subject SD is printed instead.")
    per_subject_slope = df.groupby("subject").apply(
        lambda g: g.loc[g.condition == "target", "amp_uv"].mean()
                  - g.loc[g.condition == "standard", "amp_uv"].mean(), include_groups=False)
    sd_slope = float(per_subject_slope.std(ddof=1))
    print(f"   across-subject SD of the condition effect (no pooling): {sd_slope:.3f} uV")

print(f"\nRANDOM SLOPE -- fixed effect of condition: {rs[0]:+.4f} uV, SE {rs[1]:.4f}, "
      f"z = {rs[2]:.3f}, p = {rs[3]:.3g}")
print(f"\nex-6-2 ANSWER PAIR:")
print(f"   random intercept : {ri[0]:+.4f} uV (SE {ri[1]:.4f}, p = {ri[3]:.3g})")
print(f"   random slope     : {rs[0]:+.4f} uV (SE {rs[1]:.4f}, p = {rs[3]:.3g})")
print(f"   the estimate moves by {rs[0] - ri[0]:+.4f} uV "
      f"({abs(rs[0] - ri[0]) / abs(ri[0]):.1%} of the effect);")
print(f"   the standard error moves by {rs[1] - ri[1]:+.4f} uV "
      f"({(rs[1] / ri[1] - 1):+.1%}), which is the change that matters.")

if HAVE_STATSMODELS:
    lr = 2 * (m_rs.llf - m_ri.llf)
    print(f"Likelihood comparison (REML fits, same fixed effects, so the REML likelihoods ARE comparable):")
    print(f"   random intercept  log-likelihood {m_ri.llf:10.3f}")
    print(f"   random slope      log-likelihood {m_rs.llf:10.3f}")
    print(f"   (statsmodels returns AIC = {m_ri.aic} for a REML fit, because an information criterion "
          f"computed from a REML likelihood is not comparable across models with different fixed effects; "
          f"refit with reml=False if you want AIC)")
    print(f"   2 x difference in log-likelihood: {lr:.3f} on 2 extra parameters")
    print(f"   naive chi-square p = {stats.chi2.sf(max(lr, 0), 2):.4g}")
    print()
    print("   TODO(confirm): that p-value is CONSERVATIVE and should not be quoted as if it were exact.")
    print("   Testing whether a variance is zero puts the null on the boundary of the parameter space, where")
    print("   the likelihood-ratio statistic is not chi-square with 2 df but a mixture.  The usual advice is")
    print("   to decide the random-effects structure from the design rather than from a test: if the design")
    print("   measures the effect repeatedly within subject, the slope belongs in the model whether or not a")
    print("   test says so.  A reviewer should settle which convention this course teaches.")
else:
    print("No likelihood comparison without statsmodels (the two-stage fallback has no likelihood).")

if HAVE_STATSMODELS:
    t0 = time.time()
    df["all"] = 1
    vcf = {"subject": "0 + C(subject)", "letter": "0 + C(letter)"}
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        m_cr = smf.mixedlm("amp_uv ~ condition", df, groups=df["all"], vc_formula=vcf).fit(reml=True)
    cr = fixed_effect(m_cr)
    print(f"converged: {m_cr.converged}; fitted in {time.time() - t0:.1f} s")
    if not m_cr.converged:
        print("   ^ THIS MODEL DID NOT CONVERGE.  Its numbers are printed anyway, and reported as "
              "non-converged in section 7, because a model that fails to converge is a result about the "
              "data rather than something to hide.  Read the variance components below as a description "
              "of where the optimiser stopped, not as estimates.  Section 7 lists what to do about it.")
    print(f"\nvariance components (uV^2):")
    for k, v in zip(vcf, m_cr.vcomp):
        print(f"   {k:9s} {v:9.4f}   (SD {np.sqrt(max(v, 0)):.3f} uV)")
    print(f"   residual  {m_cr.scale:9.4f}   (SD {np.sqrt(m_cr.scale):.3f} uV)")
    print(f"\nCROSSED subject + letter -- fixed effect of condition: {cr[0]:+.4f} uV, SE {cr[1]:.4f}, "
          f"z = {cr[2]:.3f}, p = {cr[3]:.3g}")
    print(f"   against the random-intercept model: estimate {cr[0] - ri[0]:+.4f} uV, "
          f"SE {(cr[1] / ri[1] - 1):+.1%}")
    print()
    print(f"With only {df.letter.nunique()} letters the item variance is estimated from very few levels, so "
          f"whatever it comes out as is imprecise.  The reason to fit it is not the number: it is that "
          f"leaving it out asserts the item variance is exactly zero.")
else:
    cr = ri
    print("No crossed model without statsmodels.")

by_letter = df.pivot_table(index="letter", columns="condition", values="amp_uv",
                           aggfunc=["mean", "count"], observed=True)
print("The condition effect letter by letter (unmodelled, for orientation):")
print(f"{'letter':>7s} {'nT':>6s} {'nS':>6s} {'target':>9s} {'standard':>9s} {'difference':>11s}")
for L in sorted(df.letter.unique()):
    g = df[df.letter == L]
    mt = g.loc[g.condition == "target", "amp_uv"]
    msd = g.loc[g.condition == "standard", "amp_uv"]
    print(f"{L:>7s} {len(mt):6d} {len(msd):6d} {mt.mean():9.3f} {msd.mean():9.3f} "
          f"{mt.mean() - msd.mean():+11.3f}")

d2 = df.dropna(subset=["rt_ms"]).copy()
d2["rt_z"] = d2.groupby("subject")["rt_ms"].transform(lambda x: (x - x.mean()) / x.std())
d2["trial_z"] = d2.groupby("subject")["trial"].transform(lambda x: (x - x.mean()) / x.std())
print(f"{len(d2):,} of {len(df):,} trials have a recorded response time "
      f"({100 * (1 - len(d2) / len(df)):.1f} % have none and are dropped from this model only)")
print(f"reaction time: median {d2.rt_ms.median():.0f} ms, IQR "
      f"{d2.rt_ms.quantile(0.25):.0f}-{d2.rt_ms.quantile(0.75):.0f} ms")

if HAVE_STATSMODELS:
    t0 = time.time()
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        m_st = smf.mixedlm("amp_uv ~ condition * trial_z + rt_z", d2, groups=d2["subject"],
                           re_formula="~condition").fit(reml=True)
    print(f"\nconverged: {m_st.converged}; fitted in {time.time() - t0:.1f} s\n")
    tab = pd.DataFrame({"estimate_uv": m_st.params, "se": m_st.bse, "z": m_st.tvalues, "p": m_st.pvalues})
    print(tab.round(4).to_string())
    st_terms = {k: fixed_effect(m_st, k) for k in
                ("condition[T.target]", "trial_z", "rt_z", "condition[T.target]:trial_z")
                if k in m_st.params.index}
else:
    m_st = two_stage(["is_target", "trial_z", "rt_z"], d2)
    print(m_st.to_string(index=False))
    st_terms = {k: fixed_effect(m_st, k) for k in ("is_target", "trial_z", "rt_z")}

print()
for k, v in st_terms.items():
    print(f"   {k:28s} {v[0]:+8.4f} uV  SE {v[1]:.4f}  z = {v[2]:+7.3f}  p = {v[3]:.3g}"
          + ("   <- significant" if v[3] < ALPHA else ""))
print()
print("Units matter for reading those: trial_z and rt_z are in standard deviations WITHIN a subject, so their")
print("coefficients are 'microvolts per within-subject SD', not microvolts per trial or per millisecond.")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))
for cond, colour in (("standard", "tab:blue"), ("target", "tab:orange")):
    g = df[df.condition == cond]
    b = np.linspace(g.trial.min(), g.trial.max(), 11)
    k = np.digitize(g.trial, b) - 1
    m = [g.amp_uv[k == i].mean() for i in range(10)]
    e = [g.amp_uv[k == i].std(ddof=1) / max(np.sqrt((k == i).sum()), 1) for i in range(10)]
    axes[0].errorbar((b[:-1] + b[1:]) / 2, m, yerr=e, marker="o", ms=4, color=colour, label=cond, capsize=2)
axes[0].axhline(0, color="gray", lw=0.7)
axes[0].set(xlabel="Trial position in the session (count)", ylabel="Mean single-trial amplitude (µV)",
            title="Amplitude over the session (±1 SEM)")
axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3)

g = d2[d2.condition == "target"]
axes[1].scatter(g.rt_ms, g.amp_uv, s=5, alpha=0.25, color="tab:orange")
qs = np.quantile(g.rt_ms, np.linspace(0, 1, 11))
mid = [(qs[i] + qs[i + 1]) / 2 for i in range(10)]
med = [g.amp_uv[(g.rt_ms >= qs[i]) & (g.rt_ms <= qs[i + 1])].mean() for i in range(10)]
axes[1].plot(mid, med, "o-", color="k", lw=1.6, ms=5, label="decile mean")
axes[1].axhline(0, color="gray", lw=0.7)
axes[1].set(xlabel="Reaction time (ms)", ylabel="Single-trial amplitude (µV)",
            title=f"Target trials: amplitude against reaction time (n = {len(g):,})")
axes[1].legend(fontsize=8)
axes[1].grid(alpha=0.3)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

print(f"{'model':>34s} {'converged':>10s} {'effect (uV)':>12s} {'SE':>8s} {'p':>10s}")
fits = [("average-then-test (no model)", True, diff.mean(), diff.std(ddof=1) / np.sqrt(n_sub), p_avg),
        ("random intercept", getattr(m_ri, "converged", None), ri[0], ri[1], ri[3])]
if HAVE_STATSMODELS:
    fits += [("random intercept + slope", m_rs.converged, rs[0], rs[1], rs[3]),
             ("crossed subject + letter", m_cr.converged, cr[0], cr[1], cr[3])]
    k_st = "condition[T.target]"
    fits += [("single-trial regression", m_st.converged, m_st.params[k_st], m_st.bse[k_st], m_st.pvalues[k_st])]
else:
    k_st = "is_target"
    r_ = m_st.loc[m_st["name"] == k_st].iloc[0]
    fits += [("single-trial (two-stage fallback)", None, r_["coef"], r_["se"], r_["p"])]
for name, conv, est, se, p in fits:
    print(f"{name:>34s} {str(conv):>10s} {est:+12.4f} {se:8.4f} {p:10.3g}")
print()
print(f"Every model puts the condition effect within "
      f"{max(abs(f[2] - fits[0][2]) for f in fits):.3f} uV of the average-then-test estimate, and every one "
      f"rejects zero.  The models disagree about PRECISION, not about the effect, which is the usual case for "
      f"a contrast this large and the reason the choice between them has to be made on the design.")

print("nb-6-2-lmm -- L6.2 numbers (draft; TODO(confirm) at author review)")
print(f"Data: ds-erpcore P3, sub-001 to sub-{SUBJECTS[-1]:03d} ({df.subject.nunique()} subjects after the "
      f"{MIN_TRIALS}-trial exclusion rule), {len(df):,} single trials; CC-BY-SA-4.0 per data/directory.yaml")
print(f"Measure: single-trial mean amplitude 300-600 ms at Pz under helpers_l6.PRESPECIFIED "
      f"({L6.describe_path(L6.PRESPECIFIED)})")
print(f"Fitted with: {'statsmodels ' + statsmodels.__version__ + ' MixedLM' if HAVE_STATSMODELS else 'the two-stage fallback (statsmodels absent)'}")
print()
print(f"1. AVERAGE-THEN-TEST (the Level 3 analysis): {diff.mean():+.4f} uV, "
      f"95% CI [{diff.mean() - ci_avg:+.4f}, {diff.mean() + ci_avg:+.4f}], t({n_sub - 1}) = {t_avg:.3f}, "
      f"p = {p_avg:.3g}")
print(f"   what it discards: trial counts ({int(counts['target'].min())}-{int(counts['target'].max())} target "
      f"trials per subject, a {(w.max() / w.min()):.1f}x spread in implied weight), the trial-level variance, "
      f"and the stimulus item.")
print()
print(f"2. ex-6-2 ANSWER PAIR -- does the random slope change the fixed effect?")
print(f"     random intercept          : {ri[0]:+.4f} uV (SE {ri[1]:.4f}, z = {ri[2]:.3f}, p = {ri[3]:.3g})")
print(f"     random intercept + slope  : {rs[0]:+.4f} uV (SE {rs[1]:.4f}, z = {rs[2]:.3f}, p = {rs[3]:.3g})")
print(f"     estimate moves {rs[0] - ri[0]:+.4f} uV; standard error moves {(rs[1] / ri[1] - 1):+.1%}")
if HAVE_STATSMODELS:
    print(f"     subject intercept SD {sd_int:.3f} uV, subject slope SD {sd_slope:.3f} uV, "
          f"correlation {corr:+.3f}, residual SD {np.sqrt(m_rs.scale):.3f} uV")
    print()
    print(f"3. CROSSED subject + letter: {cr[0]:+.4f} uV (SE {cr[1]:.4f}, p = {cr[3]:.3g}); "
          f"variance components " + ", ".join(f"{k} {v:.4f} uV^2" for k, v in zip(vcf, m_cr.vcomp))
          + f"; CONVERGED = {m_cr.converged}"
          + ("  <- do not quote its variance components as estimates" if not m_cr.converged else ""))
    print()
    print(f"4. SINGLE-TRIAL REGRESSION (amp ~ condition * trial_z + rt_z, {len(d2):,} trials with an RT):")
    for k, v in st_terms.items():
        unit = "uV" if k.startswith("condition[T.target]") and ":" not in k else "uV per within-subject SD"
        print(f"     {k:28s} {v[0]:+8.4f} {unit}, SE {v[1]:.4f}, p = {v[3]:.3g}")
print()
print("R equivalents, for a reader who works in lme4:")
print("   lmer(amp_uv ~ condition + (1 | subject), data = df, REML = TRUE)")
print("   lmer(amp_uv ~ condition + (1 + condition | subject), data = df, REML = TRUE)")
print("   lmer(amp_uv ~ condition + (1 | subject) + (1 | letter), data = df, REML = TRUE)")
print("   lmer(amp_uv ~ condition * trial_z + rt_z + (1 + condition | subject), data = df)")
print()
print("Pitfall: pf-group-demographic-confound.  No widget for this lesson.")