Average ERPs with fitgrid

For designed EEG experiments with categorical variables, a useful range of Hunt-Dawson average ERPs and ERP effects fall out of ordinary least squares regression modeling by selecting the appropriate categorical predictor variable coding with the patsy formula language. The results are identical to the addition, subtraction of average ERP waveforms without programming ad hoc algebraic manipulations.

Prepare and load epochs

import pandas as pd
import fitgrid as fg
from fitgrid import DATA_DIR, sample_data

sample_data.get_file("sub000p3.ms1500.epochs.feather")
p3_epochs_df = pd.read_feather(DATA_DIR / "sub000p3.ms1500.epochs.feather")

# select 3 types of stimulus event: standards, targets, and bioamp calibration triggers
p3_epochs_df = p3_epochs_df.query("stim in ['standard', 'target', 'cal']")

# look up the data QC flags and select the good epochs
good_epochs = p3_epochs_df.query("match_time == 0 and log_flags == 0")[
    "epoch_id"
]
p3_epochs_df = p3_epochs_df.query("epoch_id in @good_epochs")

# rename the time stamp column
p3_epochs_df.rename(columns={"match_time": "time_ms"}, inplace=True)

# select columns of interest for modeling
indices = ["epoch_id", "time_ms"]
predictors = ["stim"]  # categorical with 2 levels: standard, target
channels = ["MiPf", "MiCe", "MiPa", "MiOc"]  # midline electrodes
p3_epochs_df = p3_epochs_df[indices + predictors + channels]

# set the epoch and time column index for fg.Epochs
p3_epochs_df.set_index(["epoch_id", "time_ms"], inplace=True)

# "baseline", i.e., center each epoch on the 200 ms pre-stimulus interval
centered = []
for epoch_id, vals in p3_epochs_df.groupby("epoch_id"):
    centered.append(
        vals[channels]
        - vals[channels].query("time_ms >= -200 and time_ms < 0").mean()
    )
p3_epochs_df[channels] = pd.concat(centered)

# load data into fitgrid.Epochs
p3_epochs_fg = fg.epochs_from_dataframe(
    p3_epochs_df, epoch_id="epoch_id", time="time_ms", channels=channels
)

average ERPs by condition: \(\sim \mathsf{0 + stim}\)

Suppressing the intercept term in the patsy model formula triggers full-rank dummy (indicator) coding of the two-level categorical variable. The estimated coefficients are identical to the average ERPs in each condition. The minimal design matrix illustrates dummy coding for one categorical variable with two levels.

lmg_0_stim = fg.lm(p3_epochs_fg, RHS="0 + stim", quiet=True)

Parameter estimates = Smith & Kutas (2015) regression ERPs

beta_hats = lmg_0_stim.params
beta_hats

		MiPf	MiCe	MiPa	MiOc
time_ms
-748	stim[cal]	0.025192	0.023050	0.021799	-0.041392
	stim[standard]	2.235214	0.628217	2.202564	2.075095
	stim[target]	0.890926	0.909482	3.704045	2.489043
-744	stim[cal]	0.013173	0.029130	-0.018333	-0.015801
	stim[standard]	2.304857	0.758107	2.443701	2.423396
...	...	...	...	...	...
744	stim[standard]	-2.059429	0.405358	-0.481777	-1.198918
	stim[target]	1.159444	-3.025486	2.555300	5.525375
748	stim[cal]	0.037212	0.052357	0.018305	0.012487
	stim[standard]	-1.866571	0.507221	-0.323304	-1.001864
	stim[target]	1.733518	-2.055181	3.420878	5.769298

1125 rows × 4 columns

Parameter estimate standard errors

bses = lmg_0_stim.bse
bses

		MiPf	MiCe	MiPa	MiOc
time_ms
-748	stim[cal]	1.717181	0.850727	0.886257	0.724336
	stim[standard]	1.480024	0.733234	0.763858	0.624299
	stim[target]	3.370163	1.669648	1.739380	1.421591
-744	stim[cal]	1.750222	0.856846	0.901908	0.729748
	stim[standard]	1.508502	0.738508	0.777347	0.628964
...	...	...	...	...	...
744	stim[standard]	0.616983	0.704265	0.760638	0.676497
	stim[target]	1.404932	1.603682	1.732049	1.540452
748	stim[cal]	0.718129	0.818604	0.881669	0.780283
	stim[standard]	0.618949	0.705548	0.759903	0.672520
	stim[target]	1.409409	1.606603	1.730375	1.531395

1125 rows × 4 columns

Visualize parameter estimates +/- standard error

from matplotlib import pyplot as plt

# label index columns for pandas groupby
for attr_df in [beta_hats, bses]:
    attr_df.index.set_names(["time_ms", "beta_hats"], inplace=True)

for beta_hat, vals in beta_hats.groupby("beta_hats"):
    vals.reset_index('beta_hats', inplace=True, drop=True)
    times = vals.index.to_numpy()
    bse = bses.query("beta_hats==@beta_hat")

    fig, ax = plt.subplots(figsize=(8, 4))
    ax.set_title(beta_hat)
    ax.set(
        xlabel="Time (ms)",
        xlim=(times[0], times[-1]),
        ylabel=r"$\mu$V",
        ylim=(-15, 15),
    )
    ax.axhline(0, color="lightgray", lw=1)
    ax.axvline(0, color="gray", lw=1)

    for jdx, chan in enumerate(vals.columns):
        ax.plot(times, vals[chan], label=chan)
        ax.fill_between(
            times, vals[chan] - bse[chan], vals[chan] + bse[chan], alpha=0.2
        )
    ax.legend(loc=(1.05, 0.5))

• 
• 
• 

Why this works

Here is a small (“right hand side”) design matrix for 9 observations of a categorical variable with 3 levels. There is no intercept (constant) and when one of the 3 regressors is 1, the others are 0. The \(\hat{\beta}\) weights that minimize overall error are the means of the data at each level of the categorical variable.

from patsy import demo_data, dmatrix

cat_2 = demo_data("a", nlevels=3, min_rows=8)
dmatrix("0 + a", data=cat_2, return_type="dataframe")

	a[a1]	a[a2]	a[a3]
0	1.0	0.0	0.0
1	0.0	1.0	0.0
2	0.0	0.0	1.0
3	1.0	0.0	0.0
4	0.0	1.0	0.0
5	0.0	0.0	1.0
6	1.0	0.0	0.0
7	0.0	1.0	0.0
8	0.0	0.0	1.0

For EEG data, the “means of the data at each level of the categorical variable” are the time-domain average ERPs. In the sample data, the categorical stimulus variable has three levels: standard, target, and cal for the 10 \(\mu\mathsf{V}\) calibration square wave.

We can reach into one cell of the FitGrid at time = 0 and channel = MiPa and pull out the design matrix. The three column indicator coding is the same as the demo_data example except for the column labels and hundreds observations instead of 9.

lmg_0_stim[0, "MiPa"].model.exog_names.unstack(-2)

		MiPa
		0	1	2
time_ms
0	0	stim[cal]	stim[standard]	stim[target]

lmg_0_stim[0, "MiPa"].model.exog.unstack(-1)

		MiPa
		0	1	2
time_ms	epoch_id
0	0	0.0	0.0	1.0
	1	0.0	0.0	1.0
	2	0.0	0.0	1.0
	3	0.0	0.0	1.0
	4	0.0	0.0	1.0
	...	...	...	...
	596	1.0	0.0	0.0
	597	1.0	0.0	0.0
	598	1.0	0.0	0.0
	599	1.0	0.0	0.0
	600	1.0	0.0	0.0

542 rows × 3 columns