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 = 10 * 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(12)/sum(gausswin(12)), 'same');
    end

    for i = 1:size(gazeDataCond2,2)
        smoothedGazeCond2(:,i) = conv(gazeDataCond2{:,i}, gausswin(12)/sum(gausswin(12)), '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);

    % 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);

% 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";

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

Ploting results

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

patchX = [];
for i= 1:size(cl,1)-1
    if cl{i} > cl{end,1}
        patchX(end + 1) = cl{i,2}(1);
        patchX(end + 1) = cl{i,2}(end);
    end
end

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

% Adds shaded background for significant clusters
if (size(patchX,2) == 2)
    patch([patchX(1) patchX(1) patchX(2) patchX(2)], [min(ylim) max(ylim) max(ylim) min(ylim)], [0.8 0.8 0.8], 'FaceAlpha', 0.5)
elseif (size(patchX,2) == 4)
    patch([patchX(1) patchX(1) patchX(2) patchX(2)], [min(ylim) max(ylim) max(ylim) min(ylim)], [0.8 0.8 0.8], 'FaceAlpha', 0.5)
    hold on
    patch([patchX(3) patchX(3) patchX(4) patchX(4)], [min(ylim) max(ylim) max(ylim) min(ylim)], [0.8 0.8 0.8], 'FaceAlpha', 0.5)
end
hold on
legend('Condition 1', 'Condition 2', 'Location', 'NorthEast')

% 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']);

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® R2022a