init_unfold
Starting unfold toolbox. Adding subfolders and other toolboxes to path... Done.

Introduction

Subjects are exploring a face stimulus. Using an eyetracker and the EEG-Eye toolbox we add markers of saccade onsets into our EEG-Data. The data were cleaned for eye-movement related artefacts using ICA.
In this tutorial we try to disentangle stimulus related ERPs and (micro)saccade related ERPs
Find the data at: https://osf.io/wbz7x/
EEG = pop_loadset('C:\Users\behinger\Desktop\face_saccades_opendata_fig10.set')
pop_loadset(): loading file C:\Users\behinger\Desktop\face_saccades_opendata_fig10.set ...
EEG = struct with fields:
setname: 'face_saccades_opendata' filename: 'face_saccades_opendata_fig10.set' filepath: 'C:\Users\behinger\Desktop' subject: '' group: '' condition: '' session: [] comments: 'Data to re-create Figure 10 in Ehinger & Dimigen, 2018, BioRxiv. Please do not use this data for other purposes without asking first (olaf.dimigen@hu-berlin.de)' nbchan: 40 trials: 1 pnts: 2070402 srate: 500 xmin: 0 xmax: 4.1408e+03 times: [1×2070402 double] data: [40×2070402 single] icaact: [] icawinv: [] icasphere: [] icaweights: [] icachansind: [] chanlocs: [1×40 struct] urchanlocs: [] chaninfo: [1×1 struct] ref: 'averef' event: [1×7661 struct] urevent: [1×7661 struct] eventdescription: {'' '' '' '' '' '' '' '' '' '' '' '' '' ''} epoch: [] epochdescription: {} reject: [1×1 struct] stats: [1×1 struct] specdata: [] specicaact: [] splinefile: '' icasplinefile: '' dipfit: [] history: '↵EEG.setname='MSAFF_01';↵EEG = eeg_checkset( EEG );↵EEG = eeg_checkset( EEG );↵EEG.etc.eeglabvers = '14.1.2'; % this tracks which version of EEGLAB is being used, you may ignore it↵EEG = pop_loadset('filename','face_saccades_opendata.set','filepath','M:\\Dropbox\\shared folders_main\\deconv\\figures\\toolbox_paper\\redo_figure10_to_match_opendata\\');↵EEG = eeg_checkset( EEG );↵pop_eegplot( EEG, 1, 1, 1);↵EEG = pop_select( EEG,'nochannel',{'L-Dia-X-(px)' 'L-Dia-Y-(px)' 'R-Dia-X-(px)' 'R-Dia-Y-(px)'});↵EEG = eeg_checkset( EEG );↵EEG = pop_saveset( EEG, 'filename','face_saccades_opendata_fig10.set','filepath','M:\\Dropbox\\shared folders_main\\deconv\\figures\\toolbox_paper\\redo_figure10_to_match_opendata\\');↵EEG = eeg_checkset( EEG );↵EEG = pop_select( EEG,'nochannel',{'L-GAZE-X-(px)' 'L-GAZE-Y-(px)' 'R-GAZE-X-(px)' 'R-GAZE-Y-(px)'});↵EEG = eeg_checkset( EEG );↵EEG.etc.eeglabvers = '14.1.1'; % this tracks which version of EEGLAB is being used, you may ignore it' saved: 'justloaded' etc: [1×1 struct] datfile: []

Defining the design

cfgDesign = [];
cfgDesign.eventtypes = {{'saccade'},{'stimonset'}}; % we model the saccadeonset and the stimulus
cfgDesign.formula = {'y ~ 1 + spl(sac_amplitude,10)','y~1'}; % 10 splines are generated; One needs to be careful to not overfit the data, regularization or cross validation can help here.
EEG = uf_designmat(EEG,cfgDesign);
Multiple events with separate model-formula detected Modeling 2871 event(s) of [saccade] using formula: y~1 The spline that got removed due to collinearity in the basis set (as intended) for the effect sac_amplitude has its peak at 0.969419 This does not mean that the event-intercept represents this value! Modeling 479 event(s) of [stimonset] using formula: y~1
The term spl(sac_amplitude,10) creates 10 splines over our predictor.
We can visualise the set of splines using:
spl = EEG.unfold.splines{1};
subplot(2,1,1),histogram(spl.paramValues,100),title('Histogram of saccadic amplitudes')
splineSet = spl.splinefunction(linspace(0,5,500),spl.knots);
subplot(2,1,2),plot(linspace(0,5,500),splineSet),xlabel('saccadic amplitude')
The default setting is to put the splines at points of the quantiles of the variable. This ensures that in regions with a lot of data, lots of splines are generated and in regions with sparse data, broad splines are used. One could also put the splines in a linear or logarithmical way or define the knotsequence manually.
To understand better what basis function do, it might be helpful to have a look at the sorted designmatrix:
uf_plotDesignmat(EEG,'sort',1)
The sorted designmatrix shows overlapping predictors, these are the splines for saccade amplitude. There is one predictor removed (between sac_amp 1.5 & 2.9). This is because we have an intercept in the model and else the model would be overcomplete. Not that the empty events (all grey rows) are events that are not modeled.
cfgTimeexpand = [];
cfgTimeexpand.timelimits = [-.3,.6];
EEG = uf_timeexpandDesignmat(EEG,cfgTimeexpand);
uf_timeexpandDesignmat(): Timeexpanding the designmatrix... ...done
winrej = uf_continuousArtifactDetect(EEG,'amplitudeThreshold',250);
WARNING: boundary events were not found. So, basicrap.m will be applied over the full range of data. Please wait. This might take several seconds... 58 segments were marked. A total of 133034 samples (6.43 percent of data) was marked as bad.
EEG = uf_continuousArtifactExclude(EEG,struct('winrej',winrej)); % we find data using a simple threshold tool. More complex algorithms or manual cleaning is recommended!
Portions of data removed split up by each event in EEG.event Type: S 99 modelled eventtime: 863.10s rejected eventtime: 13.96s percent removed: 1.6% Type: S100 modelled eventtime: 0.90s rejected eventtime: 0.90s percent removed: 100.0% Type: S200 modelled eventtime: 0.90s rejected eventtime: 0.00s percent removed: 0.0% Type: buttonpress modelled eventtime: 431.10s rejected eventtime: 5.83s percent removed: 1.4% Type: fixation modelled eventtime: 766.25s rejected eventtime: 9.78s percent removed: 1.3% Type: saccade modelled eventtime: 765.77s rejected eventtime: 9.81s percent removed: 1.3% Type: stimonset modelled eventtime: 431.10s rejected eventtime: 7.29s percent removed: 1.7% Removing 6.43% of rows from design matrix (filling them with zeros)
EEG= uf_glmfit(EEG,'channel',32);
uf_glmfit(): Fitting deconvolution model...solving the equation system solving electrode 32 (of 1 electrodes in total)... 75 iterations, took 2.9s LMfit finished

Plot the results

Currently the spline-predictor is still separated in intercept + 9 spline-betas. The relationship for e.g. the P100 can be visualized. We start with a simple investigation of the raw betas, without undoing the spline-basisset
ufresult= uf_condense(EEG);
sac_amp_ix = cellfun(@(x)~isempty(x),(strfind({ufresult.param.name},'sac_amplitude'))); % get all betas related to sac_amplitude
timeix = get_min(0.115,ufresult.times); % somewhere around the p100
splinevalue = [ufresult.param(sac_amp_ix).value];
figure;
plot(splinevalue,squeeze(ufresult.beta(32,timeix,sac_amp_ix)),'-o'),xlabel('saccade amplitude')
This plot is still lacking the intercept (a DC-offset) and of course, these betas represent only weights to the above plotted set of basisfunctions.
In order to get a better picture, let's multiply the betas with the basis function to get an estimate in the predictor-value-domain.
ufresult= uf_predictContinuous(ufresult,'predictAt',{{'sac_amplitude',linspace(0.3,5,100)}});
sac_amp_ix = cellfun(@(x)~isempty(x),(strfind({ufresult.param.name},'sac_amplitude'))); % get all betas related to sac_amplitude
y = ufresult.beta(32,timeix,sac_amp_ix);
% we could add the intercept; it is a constant offset
%y = y + ufresult.beta(1,timeix,1);
plot([ufresult.param(sac_amp_ix).value],squeeze(y)),xlabel('saccade amplitude'),title(sprintf('Non-linear relationship between ERP and saccadic amplitude @ %.2f s',ufresult.times(timeix)))
It is clear from this plot, that there is a strongly non-linear relationship between saccade amplitude and ERP at around 100ms.
It's important to remember here, that above 3° saccade amplitude only very few saccades remain. It is therefore important to not overinterpret this downward trend!
To really drive home the point we add the weighted splines to the plot:
ufresult= uf_condense(EEG);
sacX =linspace(0.3,5,500); % at which points to evaluate?
splineSet = spl.splinefunction(sacX,spl.knots);
beta = squeeze(ufresult.beta(32,timeix,1:end-1)); % end-1 because of the stimulus event predictor
splineEvaluated = bsxfun(@times,splineSet,beta'); % weight the basis functions by the betas, but don't add them
plot(sacX,splineEvaluated); hold all;
plot(sacX,sum(splineEvaluated,2),'k','LineWidth',1.5), % add the basis functions at each saccade-amplitude to get the modelfit
xlabel('saccade amplitude'),title(sprintf('Non-linear relationship between ERP and saccadic amplitude @ %.2f s',ufresult.times(timeix)))
As one can also see, the modelfit is a bit rough at the edges. this is to be expected because there is not much data for very small and very large splines.
% The default parameters plot the predictor at its quantiles
ufpredict = uf_predictContinuous(ufresult); % evaluate the splines at the quantiles
uf_plotParam(ufpredict,'channel',32,'add_average',1);
plotting all parameters
Important notes to this plot:
In order to get a more representative picture, we can utilize the 'predictAt' option in the plot and put a linear spacing excluding the extremata
ufpredict= uf_predictContinuous(ufresult,'predictAt',{{'sac_amplitude',linspace(0.5,3,10)}});
ufmarginal = uf_addmarginal(ufpredict);
uf_addmarginal: working on field beta re-running uf_condense to recover unconverted splines
ufmarginal.effect_only = ufpredict.beta; % we are adding the predictor without the added intercept to show the difference
uf_plotParam(ufmarginal,'channel',32,'plotParam','sac_amplitude');
plotting selected parameters
In the top row, the output of uf_addmarginal shows what happens after we add the mean value of all other predictors (in this case only the intercept). The lower plot can be seen as 'difference' waves