Note
Click here to download the full example code
AIC and likelihood ratio model comparison
There are many approaches to comparing models and evaluating relative
goodness-of-fit. The FitGrid[times, channels] query mechanism is designed to streamline
the computations for whatever approach is deemed appropriate for the
research question at hand..
The AIC example demonstrates the fitgrid.utils.summarize function, a convenience wrapper that fits a list of models and returns key summary information for each as an indexed data frame. The summary may be passed to a function for visualizing the AIC \(\mathsf{\Delta_{min}}\) model comparison ([BurAnd2004]) as shown, or processed by custom user routines as needed.
For cases where the summary information is insufficient, the likelihood ratio example illustrates how to compute and visualize a model comparison measure derived from the original fitted grids. Other measures may be computed in a similar manner.
Prepare sample data
import pandas as pd
from matplotlib import pyplot as plt
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")
# drop calibration pulses for these examples
p3_epochs_df = p3_epochs_df.query("stim != '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", "tone"] # 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)
# slice epochs down to a shorter interval
p3_epochs_df.query("time_ms >= -200 and time_ms <= 600", inplace=True)
p3_epochs_df
Ingest as fitgrid.Epochs
p3_epochs_fg = fg.epochs_from_dataframe(
p3_epochs_df, epoch_id="epoch_id", time="time_ms", channels=channels
)
Model summaries
The fitgrid.utils.summarize function is a convenience wrapper that fits a list of models and collects a few key summary measures into a single dataframe indexed by model, beta estimate, and the measure. It supports OLS and LME model fitting and the summaries are returned are in the same format.
This experimental design is a fully crossed 2 (stim: standard, target) \(\times\) 2 (tone: hi, lo). The predictors are treatment coded (patsy default).
Here is a stack of 5 models to summarize and compare:
rhss_T = [
"1 + stim + tone + stim:tone", # long form of "stim * tone"
"1 + stim + tone",
"1 + stim",
"1 + tone",
"1",
]
from fitgrid.utils.summary import summarize
lm_summary_T = summarize(
p3_epochs_fg,
modeler="lm",
RHS=rhss_T,
LHS=channels,
parallel=False,
quiet=True,
)
lm_summary_T
Out:
/home/runner/work/fitgrid/fitgrid/fitgrid/utils/summary.py:145: FutureWarning: fitgrid summaries are in early days, subject to change
warnings.warn(
AIC model comparison: \(\Delta_{i} = \mathsf{AIC}_{i} - \mathsf{AIC_{min}}\)
Akiakie’s information criterion (AIC) increases with residual error and the number of model parameters so comparison on AIC favors the better fitting, more parsimonious models with lower AIC values.
This example visualizes the channel-wise time course of \(\Delta_{i} = \mathsf{AIC}_{i} - \mathsf{AIC_{min}}\), a measure of the AIC of model i vs. the lowest AIC of any model in the set. Burnham and Anderson propose heuristics where models with \(\Delta_{i}\) around 4 are less well supported by the data than the alternative(s) and models with \(\Delta_{i}\) > 7 substantially less so.
In the next figure, the line plots (left column) and raster plots (right column) show the same data in different ways. Higher amplitude line plots and corresponding darker shades of blue in the raster plot indicate that the model’s AIC is higher than the best candidate in the set.
Prestimulus. Prior to stimulus onset at time = 0, the more parsimonious models (bottom three rows) have systematically lower AIC values (broad regions of lighest blue) than the more complex models (top two rows). This indicates that during this interval, the predictor variables alone and in combination do not soak up enough variability to offset the penalty for increasing the model complexity. In terms of this AIC measure, none of the models appear systematically better supported by the data than the null model (bottom row) in the prestimulus interval.
Post-stimulus. In the interval between around 300 - 375 ms poststimulus, the full model that includes stim and tone predictors and their interaction has the minimum AIC among the models compared at all channels except the prefrontal MiPf. The sharp transients in the magntiude of the AIC differences (> 7) at these channels in this interval indicates substantially less support for the alternative models.
from fitgrid.utils.summary import plot_AICmin_deltas
fig, axs = plot_AICmin_deltas(lm_summary_T, figsize=(12, 12))
fig.tight_layout()
for ax_row in range(len(axs)):
axs[ax_row, 0].set(ylim=(0, 50))
Likelihood Ratio
This example fits a full and reduced model, computes and then visualizes the time course of likelihood ratios with a few lines of code.
Fit the full model. The log likelihood dataframe is returned by
querying the FitGrid[times, channels].
lmg_full = fg.lm(p3_epochs_fg, RHS="1 + stim + tone + stim:tone", LHS=channels)
lmg_full.llf
Out:
0%| | 0/201 [00:00<?, ?it/s]
1%|1 | 3/201 [00:00<00:07, 27.12it/s]
3%|2 | 6/201 [00:00<00:06, 28.43it/s]
5%|4 | 10/201 [00:00<00:06, 29.23it/s]
6%|6 | 13/201 [00:00<00:06, 29.45it/s]
8%|7 | 16/201 [00:00<00:06, 28.70it/s]
9%|9 | 19/201 [00:00<00:06, 28.12it/s]
11%|# | 22/201 [00:00<00:06, 28.35it/s]
12%|#2 | 25/201 [00:00<00:06, 28.46it/s]
14%|#3 | 28/201 [00:00<00:06, 28.29it/s]
15%|#5 | 31/201 [00:01<00:06, 28.27it/s]
17%|#6 | 34/201 [00:01<00:05, 28.08it/s]
18%|#8 | 37/201 [00:01<00:05, 28.23it/s]
20%|#9 | 40/201 [00:01<00:05, 28.51it/s]
21%|##1 | 43/201 [00:01<00:05, 28.70it/s]
23%|##2 | 46/201 [00:01<00:05, 28.30it/s]
24%|##4 | 49/201 [00:01<00:05, 28.61it/s]
26%|##5 | 52/201 [00:01<00:05, 28.75it/s]
27%|##7 | 55/201 [00:01<00:05, 28.72it/s]
29%|##8 | 58/201 [00:02<00:05, 28.50it/s]
30%|### | 61/201 [00:02<00:04, 28.76it/s]
32%|###1 | 64/201 [00:02<00:04, 28.80it/s]
33%|###3 | 67/201 [00:02<00:04, 28.89it/s]
35%|###4 | 70/201 [00:02<00:04, 28.99it/s]
36%|###6 | 73/201 [00:02<00:04, 27.99it/s]
38%|###7 | 76/201 [00:02<00:04, 28.33it/s]
40%|###9 | 80/201 [00:02<00:04, 29.13it/s]
42%|####1 | 84/201 [00:02<00:03, 29.56it/s]
43%|####3 | 87/201 [00:03<00:03, 29.08it/s]
45%|####4 | 90/201 [00:03<00:03, 29.22it/s]
46%|####6 | 93/201 [00:03<00:03, 29.26it/s]
48%|####7 | 96/201 [00:03<00:03, 29.35it/s]
49%|####9 | 99/201 [00:03<00:03, 29.32it/s]
51%|##### | 102/201 [00:03<00:03, 29.29it/s]
52%|#####2 | 105/201 [00:03<00:03, 29.15it/s]
54%|#####3 | 108/201 [00:03<00:03, 29.27it/s]
55%|#####5 | 111/201 [00:03<00:03, 29.30it/s]
57%|#####6 | 114/201 [00:03<00:02, 29.14it/s]
58%|#####8 | 117/201 [00:04<00:07, 11.58it/s]
60%|#####9 | 120/201 [00:04<00:05, 14.10it/s]
61%|######1 | 123/201 [00:04<00:04, 16.69it/s]
63%|######2 | 126/201 [00:04<00:03, 19.19it/s]
64%|######4 | 129/201 [00:04<00:03, 21.38it/s]
66%|######5 | 132/201 [00:05<00:02, 23.19it/s]
67%|######7 | 135/201 [00:05<00:02, 24.68it/s]
69%|######8 | 138/201 [00:05<00:02, 25.97it/s]
70%|####### | 141/201 [00:05<00:02, 26.89it/s]
72%|#######1 | 144/201 [00:05<00:02, 27.56it/s]
73%|#######3 | 147/201 [00:05<00:01, 28.09it/s]
75%|#######4 | 150/201 [00:05<00:01, 28.32it/s]
76%|#######6 | 153/201 [00:05<00:01, 28.48it/s]
78%|#######7 | 156/201 [00:05<00:01, 28.87it/s]
79%|#######9 | 159/201 [00:06<00:01, 28.83it/s]
81%|######## | 162/201 [00:06<00:01, 28.80it/s]
82%|########2 | 165/201 [00:06<00:01, 28.87it/s]
84%|########3 | 168/201 [00:06<00:01, 28.35it/s]
85%|########5 | 171/201 [00:06<00:01, 28.76it/s]
87%|########6 | 174/201 [00:06<00:00, 28.91it/s]
89%|########8 | 178/201 [00:06<00:00, 29.26it/s]
91%|######### | 182/201 [00:06<00:00, 29.60it/s]
93%|#########2| 186/201 [00:06<00:00, 29.85it/s]
94%|#########4| 189/201 [00:07<00:00, 29.73it/s]
96%|#########5| 192/201 [00:07<00:00, 29.74it/s]
97%|#########7| 195/201 [00:07<00:00, 29.76it/s]
99%|#########9| 199/201 [00:07<00:00, 29.82it/s]
100%|##########| 201/201 [00:07<00:00, 27.00it/s]
Fit the reduced model likewise.
lmg_reduced = fg.lm(p3_epochs_fg, RHS="1 + stim + tone", LHS=channels)
lmg_reduced.llf
Out:
0%| | 0/201 [00:00<?, ?it/s]
1%|1 | 3/201 [00:00<00:06, 29.17it/s]
3%|3 | 7/201 [00:00<00:06, 30.10it/s]
5%|5 | 11/201 [00:00<00:06, 29.53it/s]
7%|7 | 15/201 [00:00<00:06, 29.98it/s]
9%|9 | 19/201 [00:00<00:05, 30.39it/s]
11%|#1 | 23/201 [00:00<00:05, 30.55it/s]
13%|#3 | 27/201 [00:00<00:05, 30.84it/s]
15%|#5 | 31/201 [00:01<00:05, 31.12it/s]
17%|#7 | 35/201 [00:01<00:05, 31.15it/s]
19%|#9 | 39/201 [00:01<00:05, 31.24it/s]
21%|##1 | 43/201 [00:01<00:04, 31.72it/s]
23%|##3 | 47/201 [00:01<00:04, 31.97it/s]
25%|##5 | 51/201 [00:01<00:04, 31.99it/s]
27%|##7 | 55/201 [00:01<00:04, 32.03it/s]
29%|##9 | 59/201 [00:01<00:04, 32.12it/s]
31%|###1 | 63/201 [00:02<00:04, 32.03it/s]
33%|###3 | 67/201 [00:02<00:04, 32.08it/s]
35%|###5 | 71/201 [00:02<00:04, 32.21it/s]
37%|###7 | 75/201 [00:02<00:04, 31.47it/s]
39%|###9 | 79/201 [00:02<00:03, 30.92it/s]
41%|####1 | 83/201 [00:02<00:03, 31.02it/s]
43%|####3 | 87/201 [00:02<00:03, 29.87it/s]
45%|####5 | 91/201 [00:02<00:03, 30.31it/s]
47%|####7 | 95/201 [00:03<00:03, 30.63it/s]
49%|####9 | 99/201 [00:03<00:03, 30.83it/s]
51%|#####1 | 103/201 [00:03<00:03, 31.17it/s]
53%|#####3 | 107/201 [00:03<00:03, 31.30it/s]
55%|#####5 | 111/201 [00:03<00:02, 31.40it/s]
57%|#####7 | 115/201 [00:03<00:02, 31.53it/s]
59%|#####9 | 119/201 [00:03<00:02, 31.74it/s]
61%|######1 | 123/201 [00:03<00:02, 31.19it/s]
63%|######3 | 127/201 [00:04<00:02, 31.27it/s]
65%|######5 | 131/201 [00:04<00:02, 31.26it/s]
67%|######7 | 135/201 [00:04<00:02, 31.30it/s]
69%|######9 | 139/201 [00:04<00:01, 31.42it/s]
71%|#######1 | 143/201 [00:04<00:01, 31.42it/s]
73%|#######3 | 147/201 [00:04<00:01, 30.99it/s]
75%|#######5 | 151/201 [00:04<00:01, 31.15it/s]
77%|#######7 | 155/201 [00:04<00:01, 31.21it/s]
79%|#######9 | 159/201 [00:05<00:01, 31.05it/s]
81%|########1 | 163/201 [00:05<00:01, 31.09it/s]
83%|########3 | 167/201 [00:05<00:01, 31.20it/s]
85%|########5 | 171/201 [00:05<00:00, 31.39it/s]
87%|########7 | 175/201 [00:05<00:00, 31.62it/s]
89%|########9 | 179/201 [00:05<00:00, 31.47it/s]
91%|#########1| 183/201 [00:05<00:00, 31.39it/s]
93%|#########3| 187/201 [00:05<00:00, 31.42it/s]
95%|#########5| 191/201 [00:06<00:00, 31.56it/s]
97%|#########7| 195/201 [00:06<00:00, 31.29it/s]
99%|#########9| 199/201 [00:06<00:00, 31.22it/s]
100%|##########| 201/201 [00:06<00:00, 31.22it/s]
Calculate. The likelihood ratio is the difference of the log likelihoods.
Visualize. This comparison shows that stimulus x tone interaction term in the model has little systematic impact on the goodness-of-fit as given by the likelihood except around 300 - 375 ms poststimulus, largest over central scalp (MiCe).
fig, ax = plt.subplots(figsize=(12, 3))
# render
im = ax.imshow(likelihood_ratio.T, interpolation="none", aspect=16)
cb = fig.colorbar(im, ax=ax)
# label
ax.set_title("Likelihood ratio")
ax.set(xlabel="Time (ms)", ylabel="Channel")
xtick_labels = range(-200, 600, 100)
ax.set_xticks([likelihood_ratio.index.get_loc(tick) for tick in xtick_labels])
ax.set_xticklabels(xtick_labels)
ax.set_yticks(range(len(likelihood_ratio.columns)))
ax.set_yticklabels(likelihood_ratio.columns)
fig.tight_layout()
fig, ax = plt.subplots(figsize=(8, 3))
ax.set_title("Likelihood Ratio")
_ = (lmg_full.llf - lmg_reduced.llf).plot(ax=ax)
Total running time of the script: ( 0 minutes 58.258 seconds)