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

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

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))
  • stim[cal]
  • stim[standard]
  • stim[target]

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



Total running time of the script: ( 0 minutes 15.773 seconds)

Gallery generated by Sphinx-Gallery