Smoothed Proportion of Looks Over Time (SPLOT) method

This script is based on the paper by Belopolsky (2020) and implements the SPLOT method in Matlab. This method uses cluster based permutation testing to analyze gaze data with high temporal resolution and without risking the inflation of alpha error. The script uses files (.csv) with binary gaze data over the time course to be analyzed as input. Binary data reflect dwell times or fixations directed to a specific Area of Interest (AOI). In a first step the script smoothes the binary data by means of Gaussian kernel. In a next step data is averaged across trials. Lastly, the actual cluster based permutation testing is performed.

Contents

Loading of binary data for AOI

% Select directory with .csv files for each condition seperately
% containing binary gaze data over time
srcGazeCond1 = uigetdir('.\SPLOT\Data\', 'Select folder for files in Condition 1'); % condition 1 (e.g. negative affect)
srcGazeCond2 = uigetdir('.\SPLOT\Data\', 'Select folder for files in Condition 2'); % condition 2 (e.g. positive affect)

% Loads data from the directories
filesGazeCond1 = dir([srcGazeCond1, '\*.csv']);
filesGazeCond2 = dir([srcGazeCond2, '\*.csv']);

% Output will be saved in a subfolder of the SPLOT-folder
dest = '.\SPLOT\Out\';

% Set changeable variables
fps = 120;              % Set resolution of data
duration = 9 * fps;    % Set duration of analysis sequence
iterations = 1000;      % Set number of iterations of permutation

% Creates empty cells for the two conditions (e.g. affect types)
cond1Gaze = cell(length(filesGazeCond1),1);
cond2Gaze = cell(length(filesGazeCond2),1);

% Warning for unequal number of files to be compared
if (length(filesGazeCond1) ~= length(filesGazeCond2))
    disp('Warning: Unequal number of files in folders');
end

Smoothing of binary data

% Interates over files
for j=1:length(filesGazeCond1)
    gazeDataCond1 = readtable([srcGazeCond1,'\', filesGazeCond1(j).name]); % Loads individual file for Condition 1
    gazeDataCond2 = readtable([srcGazeCond2,'\', filesGazeCond2(j).name]); % Loads individual file for Condition 2

    % Applies gaussian kernel
    smoothedGazeCond1 = zeros(size(gazeDataCond1,1), size(gazeDataCond1,2));
    smoothedGazeCond2 = zeros(size(gazeDataCond2,1), size(gazeDataCond2,2));
    for i = 1:size(gazeDataCond1,2)
        smoothedGazeCond1(:,i) = conv(gazeDataCond1{:,i}, gausswin(98)/sum(gausswin(98)), 'same');
    end

    for i = 1:size(gazeDataCond2,2)
        smoothedGazeCond2(:,i) = conv(gazeDataCond2{:,i}, gausswin(98)/sum(gausswin(98)), 'same');
    end

    % Adds smoothed data to empty data container
    cond1Gaze{j,1} = smoothedGazeCond1;
    cond2Gaze{j,1} = smoothedGazeCond2;

end

Creating permutation cell

% Adds all data in one matrix
totalPermTrials = size(cond1Gaze{2,1},2) + size(cond2Gaze{1,1},2);
permCell = cell(size(cond1Gaze,1),totalPermTrials);
for j=1:size(cond1Gaze,1)
    for i=1:size(cond1Gaze{j,1},2)
        t = cond1Gaze{j,1};
        permCell{j,i} = t(:,i);
    end
    if size(cond1Gaze{j,1},2) == 3
        for i=(size(cond1Gaze{j,1},2)+1) + 1:(size(cond1Gaze{j,1},2)+1) + size(cond2Gaze{1,1},2)
            t = cond2Gaze{j,1};
            permCell{j,i} = t(:,i-(size(cond1Gaze{j,1},2)+1));
        end
    else
        for i=size(cond1Gaze{j,1},2) + 1:size(cond1Gaze{j,1},2) + size(cond2Gaze{1,1},2)
            t = cond2Gaze{j,1};
            permCell{j,i} = t(:,i-size(cond1Gaze{j,1},2));
        end
    end
end

Permutation iteration

% Secures original permutation cell
original = permCell;
tValues = [];

for k=1:iterations

    % Creates random permuation cell
    for j=1:size(original,1)
        t = randperm(size(original,2));
        tmpCell = original(j,:);
        for i=1:length(t)
            permCell{j,i} = tmpCell{t(i)};
        end
    end

    % Divides created permutation cell into two artifical conditions
    permCondTrials1 = permCell(:,1:floor(size(permCell,2)/2));
    permCondTrials2 = permCell(:,floor(size(permCell,2)/2+1):size(permCell,2));

    % Averages across trials
    permCond1 = averagePermTrials(permCondTrials1, duration);
    permCond2 = averagePermTrials(permCondTrials2, duration);

    % Perform paired t-test on each point in time-series
    [h,p,ci,stats] = ttest(permCond1,permCond2);


    % Defines clusters consisting of two or more adjacent significant
    % time points
    cluster = clusterValues(p, 0.05);

    % Defines cluster strengths
    [t, cl] = addTValues(cluster, stats);
    t = max(t);

    % Adds cluster strengths to distribution
    tValues = [tValues; t];
end

Defining threshold of significance ´

% Sorts t-values in distribution
tValues = sort(tValues);

% Gets the 95th percentil of distribution
tThres = prctile(tValues,95);

Preparing & analyzing comparison of original time-series

% Divides original permutation cell
orig1 = original(:,1:floor(size(original,2)/2));
orig2 = original(:,floor(size(original,2)/2+1):size(original,2));

% Averages across trials
origCond1 = averagePermTrials(orig1, duration);
origCond2 = averagePermTrials(orig2, duration);

% Compares time points in time series between conditions by means of paired
% t-tests
[h,p,ci,stats] = ttest(origCond1,origCond2);

% Computes mean cohen's d for each cluster
% Ruggero G. Bettinardi (2022). computeCohen_d(x1, x2, varargin)
%(https://www.mathworks.com/matlabcentral/fileexchange/62957-computecohen_d-x1-x2-varargin),
% MATLAB Central File Exchange. Retrieved December 20, 2022.
% Small changes have been made to the function to handle and return vectors:
% line 40: mean_x1  = mean(x1, 1);
% line 41: mean_x2  = mean(x2, 1);
% line 69: d        = meanDiff./ s;
d = computeCohen_d(origCond1, origCond2, 'paired');

% Defines clusters
cluster = clusterValues(p, 0.05);

% Defines cluster strengths
[t, cl] = addTValues(cluster, stats);

% Adds threshold to matrix
cl{end + 1,1} = tThres;
cl{end, 2} = "Threshold";

% Calculated p-values
PERCENTRANK = @(YourArray, TheProbes) reshape( mean( bsxfun(@le, YourArray(:), TheProbes(:).') ) * 100, size(TheProbes) );
percentiles = PERCENTRANK(tValues, t);
p_values = 1 - percentiles/100;

% Adds p-values to matrix
for i= 1:length(t)
    cl{i,3} = p_values(i);
end

% Adds cohen's d value to matrix
for i=1:length(t)
    cl{i,4} = mean(d(cluster{i}));
end

% Converts matrix to table and exports as .csv file
clTable = array2table(cl);
writetable(clTable, [dest, 'SPLOT_clusters.csv']); % Save table in output folder

Ploting results

% Averages across participants
avgCond1 = mean(origCond1,1);
avgCond2 = mean(origCond2,1);


% Creates data container for confidence intervals
lowConfCond1 = zeros(1, length(origCond1));
upperConfCond1 = zeros(1, length(origCond1));

lowConfCond2 = zeros(1, length(origCond2));
upperConfCond2 = zeros(1, length(origCond2));

% Defines confidence intervals
for k =1:length(origCond1)
    SEM = std(origCond1(:,k))/sqrt(length(origCond1(:,k)));         % Standard Error
    ts = tinv([0.025  0.975],length(origCond1(:,k))-1);             % T-Score
    CI = mean(origCond1(:,k)) + ts.*SEM;                            % Confidence Intervals
    lowConfCond1(1,k) = CI(1);
    upperConfCond1(1,k) = CI(2);

    SEM = std(origCond2(:,k))/sqrt(length(origCond2(:,k)));         % Standard Error
    ts = tinv([0.025  0.975],length(origCond2(:,k))-1);             % T-Score
    CI = mean(origCond2(:,k)) + ts.*SEM;                            % Confidence Intervals
    lowConfCond2(1,k) = CI(1);
    upperConfCond2(1,k) = CI(2);
end

% Transforms time series data for export & adds confidence intervals to data
avgCond1Table = array2table(avgCond1');
avgCond1Table.Properties.VariableNames{'Var1'} = 'Look_probability';
lowConfCond1Table = array2table(lowConfCond1');
lowConfCond1Table.Properties.VariableNames{'Var1'} = 'Lower_bound';
upperConfCond1Table = array2table(upperConfCond1');
upperConfCond1Table.Properties.VariableNames{'Var1'} = 'Upper_bound';
tableCond1 = [avgCond1Table,upperConfCond1Table, lowConfCond1Table];

avgCond2Table = array2table(avgCond2');
avgCond2Table.Properties.VariableNames{'Var1'} = 'Look_probability';
lowConfCond2Table = array2table(lowConfCond2');
lowConfCond2Table.Properties.VariableNames{'Var1'} = 'Lower_bound';
upperConfCond2Table = array2table(upperConfCond2');
upperConfCond2Table.Properties.VariableNames{'Var1'} = 'Upper_bound';
tableCond2 = [avgCond2Table, upperConfCond2Table, lowConfCond2Table];

% Exports data to .csv
writetable(tableCond1, [dest, 'Timeseries_negative.csv']);
writetable(tableCond2, [dest, 'Timeseries_positive.csv']);

% Defines areas of clusters and p-values
patchX = [];
txtPvalues = cell(length(t),1);

for i= 1:size(cl,1)-1
    if length(cl{i,2}) >= 10
        patchX(end + 1) = cl{i,2}(1);
        patchX(end + 1) = cl{i,2}(end);
        txtPvalues{i,1} = strcat('p = ',strrep(num2str(cl{i,3}), '0.', '.'));
    end
end

txtP= txtPvalues(~cellfun('isempty',txtPvalues));

% Sets colors and formating of figure
str1 = '#FFAE49';
str2 = '#024B7A';
color1 = sscanf(str1(2:end),'%2x%2x%2x',[1 3])/255;
color2 = sscanf(str2(2:end),'%2x%2x%2x',[1 3])/255;
figure('DefaultAxesFontSize',16)

% Plots shaded area for time of affect modulation
patch([360 360  720 720], [min(ylim) max(ylim) max(ylim) min(ylim)], [0.5 0.5 0.5], 'FaceAlpha', 0.2, 'EdgeColor','none')
hold on

% Plots time series
plot(1:duration, avgCond1(1:duration), 'LineWidth',1.5, 'Color', color2)
xlabel('Time in frames')
ylabel('Look probability')
hold on
plot(1:duration, avgCond2(1:duration), 'LineWidth',1.5, 'Color', color1)
axis([0,duration,0,1])
hold on

% Plots confidence intervals
patch([1:(duration) fliplr(1:(duration))], [lowConfCond1 fliplr(upperConfCond1)], color2, 'FaceAlpha',0.2, 'EdgeColor','none')
hold on
patch([1:(duration) fliplr(1:(duration))], [lowConfCond2 fliplr(upperConfCond2)], color1, 'FaceAlpha',0.2, 'EdgeColor','none')
hold on

% Plots thicker lines for clusters
for i= 1:2:length(patchX)
    plot(patchX(i):patchX(i+1), avgCond1(patchX(i):patchX(i+1)), 'LineWidth',4, 'Color', color2)
    hold on
    plot(patchX(i):patchX(i+1), avgCond2(patchX(i):patchX(i+1)), 'LineWidth',4, 'Color', color1)
    hold on
end

% Plots p-values of clusters
jitterY = repmat([0.05, -0.1], 1, 30);
jitterX = repmat([-20, 20], 1, 30);
counter = 1;
for i=1:length(txtP)
    text((patchX(counter)+ jitterX(i)),(avgCond1(patchX(counter))+ jitterY(i)),txtP{i})
    counter = counter +2;
    hold on
end

% Adds legend
l = legend('Affect modulation', 'Negative', 'Positive', 'Location', 'NorthEast');
l.FontSize = 12;

% Saves plot
imagewd = getframe(gcf);
plotdest = strcat(dest, 'fig.tiff');
imwrite(imagewd.cdata,plotdest);

About

Important functions

Averaging over trials

function [avgTrial] = averagePermTrials(trials, duration)
    avgTrial = zeros(size(trials,1), duration);
    divider = zeros(size(trials,1),1);
    for j=1:size(trials,1)
        for i=1:size(trials,2)
            try
                avgTrial(j,:) = avgTrial(j,:) + trials{j,i}';
                divider(j,1) = 4;
            catch
                divider(j,1) = 3;
            end
        end
        avgTrial(j,:) = avgTrial(j,:) / divider(j,1);
    end
end

% Defining clusters
function [cluster] = clusterValues(pArray, threshold)
    cluster = {};
    newCluster = 1;
    for i=1:length(pArray)
        if pArray(i) < threshold
            if newCluster
                cluster{end+1} = i;
            else
                cluster{end}(end+1) = i;
            end
            newCluster = 0;
        else
            newCluster = 1;
        end
    end
    % Removes cluster with only one entry
    for i=length(cluster):-1:1
        if length(cluster{i}) <= 1
            cluster(i) = [];
        end
    end
end

% Defining cluster strength
function [tValues, retCluster] = addTValues(cluster, stats)
tValues = zeros(length(cluster),1);
retCluster = cell(length(cluster),2);

if isempty(cluster)
    strength = max(stats.tstat);
    index = find(stats.tstat==strength);
    retCluster{end+1,1} = strength;
    retCluster{end,2} = index;
    tValues(end+1) = abs(strength);
else
    for i=1:length(cluster)
        strength = sum(stats.tstat(cluster{i}));
        retCluster{i,1} = strength;
        retCluster{i,2} = cluster{i};
        tValues(i) = abs(strength);
    end
end
end

Published with MATLAB® R2018a