Note
Go to the end to download the full example code.
Estimate the empirical confidence interval#
This example illustrates how to estimate the confidence interval
import numpy as np
from frites.simulations import sim_local_cc_ms
from frites.dataset import DatasetEphy
from frites.workflow import WfMi
from frites import set_mpl_style
import matplotlib.pyplot as plt
set_mpl_style()
Plotting functions#
First, we define the function that is then going to be used for plotting the results
def plot(mi, pv, ci, color='C0', p=0.05, title='', units='MI (bits)'):
# figure definition
n_cis, n_rois = len(ci['ci']), len(mi['roi'])
width, height = int(np.round(4 * n_rois)), int(np.round(4 * n_cis))
fig, axs = plt.subplots(
nrows=n_cis, ncols=n_rois, sharex=True, sharey=True,
figsize=(width, height))
fig.suptitle(title, fontweight='bold')
# select significant results
mi_s = mi.copy()
mi_s.data[pv.data >= p] = np.nan
# plot the results
for n_r, r in enumerate(mi['roi'].data):
for n_c, c in enumerate(ci['ci'].data):
plt.sca(axs[n_c, n_r])
plt.plot(mi['times'].data, mi.sel(roi=r).data, color='C3',
linestyle='--')
plt.plot(mi['times'].data, mi_s.sel(roi=r).data, color=color, lw=3)
plt.fill_between(
mi['times'].data, ci.sel(ci=c, roi=r, bound='high'),
ci.sel(ci=c, roi=r, bound='low'), alpha=.5, color=color)
plt.title(f"ROI={r}; CI={c}%")
plt.ylabel(units)
Data simulation#
Let’s simulate some data with 10 subjects, 100 epochs per subject and 2 brain regions. As a result, we get a variable x representing the simulated neural data coming from the 10 subjects and y, the task-related variable.
n_subjects, n_epochs, n_roi = 10, 100, 2
x, y, roi, times = sim_local_cc_ms(n_subjects, n_epochs=n_epochs, n_roi=n_roi,
random_state=0)
dt = DatasetEphy(x.copy(), y=y, roi=roi, times=times)
Empirical confidence interval with FFX models#
Then, we estimate the confidence interval when using a fixed-effect model
# computes mi
wf = WfMi(mi_type='cc', inference='ffx')
mi, pv = wf.fit(dt, n_perm=200, n_jobs=1, random_state=0)
# computes confidence interval
ci = wf.confidence_interval(dt, n_boots=200, ci=[95, 99.9], n_jobs=1,
random_state=0)
print(ci)
# plot the results
# sphinx_gallery_thumbnail_number = 1
plot(mi, pv, ci, title='CI - FFX model')
plt.show()
0%| | Estimating MI : 0/2 [00:00<?, ?it/s]
50%|█████ | Estimating MI : 1/2 [00:00<00:00, 5.30it/s]
100%|██████████| Estimating MI : 2/2 [00:00<00:00, 5.62it/s]
100%|██████████| Estimating MI : 2/2 [00:00<00:00, 5.61it/s]
0%| | Estimating CI : 0/2 [00:00<?, ?it/s]
50%|█████ | Estimating CI : 1/2 [00:00<00:00, 4.14it/s]
100%|██████████| Estimating CI : 2/2 [00:00<00:00, 4.16it/s]
100%|██████████| Estimating CI : 2/2 [00:00<00:00, 4.15it/s]
<xarray.DataArray (ci: 2, bound: 2, times: 100, roi: 2)> Size: 6kB
array([[[[-7.23073434e-04, -7.16491311e-04],
[-7.22913886e-04, -7.21564225e-04],
[-7.12451641e-04, -3.77852906e-04],
[-7.21695018e-04, -7.23084144e-04],
[-7.20416487e-04, -7.22600962e-04],
[-7.20839773e-04, -7.19997857e-04],
[-7.21475692e-04, -7.16736191e-04],
[-7.21381162e-04, -7.22534431e-04],
[-7.20533659e-04, -7.21102988e-04],
[-7.21385644e-04, -7.19823409e-04],
[-7.23042409e-04, -7.21751887e-04],
[-7.21119577e-04, -7.21950026e-04],
[-7.17653544e-04, -5.43816248e-04],
[-7.22834957e-04, -7.22015509e-04],
[-6.43601292e-04, -7.22082099e-04],
[-7.20890472e-04, -7.01435260e-04],
[-7.21521094e-04, -7.06132210e-04],
[-7.20710494e-04, -7.21843622e-04],
[-7.22542172e-04, -7.21560733e-04],
[-7.22901896e-04, -7.22991012e-04],
...
[ 4.94244881e-03, 1.07339406e-02],
[ 9.79310367e-03, 4.22542682e-03],
[ 6.94442075e-03, 8.52262974e-03],
[ 8.93999357e-03, 1.42579032e-02],
[ 6.89398730e-03, 1.20471157e-02],
[ 1.16844224e-02, 4.71777236e-03],
[ 6.02946244e-03, 5.92530705e-03],
[ 9.03900806e-03, 4.65719309e-03],
[ 1.27212107e-02, 6.39466289e-03],
[ 1.15400851e-02, 1.47285946e-02],
[ 6.34947699e-03, 4.47972724e-03],
[ 4.56949603e-03, 4.61853016e-03],
[ 6.94840448e-03, 4.63257404e-03],
[ 1.17303263e-02, 9.08010732e-03],
[ 4.93415445e-03, 1.25797410e-02],
[ 9.11222771e-03, 6.01635734e-03],
[ 6.81227539e-03, 1.18603073e-02],
[ 5.50806243e-03, 1.50310928e-02],
[ 9.75895487e-03, 1.13908155e-02],
[ 1.10290460e-02, 7.38632120e-03]]]])
Coordinates:
* ci (ci) float64 16B 95.0 99.9
* bound (bound) <U4 32B 'low' 'high'
* times (times) int64 800B 0 1 2 3 4 5 6 7 8 ... 91 92 93 94 95 96 97 98 99
* roi (roi) <U5 40B 'roi_0' 'roi_1'
Empirical confidence interval with RFX models#
When using the random-effect model, it’s either possible to estimate the confidence interval on the returned mutual-information or on t-values. To do the switch, you can use the parameter rfx_es for choosing between ‘mi’ or ‘tvalues’
0%| | Estimating MI : 0/2 [00:00<?, ?it/s]
50%|█████ | Estimating MI : 1/2 [00:00<00:00, 1.05it/s]
100%|██████████| Estimating MI : 2/2 [00:01<00:00, 1.07it/s]
100%|██████████| Estimating MI : 2/2 [00:01<00:00, 1.07it/s]
0%| | Estimating CI : 0/2 [00:00<?, ?it/s]
50%|█████ | Estimating CI : 1/2 [00:01<00:01, 1.52s/it]
100%|██████████| Estimating CI : 2/2 [00:03<00:00, 1.51s/it]
100%|██████████| Estimating CI : 2/2 [00:03<00:00, 1.51s/it]
confidence interval on t-values
wf = WfMi(mi_type='cc', inference='rfx')
_, pv = wf.fit(dt, n_perm=200, n_jobs=1, random_state=0)
tv = wf.tvalues
ci = wf.confidence_interval(dt, n_boots=200, ci=[95, 99.9], n_jobs=1,
random_state=0, rfx_es='tvalues')
plot(tv, pv, ci, title='CI - RFX model / T-values', units='T-values')
plt.show()
0%| | Estimating MI : 0/2 [00:00<?, ?it/s]
50%|█████ | Estimating MI : 1/2 [00:00<00:00, 1.06it/s]
100%|██████████| Estimating MI : 2/2 [00:01<00:00, 1.08it/s]
100%|██████████| Estimating MI : 2/2 [00:01<00:00, 1.08it/s]
0%| | Estimating CI : 0/200 [00:00<?, ?it/s]
22%|██▎ | Estimating CI : 45/200 [00:00<00:00, 2785.23it/s]
45%|████▌ | Estimating CI : 90/200 [00:00<00:00, 2784.64it/s]
68%|██████▊ | Estimating CI : 135/200 [00:00<00:00, 2786.72it/s]
90%|█████████ | Estimating CI : 180/200 [00:00<00:00, 2787.71it/s]
100%|██████████| Estimating CI : 200/200 [00:00<00:00, 2651.64it/s]
Total running time of the script: (0 minutes 10.439 seconds)
Estimated memory usage: 493 MB