Demo
An example dataset and a demo script are available in the demo directory.
Two .mat files (data_face, data_scrambled) contain preprocessed single-trial MEG data for one subject (two conditions). These have been obtained by reading in data corresponding to trials in which faces and scrambled stimuli were shown on the screen (90 trials for each condition). These were baselined and centered around stimulus onset (from 0.5 s prestimulus to 1 s post-stimulus), low-pass filtered at 100 Hz and resampled to 300 Hz.
For spatial/temporal feature selection, a neighbours structure has also been created, containing the sensor layout and the time axis after resampling (neighbours.mat).
These files were obtained using read_meg_data and get_sensor_info respectively.
The following steps can be tested by running the mvpa_demo script:
To prepare the data for classification, we can now create a data matrix (272 channels x 450 time points x 180 trials) and a label vector.
%after navigating to the demo directory
addpath('..')
mvpa_setup
[data, labels] = prepare_binary_data('data_face.mat','data_scrambled.mat');
load('neighbours.mat')
Next, we run time-resolved decoding on all channels using time windows of 5 samples each, selecting a shorter epoch for decoding (we will start at 100 ms before stimulus onset). We will use 5-fold CV, and we will try this with and without trial averaging (100 permutations of 3-trial groups) and multivariate noise normalization (MNN):
sens_result1 = time_resolved_kfold(data, labels, 'mnn', false, 'kfold', 5, 'weights', true, ....
'window_length', 5, 'decoding_window', [-0.1 1], 'time', neighbours(1).time);
sens_result2 = time_resolved_kfold(data, labels, 'mnn', true, 'pseudo', [3 100], 'kfold', 5, 'weights', true,...
'window_length', 5, 'decoding_window', [-0.1 1], 'time', neighbours(1).time);
if ~exist('results','dir'), mkdir ('results'); end
save(fullfile('results','sens_result.mat'), 'sens_result*')
Next, we can plot the accuracy (or any other metric) over time, making sure we get the right time axis:
%get time axis for decoding window and save indices to timepoints at -0.1 and 0.1 s
time = neighbours(1).time;
tpstart = find(round(time,3)==-0.1);
dec_time = neighbours(1).time(tpstart:end); %time axis for decoding window
time = dec_time(1:5:end); %downsampled time axis for decoding window
tp100 = find(round(time,3)==0.1); %100 ms time point
%calculate the standard error of the mean using fold-wise accuracy
sem1 = std(sens_result1.AccuracyFold,1)/sqrt(5);
sem2 = std(sens_result2.AccuracyFold,1)/sqrt(5);
%plot time-resolved accuracy for both analyses in one plot
figure('color','w')
plot_time_results(sens_result1.Accuracy, sem1,'time',time,'smooth',5, 'legend', 'No MNN/avg','ylim',[35 105]); hold on
plot_time_results(sens_result2.Accuracy, sem2,'time',time,'smooth',5, 'legend', 'MNN+avg', 'color', 'b','ylim',[35 105])
%save figure
if ~exist('figures',dir), mkdir('figures'); end
saveas(gcf, fullfile('figures','time-resolved-accuracy.png'))
Time-resolved decoding results (accuracy over time). We can clearly see the benefit of using MNN and trial averaging in terms of increasing SNR.
In the decoding analysis, we set weights to true, so the results structure contains the raw classifier weights (results.Weights), the activation patterns obtained using the data covariance matrix as per the Haufe method (results.WeightPatterns) and the min-max normalized patterns (results.WeightPatternsNorm). We can compare the average weight-based patterns with and without MNN and averaging:
%plot averaged weight-based patterns from 100 ms onwards
ave_weights1 = mean(sens_result1.WeightPatternsNorm(:,tp100:end),2);
ave_weights2 = mean(sens_result2.WeightPatternsNorm(:,tp100:end),2);
figure;
subplot(1,2,1)
plot_sensor_results(ave_weights1, neighbours, 'window_length', 1, 'time_label', false, 'colorlim', [0 1], 'colormap', 'parula'); colorbar; title('No MNN/Ave', 'FontWeight', 'normal')
subplot(1,2,2)
plot_sensor_results(ave_weights2, neighbours, 'window_length', 1, 'time_label', false, 'colorlim', [0 1], 'colormap', 'parula'); colorbar; title('MNN + Ave', 'FontWeight', 'normal')
saveas(gcf, fullfile('figures','average-weights.png'))
Normalized weight-based patterns averaged over time. The strongest weights are localized to a right occipital cluster. Since we used broadband MEG data, this suggests that we're picking up the evoked response to faces.
We can also plot these patterns on the sensor layout over time (here only for the whitened and averaged data) as a figure (weights2.png) or as a movie (weights2.avi)):
%time-resolved plots (only for MNN/averaged data)
plot_sensor_results(sens_result2.WeightPatternsNorm, neighbours, 'time', time, 'window_length', 0.1, 'colorlim', [0 1], 'colormap', 'parula')
saveas(gcf, fullfile('figures','weights2.png'))
movie_sensor_results(sens_result2.WeightPatternsNorm, neighbours, fullfile('figures','weights2.avi'), 'colorlim', [0 1], 'colormap', 'parula', 'result_type', 'Weight patterns')
Using whole-head data, we can also look at the generalizability of the decoding model across time, which can give us an idea of how transient or sustained the underlying information is. We'll run this with the default window length of 1 (a decoder at every time-point to get the maximal temporal resolution). Then, we'll plot the two time generalization matrices (with/without MNN and averaging). Just as an example, we'll highlight the time clusters with accuracies higher than 85% and containing at least 5 time-points:
tg_accuracy1 = temporal_generalization(data,labels,'mnn',false, 'decoding_window', [-0.1 1], 'time', neighbours(1).time); %temporal generalization without MNN/averaging
tg_accuracy2 = temporal_generalization(data,labels,'mnn',true,'pseudo', [3 100], 'decoding_window', [-0.1 1], 'time', neighbours(1).time); %with MNN/averaging
save(fullfile('results','tg_accuracy.mat'),'tg_accuracy*')
%we'll highlight clusters (size>=5) with accuracies over 85
mask1 = logical(tg_accuracy1>85); %create mask
mask2 = logical(tg_accuracy2>85);
colorlim = [min(min([tg_accuracy1;tg_accuracy2])) max(max([tg_accuracy1;tg_accuracy2]))]; %common color axis
figure('color','w')
subplot(1,2,1)
plot_temporal_generalization(tg_accuracy1,'time',dec_time, 'mask',mask1, 'clustersize',5,'colormap', 'parula','colorbar_label',[],'colorlim', colorlim, 'title', 'No MNN/Ave')
subplot(1,2,2)
plot_temporal_generalization(tg_accuracy2,'time',dec_time, 'mask',mask2, 'clustersize',5,'colormap', 'parula','colorlim', colorlim, 'title', 'MNN + Ave')
saveas(gcf, fullfile('figures','tg-accuracy.png'))
Temporal generalization results obtained with and without MNN and trial averaging. Accuracies are plotted against the training and test timepoints. The results show fairly sustained patterns, with a short-lived dip in accuracy after M170 latency.
Searchlight classification can help with feature reduction while maintaining whole-brain coverage (at extra computational cost). To use a different and less expensive hold-out validation approach, we will simply divide our data into a training and a test set. To do this, we prepare our SVM-ready datasets separately for two halves of the data. Next, we run the searchlight_holdout function with and without MNN and trial averaging, as above. You can't control the searchlight spatial resolution at the sensor space with this function (it's based on immediate neighbours of each sensors according to the Fieldtrip template), but you can control the time window length - we'll set this to 10 to speed things up a bit.
Finally, we will plot the results as movies (srcl_accuracy1.avi - no MNN or averaging, srcl_accuracy2.avi - MNN and averaging), and as a figure showing the average searchlight performance (100 ms onwards) for each analysis. Note that the validation scheme and the reduced number of features/increased time window length can reduce performance.
%create separate training & test datasets and use a speedier (but not exhaustive) approach
load('data_face.mat'); load('data_scrambled.mat');
[train_data, train_labels] = prepare_binary_data(data_face(:,:,1:45), data_scrambled(:,:,46:90)); %divide data into halves
[test_data, test_labels] = prepare_binary_data(data_face(:,:,46:90), data_scrambled(:,:,1:45));
%searchlight decoding with longer time windows for speed (without and with MNN/trial averaging)
[srcl_accuracy1, srcl_Fscore1] = searchlight_holdout(train_data, train_labels, test_data, test_labels, neighbours, 'window_length', 10, 'mnn', false);
[srcl_accuracy2, srcl_Fscore2] = searchlight_holdout(train_data, train_labels, test_data, test_labels, neighbours, 'window_length', 10, 'mnn', true, 'pseudo', [3 100]);
save(fullfile('results','srcl_sens_result.mat'),'srcl*')
%the time axis has been downsampled further
srcl_time = neighbours(1).time(1:10:end); %get time axis in steps of 10 samples
tp100 = find(round(srcl_time,3)==0.1); %find the 100 ms time point
%plot movies showing searchlight accuracies over time
colorlim = [min(min([srcl_accuracy1;srcl_accuracy2])) max(max([srcl_accuracy1;srcl_accuracy2]))]; %get common color axis
movie_sensor_results(srcl_accuracy1, neighbours, fullfile('figures','srcl_accuracy1.avi'), 'colormap', 'parula', 'colorlim', colorlim)
movie_sensor_results(srcl_accuracy2, neighbours, fullfile('figures','srcl_accuracy2.avi'), 'colormap', 'parula', 'colorlim', colorlim)
% plot average searchlight accuracy after 100 ms for both analyses
mean_acc1 = mean(srcl_accuracy1(:,tp100:end),2);
mean_acc2 = mean(srcl_accuracy2(:,tp100:end),2);
colorlim = [min([mean_acc1;mean_acc2]) max([mean_acc1;mean_acc2])];
figure('color','w')
subplot(1,2,1)
plot_sensor_results(mean_acc1, neighbours, 'colormap', 'parula', 'colorlim', colorlim, 'window_length', 1, 'time_label', false); colorbar; title('No MNN/Ave', 'FontWeight', 'normal')
subplot(1,2,2)
plot_sensor_results(mean_acc2, neighbours, 'colormap', 'parula', 'colorlim', colorlim, 'window_length', 1, 'time_label', false); colorbar; title('MNN + Ave', 'FontWeight', 'normal')
saveas(gcf,fullfile('figures','average-srcl-accuracy.png'))
Searchlight results averaged across time. With MNN and trial averaging, high accuracies are more focal (and coincide better with the weights obtained in whole-brain classification), but we don't see the same increase in performance as when using the entire sensor set.