Exercise 1 (4 points):

Load the file Ch2-EEG-4.mat into MATLAB (Kramer & Eden, Chapter 2). You’ll find twovariables:

EEG = the EEG data, consisting of 1,000 trials, each observed for 1 s;
t = the time axis, which ranges from 0 s to 1 s.

These data have a similar structure to the data studied in chapter 2. To collect these data, a stimulus was presented to the subject at 0.25 s. Analyze these data for the presence of an evoked response. To do so, answer the following questions.

Questions:

What is the sampling rate (samples/s)? Is the sampling interval the same across the recording interval? Plot all the individual trials of these data using the function "imagesc.m" (similarly to Fig. 2.4, Chapter 2). Explain what you observe in words. From your visual inspection, do you expect to find an ERP in these data?
Compute the ERP for these data, and plot the results. Do you observe an ERP (i.e., times at which the 95% confidence intervals do not include zero)? Include 95% confidence intervals in your ERP plot, and label the axes. Use the nonparametric bootstrap (Chapter 2, pages 34 – 37). Explain in a few sentences the results of your analysis, as you would to a collaborator who collected these data.
Compute the ensemble variance for these data, and plot the results. Show in the same plot the corresponding 95% confidence interval obtained via nonparametric bootstrap.
Is the process mean stationary? Is the process variance stationary? Explain what you see.



In [1]:

    
# Pkg.add("MAT")
# Pkg.add("Seaborn")
# Pkg.add("StatsBase")
# Pkg.checkout("MAT") - if issues, checkout master branch



In [2]:

    
using MAT
# using Seaborn
using PlotlyJS
import PyPlot.imshow
using StatsBase









    






 
 Plotly javascript loaded.
 To load again call 
init_notebook(true)



In [3]:

    
#import data
mat_vars = matread("Ch2-EEG-4.mat")
eeg = mat_vars["EEG"]
time = mat_vars["t"]
ntrials = size(eeg, 1)

println("Loaded EEG signal with dimensions: ", size(eeg))
println("Number of trials: ", ntrials)
println("Loaded 1 second of time with samples: ", size(time))









    



Loaded EEG signal with dimensions: (1000, 500)
Number of trials: 1000
Loaded 1 second of time with samples: (1, 500)

Question 1

What is the sampling rate (samples/s)?

500 samples per second

Is the sampling interval the same across the recording interval?

Yes.

Plot all the individual trials of these data using the function "imagesc.m" (similarly to Fig. 2.4, Chapter 2). Explain what you observe in words. From your visual inspection, do you expect to find an ERP in these data?

See figure below. The signal seems stationary except for the betweeb 150~220 samples where rate of change/frequency of the signal is very different



In [4]:

    
#plot data as an image using the Images.jl package
imshow(eeg); #if using colormap use additional input cmap = "jet"

Question 2

Compute the ERP for these data, and plot the results. Do you observe an ERP (i.e., times at which the 95% confidence intervals do not include zero)? Include 95% confidence intervals in your ERP plot, and label the axes. Use the nonparametric bootstrap (Chapter 2, pages 34 – 37). Explain in a few sentences the results of your analysis, as you would to a collaborator who collected these data.

For the plot and computations see code below.

The mean appears stationary until the stimulus was presented. When the stimulus is presented there appears a clear activity evoked by the stimulus.



In [19]:

    
function bootstrapped_erp_sample(eeg)
    #bootstraped eeg samples
    ntrials = size(eeg, 1)
    eeg_i = zeros(Int, ntrials)
    sample!(1:ntrials, eeg_i; replace=true, ordered=false)
    return vec(mean(eeg[eeg_i,:], 1))    
end

function plot_erp(eeg_mat, t; ci="sd", nboot=1000)
    
    t = vec(t) #as column vector
    nsamples = length(t)
    
    #stats - as column vectors    
    eeg_mean = vec(mean(eeg_mat, 1))
    eeg_sd = vec(std(eeg_mat, 1))
    
    if ci=="sd"
        
        #return 95% confidence interval as given by the Central Limit Theorem ( 2 STD from the mean)
        eeg_sd_mean = eeg_sd./sqrt(ntrials)

        #traces
        ci_bottom_trace = PlotlyJS.scatter(; x = t, y= eeg_mean - 2*eeg_sd_mean.*ones(nsamples), marker_color="blue", line_width=0.2, name="Lower CI")
        mean_trace = PlotlyJS.scatter(;x = t, y=eeg_mean, marker_color="blue", name="Mean ERP")
        ci_top_trace = PlotlyJS.scatter(; x = t, y= eeg_mean + 2*eeg_sd_mean.*ones(nsamples), marker_color="blue", line_width=0.2, name="Upper CI")

        data = [ci_bottom_trace, mean_trace, ci_top_trace]
        layout = Layout(;title="ERP with CLT confidence intervals",
                         legend=attr(family="Arial, sans-serif", size=20,
                                     color="grey"))
        #plot
        PlotlyJS.plot(data, layout)
    else
        assert(ci <= 100 && ci >=0)
        #for ci in [0,100] plot bootstrap confidence intervals
        ci_float  = Float64(100-ci)/200.0
        erp_boot = zeros(Float64, nboot, nsamples)
        
        for b=1:nboot
            erp_boot[b, :] = bootstrapped_erp_sample(eeg)
        end
        
        erp_boot_sorted = sort(erp_boot, 1)
        ci_side = (100-ci)/200
        ci_lower_idx = nboot * ci_side 
        ci_lower = vec(erp_boot_sorted[Int(round(ci_lower_idx)),:])
        ci_upper_idx = nboot * (1.0 - ci_side)
        ci_upper = vec(erp_boot_sorted[Int(round(ci_upper_idx)),:])
        
        
        #traces
        ci_bottom_trace = PlotlyJS.scatter(; x = t, y= ci_lower, marker_color="blue", line_width=0.2, name="Lower CI")
        mean_trace = PlotlyJS.scatter(;x = t, y=eeg_mean, marker_color="blue", name="Mean ERP")
        ci_top_trace = PlotlyJS.scatter(; x = t, y= ci_upper, marker_color="blue", line_width=0.2, name="Upper CI")

        data = [ci_bottom_trace, mean_trace, ci_top_trace]
        layout = Layout(;title="ERP with with Bootstrap Confidence Intervals",
                         legend=attr(family="Arial, sans-serif", size=20,
                                     color="grey"))
        #plot
        PlotlyJS.plot(data, layout)
    end
end









    Out[19]:





plot_erp (generic function with 1 method)



In [6]:

    
plot_erp(eeg, time, ci="sd")









    Out[6]:



In [7]:

    
plot_erp(eeg, time, ci=95, nboot=1000)









    Out[7]:



In [8]:

    
# Note: Instead plotting the ERP manually, Seaborn library could be used
# The defaults are not as nice looking as PlotlyJS, but it is often convinient

# plot the time series - can change err_style: {ci_band, ci_bars, boot_traces, boot_kde, unit_traces, unit_points}
# Seaborn.tsplot(data=eeg, time=time, err_style="ci_band", ci=95)

Question 3:

Compute the ensemble variance for these data, and plot the results. Show in the same plot the corresponding 95% confidence interval obtained via nonparametric bootstrap.



In [18]:

    
function bootstrap_ensamble_variance_sample(eeg)
    #bootstraped eeg samples
    ntrials = size(eeg, 1)
    eeg_i = zeros(Int, ntrials)
    sample!(1:ntrials, eeg_i; replace=true, ordered=false)
    return vec(var(eeg[eeg_i,:], 1))    
end

function plot_ensemble_var(eeg_mat, t; ci=95, nboot=1000)
    
    t = vec(t) #as column vector
    nsamples = length(t)
    
    #stats - as column vectors    
    eeg_var = vec(var(eeg_mat, 1))
    
    assert(ci <= 100 && ci >=0)
    #for ci in [0,100] plot bootstrap confidence intervals
    ci_float  = Float64(100-ci)/200.0
    var_boot = zeros(Float64, nboot, nsamples)

    for b=1:nboot
        var_boot[b, :] = bootstrap_ensamble_variance_sample(eeg)
    end

    var_boot_sorted = sort(var_boot, 1)
    ci_side = (100-ci)/200
    ci_lower_idx = nboot * ci_side 
    ci_lower = vec(var_boot_sorted[Int(round(ci_lower_idx)),:])
    ci_upper_idx = nboot * (1.0 - ci_side)
    ci_upper = vec(var_boot_sorted[Int(round(ci_upper_idx)),:])


    #traces
    ci_bottom_trace = PlotlyJS.scatter(; x = t, y= ci_lower, marker_color="blue", line_width=0.2, name="Lower CI")
    var_trace = PlotlyJS.scatter(;x = t, y=eeg_var, marker_color="blue", name="Ensemble Variance")
    ci_top_trace = PlotlyJS.scatter(; x = t, y= ci_upper, marker_color="blue", line_width=0.2, name="Upper CI")

    data = [ci_bottom_trace, var_trace, ci_top_trace]
    layout = Layout(;title="Ensemble Variance with Bootstrap Confidence Intervals",
                     legend=attr(family="Arial, sans-serif", size=20,
                                 color="grey"))
    #plot
    PlotlyJS.plot(data, layout)
    
end

plot_ensemble_var(eeg, time, ci=95, nboot=1000)









    Out[18]: