NeuroAnalyzer tutorials: Statistics (1)

Load data:

using NeuroAnalyzer
using NeuroStats
using Plots
eeg = load("files/eeg.hdf");
e10 = epoch(eeg, ep_len=10);
delete_epoch!(e10, ep=1:10);
[ Info: Loaded: EEG (24 × 308480 × 1; 1204.996 s)

(!) Most statistical functions used by the NeuroAnalyzer are in the separate package: NeuroStats.jl.

Calculate and plot CI using bootstrapping

Get the first 2 seconds of the channel 1 of epochs 1:10 and calculate its 95%CI:

s = e10.data[1, 1:512, 1:10]
t = e10.epoch_time[1:512]
s_avg, s_l, s_u = bootstrap_ci(s, cl=0.95)
(s_avg = [11.654898139351815, 10.262248215819774, 11.97869012655289, 12.630062992707838, 11.570191509089227, 11.540959091716795, 10.541115219980027, 10.531517956966704, 10.421700762819713, 10.740942023467158  …  -14.93467142418011, -17.482067855648076, -16.37788284377348, -12.132211640586203, -9.671487207031488, -12.30071929298087, -15.216810566070887, -15.194893417448966, -12.842532095529629, -11.009415842579445],
 s_ci_l = [10.822860742097781, 9.408763518405909, 11.092509420219539, 11.649015715837571, 10.657671056665913, 10.635817998559563, 9.682447968681974, 9.764797992741789, 9.582581624381268, 9.847133243994605  …  -17.115832768623285, -20.197938065265543, -19.033876224269644, -14.177775578676833, -11.463232853913574, -14.505658567550112, -17.938809501531136, -17.867950175936876, -14.920114352300565, -12.864605422159261],
 s_ci_h = [12.512148422280918, 11.109775921957953, 12.871125723524765, 13.63360775095151, 12.480159298981198, 12.455200730221277, 11.419837150051526, 11.32791623654601, 11.311718997478932, 11.67202688837758  …  -12.866814069530454, -14.88802574341585, -13.869626726889976, -10.211028221997651, -7.9874481181035675, -10.231731430775788, -12.669568623887448, -12.68246561037163, -10.854067971247654, -9.236188498148802],)

Plot CI:

plot_ci(s_avg, s_l, s_u, t)

Calculate statistic using bootstrapping

Any statistical function may be calculated using bootstrapping method. The function must be provided as f="function_name(obj)", obj will be replaced with the signal data, e.g. f="mean(obj)".

s = e10.data[1, :, 1:10]
s_stat = mean(s)
s_dist  = bootstrap_stat(s, f="mean(obj)")
3000-element Vector{Float64}:
  0.20066477399049693
  0.14394971467457118
  0.07324944491156433
  0.03000076203671398
  0.1515606187996063
  0.0641443042164969
 -0.11203781864159464
  0.13498272826955598
  0.10262449640463948
  0.05463395196986447
  0.03215512833206482
  0.08024724075223144
  0.13264935639988237
  ⋮
  0.026014188132643133
  0.17937168246312005
  0.07191929213262852
  0.2452989615234543
  0.2514382077903489
  0.001506255250420807
 -0.13471432149644222
  0.2210818992910987
  0.11510522773132888
  0.021775359510878545
 -0.04079052021458942
 -0.0016443799952444494
plot_histogram(s_dist, s_stat, draw_median=false, draw_mean=true)
[ Info: Proportion of values > 0.09340081426013924: 0.529
[ Info: Proportion of values < 0.09340081426013924: 0.471