analysisPsychometric.m

%% Script for all the psychometric analysis
%
%    If this code is used in a publication, please cite the manuscript:
%    "H Huang, KN Kay, NM Gregg, G Ojeda Valencia, M In, C Kapeller, Y Shu, GA Worrell, KJ Miller, and D Hermes.
%    Single pulse electrical stimulation in white matter modulates iEEG visual responses in human early visual cortex. (Under Review)"
%
%    A preprint is available currently at doi: https://doi.org/10.1101/2025.05.05.652264.
%
%    The dataset corresponding to this code and manuscript is in BIDS format (version 1.10.0) on OpenNeuro (ds006485),
%    and it will be made publicly available upon manuscript acceptance.
%
%    Copyright (C) 2025 Harvey Huang
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
set(0, 'DefaultFigureRenderer', 'painters');

isiCols = [0, 0, 0; 0.9290 0.6940 0.1250; 0.8500 0.3250 0.0980; 0 0.4470 0.7410]; % matches the rest of the analyses. here we flipped order so yellow (first) = 0ms

%sub = '1';
sub = '2';

saveOutputs = false;

switch upper(sub)
    case '2'
        runs = 1:3; % run 3 has the ROC stim site (1 run only) (but sham)
        imgType = 'Coh';
        stimsites = {'ROC6-ROC7', 'RMO3-RMO4'};
    case '1'
        runs = 1:4;
        imgType = 'Coh';
        stimsites = {'LOC4-LOC5', 'LG6-LG7'};
    otherwise
        error('Error: choose subjects 1 or 2');
end
outdir = fullfile('output', sprintf('sub-%s', sub), 'psychometric');
mkdir(outdir);

% load the event annotations for button press information
events = [];
for ii = 1:length(runs)
    fprintf('Loading run %d...\n', runs(ii));

    eventsCurr = ieeg_readtableRmHyphens(fullfile('data', 'derivatives', 'event_annotations', sprintf('sub-%s', sub), sprintf('sub-%s_eventsCcepVisualRun%02d.tsv', sub, runs(ii))), 'electrical_stimulation_site', 1);
    eventsCurr.run = repmat(ii, height(eventsCurr), 1);
    
    events = [events; eventsCurr];
    
    clear eventsCurr;
end

% make trial types easier to sort on
events.trial_type = strrep(events.trial_type, 'high', 'zhigh');

% keep a column with the 0-coh separate, for accuracy analysis later
events.trial_type_sep = events.trial_type;

% combine 0-coh together
events.trial_type(contains(events.trial_type, 'zero')) = {'0Coh'};


% correctness vector: nan = no response (in time), 0 = incorrect, 1 = correct, -1 = zero coherence
events.respCorrect = nan(height(events), 1);
events.respCorrect(contains(events.trial_type, 'img1') & events.respButton == 3 | contains(events.trial_type, 'img2') & events.respButton == 4) = 1;
events.respCorrect(contains(events.trial_type, 'img1') & events.respButton == 4 | contains(events.trial_type, 'img2') & events.respButton == 3) = 0;
events.respCorrect(strcmp(events.trial_type, '0Coh') & ~isnan(events.respButton)) = -1;
events.respCorrectOrig = nan(height(events), 1); % keep track with the 0-coh conditions accuracies
events.respCorrectOrig(contains(events.trial_type_sep, 'img1') & events.respButton == 3 | contains(events.trial_type_sep, 'img2') & events.respButton == 4) = 1;
events.respCorrectOrig(contains(events.trial_type_sep, 'img1') & events.respButton == 4 | contains(events.trial_type_sep, 'img2') & events.respButton == 3) = 0;

% get indices with a tolerance around the desired E2V conditions. Also add e2vType column to events
events.e2vType = nan(height(events), 1);
srate = 4800;
e2vs = struct();
e2vProg = [nan, 0, 0.1, 0.2]; % Programmed E2V durations, in s
e2vs(1).mode = nan; e2vs(1).idx = isnan(events.e2v); events.e2vType(e2vs(1).idx) = 1; % the sham condition. Assign as type 1
for ii = 2:4
    [e2vs(ii).mode, e2vs(ii).idx] = getE2vs(e2vProg(ii), events.e2v/srate, 0.05);
    events.e2vType(e2vs(ii).idx) = ii; % assign to events table as corresponding e2v index
end

% image categories, in same format as the e2v struct
imgs = struct();
imglabs = cellfun(@(x) sprintf('%s%s', x, imgType), ...
    {'0', 'img1-low', 'img1-med', 'img1-zhigh', 'img2-low', 'img2-med', 'img2-zhigh'}, 'UniformOutput', false);
for ii = 1:7
    imgs(ii).label = imglabs{ii};
    imgs(ii).idx = contains(events.trial_type, imglabs{ii});
end
clear imglabs;

%% 1) Summarize and plot response times for each image category

% Colormap, markers, line style matching the 9 stimulus conditions
cmap = [0.5*ones(1, 3);
        1, 0.8, 0.8; % img1 low
        1, 0.5, 0.5; % img1 med
        1, 0, 0; % img1 high
        0.8, 0.8, 1; % img2 low
        0.5, 0.5, 1; % img2 med
        0, 0, 1]; % img2 high

[grp, ix] = sort(events.trial_type);
rt = events.respTime(ix);
                
% 1) raw response times, ignoring accuracy info, run, etc
figure('Position', [200, 200, 600, 800]); hold on
% no ID trial has stim code 1, so that's skipped
distributionPlot(rt, 'groups', grp, 'color', num2cell(cmap, 2), 'showMM', 2);
hold on
plotSpread(rt, 'distributionIdx', grp, 'spreadWidth', 1, 'distributionColors', [0.2, 0.2, 0.2]);
ylim([0, 2]);
xtickangle(30);
title('All response times regardless of accuracy, all runs')
xlabel('Stim Condition (Img-Coh)'); ylabel('Response Time (s)');
hold off
if saveOutputs
    saveas(gcf, fullfile(outdir, 'responsetimes_raw'), 'png');
    saveas(gcf, fullfile(outdir, 'responsetimes_raw'), 'svg');
end

% for raw response times calculated ANOVA and then posthoc Dunnett's test
figure
Tsimple = events(~isnan(events.respTime) & ~isnan(events.e2vType), {'respTime', 'trial_type'}); % exclude the few non-e2vType ones to be consistent w table 1
Tsimple.trial_type = categorical(Tsimple.trial_type);
[p, ~, stats] = anova1(Tsimple.respTime, Tsimple.trial_type, 'off');
[results, ~, ~, gnames] = multcompare(stats, 'CriticalValueType', 'dunnett');
tbl = array2table(results,"VariableNames", ...
    ["Group","Control Group","Lower Limit","Difference","Upper Limit","P-value"]);
tbl.("Group") = gnames(tbl.("Group"));
disp(tbl);

% 2) Response times of correct presses only. Keep 0-coh in on the left for comparison
respCorrect = events.respCorrect; % sort also by the trial_type groups
respCorrect = respCorrect(ix);
figure('Position', [200, 200, 600, 800]); hold on
% no ID trial has stim code 1, so that's skipped
distributionPlot(rt(respCorrect == 1 | respCorrect == -1), 'groups', grp(respCorrect == 1 | respCorrect == -1), 'color', num2cell(cmap, 2), 'showMM', 2);
hold on
plotSpread(rt(respCorrect == 1 | respCorrect == -1), 'distributionIdx', grp(respCorrect == 1 | respCorrect == -1), 'spreadWidth', 1, 'distributionColors', [0.2, 0.2, 0.2]);
ylim([0, 2]);
xtickangle(30);
title('RT correct presses')
xlabel('Stim Condition (Img-Coh)'); ylabel('Response Time (s)');
hold off
if saveOutputs, saveas(gcf, fullfile(outdir, 'responsetimes_correctonly.png')); end

%% 1b) Breakdown violins by stim condition

% logical ordering of stim sites by subject (sham, main (200-100-0), control). True for both subjects
stim_cond_order = 7:-1:1;

% matching the sham, 200-100-0-200-100-0 order
stim_cond_cols = [0, 0, 0;
           0 0.4470 0.7410;
           0.8500 0.3250 0.0980;
           0.9290 0.6940 0.1250;
           0 0.4470 0.7410;
           0.8500 0.3250 0.0980;
           0.9290 0.6940 0.1250];

% stim site_e2v combo
stim_cond = cell(height(events), 1);
for ii = 1:height(events)
    if isnan(events.e2vType(ii))
        stim_cond{ii} = ''; % stim mistimed events are empty
    else
        stim_cond{ii} = sprintf('%s_%d', events.electrical_stimulation_site{ii}, events.e2vType(ii));
    end
end

% iterate through each trial-type (image)
trial_types = unique(events.trial_type);
grps_all = [];
for ii = 1:length(trial_types)
    trial_type_name = trial_types{ii};
    ix_trial_type = strcmpi(trial_type_name, events.trial_type);

    % response times and stim conds for this image type
    rt_curr = events.respTime(ix_trial_type);
    stim_cond_curr = stim_cond(ix_trial_type);

    grps_curr = groupby(stim_cond_curr); % get indices of each unique stim cond for this image type
    for jj = 1:size(grps_curr, 1) % store
        grps_curr{jj, 3} = rt_curr(grps_curr{jj, 2});
    end

    % delete the empty group from mistimed stim
    grps_curr(cellfun(@isempty, (grps_curr(:, 1))), :) = []; 
    assert(size(grps_curr, 1) == 7, 'Error: not exactly 7 stim conditions');

    grps_curr = grps_curr(stim_cond_order, :);
    grps_curr(:, 4) = {trial_type_name};

    grps_all = [grps_all; grps_curr];

end

% Plot compound jitter plot

% set space between image types
img_type_space = 4;
x_pos = 1:size(grps_all, 1);
for ii = 1:length(x_pos)
    x_pos(ii) = x_pos(ii) + floor((ii-1) / 7)*img_type_space;
end
x_pos = x_pos + 1; % right-shift a bit

% make image 1 upward triangles, image 2 downward triangles
figure;
h = jitterplot(grps_all(:, 3), stim_cond_cols, 0.1);
for ii = 1:length(h)
    set(h(ii), 'XData', h(ii).XData - ii + x_pos(ii)); % switch the center value
    if any(mod(ii, 7) == [2, 3, 4])
        set(h(ii), 'Marker', '^');
    elseif any(mod(ii, 7) == [5, 6, 0])
        set(h(ii), 'Marker', 'v');
    else
        set(h(ii), 'Marker', 'o');
    end
end
set(h, 'MarkerSize', 4);
set(h, 'MarkerEdgeColor', 'none');
hold on
plot(x_pos, cellfun(@(x) mean(x, 'omitnan'), grps_all(:, 3)), 'r_', 'LineWidth', 2);
xlim([0, x_pos(end) + 1]);
ylim([0, 2]);
if saveOutputs
    saveas(gcf, fullfile(outdir, 'responsetimes_raw_stimcond'), 'png');
    saveas(gcf, fullfile(outdir, 'responsetimes_raw_stimcond'), 'svg');
end


%% 2) Summary and plot accuracy for each image category

% accuracy for each image condition
hitsMissesImg = zeros(3, 7); % first row hits, second row misses, third row binomial p
hitsMissesRuns = zeros(2*length(runs), 7); % hits/misses per run
for ii = 1:7

    % for zero-coh condition we use the original img1 img2 labels to judge hit/miss
    if ii == 1
        hits = sum(events.respCorrectOrig == 1 & imgs(ii).idx);
        misses = sum(events.respCorrectOrig == 0 & imgs(ii).idx);
    else % otherwise we use the 1s and 0s
        hits = sum(events.respCorrect == 1 & imgs(ii).idx);
        misses = sum(events.respCorrect == 0 & imgs(ii).idx);
    end
    hitsMissesImg(1, ii) = hits;
    hitsMissesImg(2, ii) = misses;
    hitsMissesImg(3, ii) = binocdf(hits, hits + misses, 0.5, 'upper');

    % calculate per run
    for rr = 1:length(runs)
        if ii == 1
            hitsRun = sum(events.respCorrectOrig == 1 & imgs(ii).idx & events.run == runs(rr));
            missesRun = sum(events.respCorrectOrig == 0 & imgs(ii).idx & events.run == runs(rr));
        else
            hitsRun = sum(events.respCorrect == 1 & imgs(ii).idx & events.run == runs(rr));
            missesRun = sum(events.respCorrect == 0 & imgs(ii).idx & events.run == runs(rr));
        end
        hitsMissesRuns(rr*2-1, ii) = hitsRun;
        hitsMissesRuns(rr*2, ii) = missesRun;
    end
    
end
acc = 100*hitsMissesImg(1, :) ./ sum(hitsMissesImg([1, 2], :));
fprintf('Accuracies: '); fprintf('%0.2f, ', acc); fprintf('\n');

% save bar graph of accuracy by condition
figure('Position', [200, 200, 400, 400]); hold on
xPos = [1, 2.5:4.5, 6:8];
bar(xPos, acc, 'FaceColor', [0.5, 0.5, 0.5]);
yline(50, 'Color', 'r');
ylim([0, 100]);
set(gca, 'xtick', xPos, 'xticklabels', {imgs.label});
if saveOutputs
    saveas(gcf, fullfile(outdir, 'response_accuracy'), 'png');
    saveas(gcf, fullfile(outdir, 'response_accuracy'), 'svg');
end

% save output to text file
if saveOutputs
    f = fopen(fullfile(outdir, 'response_accuracy.tsv'), 'w');
    fprintf(f, '\t'); fprintf(f, '%s\t', imgs.label); fprintf(f, '\n');
    fprintf(f, 'Overall hits:\t'); fprintf(f, '%d\t', hitsMissesImg(1, :)); fprintf(f, '\n');
    fprintf(f, 'Overall misses:\t'); fprintf(f, '%d\t', hitsMissesImg(2, :)); fprintf(f, '\n');
    fprintf(f, 'Overall accuracy (%%):\t'); fprintf(f, '%0.2f\t', acc); fprintf(f, '\n');
    fprintf(f, 'Overall binominial p:\t');fprintf(f, '%0.2e\t', hitsMissesImg(3, :)); fprintf(f, '\n');
    for rr = 1:length(runs)
        fprintf(f, '\n');
        fprintf(f, 'Run %d hits:\t', runs(rr)); fprintf(f, '%d\t', hitsMiswsesRuns(rr*2-1, :)); fprintf(f, '\n');
        fprintf(f, 'Run %d misses:\t', runs(rr)); fprintf(f, '%d\t', hitsMissesRuns(rr*2, :)); fprintf(f, '\n');
    end
    fclose(f);
end


%% 3a) Calculate univariate predictors of response time

predsRT = {'run', 'onset', 'img', 'site1-stim200', 'site1-stim100', 'site1-stim0', 'site2-stim200', 'site2-stim100', 'site2-stim0'}; % predictors
Bpreds = array2table(nan(length(predsRT), 7), 'VariableNames', {'Coef', 'Conf1', 'Conf2', 'dfnum', 'dfden', 'F', 'p'});

% which indices to use: remove ones with no responses and ones not assigned an e2v
validIdxes = ~isnan(events.respTime) & ~isnan(events.e2vType);

% response times
Y = events.respTime(validIdxes);

% dummy intercept column
Xint = ones(sum(validIdxes), 1);

% regress by run #
% df: https://math.stackexchange.com/questions/617735/multiple-regression-degrees-of-freedom-f-test#:~:text=The%20correct%20approach%20is%20to,freedom%2C%20i.e.%20n%E2%88%921.
Xrun = [Xint, events.run(validIdxes)]; rr = 1;
[b, bint, ~, ~, stats] = regress(Y, Xrun);
Bpreds.Coef(rr) = b(2); Bpreds{rr, {'Conf1', 'Conf2'}} = bint(2, :); Bpreds{rr, {'dfnum', 'dfden'}} = [1, length(Y) - 2]; Bpreds{rr, {'F', 'p'}} = stats(2:3);

% regress by onset
Xonset = [Xint, events.onset(validIdxes)]; rr = 2;
[b, bint, ~, ~, stats] = regress(Y, Xonset);
Bpreds.Coef(rr) = b(2); Bpreds{rr, {'Conf1', 'Conf2'}} = bint(2, :); Bpreds{rr, {'dfnum', 'dfden'}} = [1, length(Y) - 2]; Bpreds{rr, {'F', 'p'}} = stats(2:3);

% regress by image (use anova)
Ximg = events.trial_type(validIdxes); rr = 3;
[p, tbl, stats] = anova1(Y, Ximg, 'off');
Bpreds{rr, {'dfnum', 'dfden'}} = [6, length(Y) - 7]; Bpreds.F(rr) = tbl{2, 5}; Bpreds.p(rr) = tbl{2, 6};

% regress by each stimsite-ISI binary predictor
for ss = 1:2 % stim site

    % Sham vs 200msISI (at the current stimsite). Pull out indices that are either sham or 200ms at current stimsite
    validIdxes200 = ~isnan(events.respTime) & (events.e2vType == 1 | events.e2vType == 4 & strcmpi(events.electrical_stimulation_site, stimsites{ss}));
    Ystim = events.respTime(validIdxes200);
    Xstim = events.e2vType(validIdxes200); Xstim(Xstim == 1) = 0; Xstim(Xstim == 4) = 1; % pull out e2vtype, but dummycode sham as 0 and stim as 1
    assert(length(unique(Xstim)) == 2, 'Error: Unaccounted e2vTypes');
    Xstim = [ones(sum(validIdxes200), 1), Xstim];
    [b, bint, ~, ~, stats] = regress(Ystim, Xstim);
    rr = 3 + (ss-1)*3 + 1; % stimsite1 = rows 4:6, stimsite2 = rows 7:9
    Bpreds.Coef(rr) = b(2); Bpreds{rr, {'Conf1', 'Conf2'}} = bint(2, :); Bpreds{rr, {'dfnum', 'dfden'}} = [1, length(Ystim) - 2]; Bpreds{rr, {'F', 'p'}} = stats(2:3);

    % Sham vs 100msISI (at the current stimsite). Pull out indices that are either sham or 100ms at current stimsite
    validIdxes100 = ~isnan(events.respTime) & (events.e2vType == 1 | events.e2vType == 3 & strcmpi(events.electrical_stimulation_site, stimsites{ss}));
    Ystim = events.respTime(validIdxes100);
    Xstim = events.e2vType(validIdxes100); Xstim(Xstim == 1) = 0; Xstim(Xstim == 3) = 1; % pull out e2vtype, but dummycode sham as 0 and stim as 1
    assert(length(unique(Xstim)) == 2, 'Error: Unaccounted e2vTypes');
    Xstim = [ones(sum(validIdxes100), 1), Xstim];
    [b, bint, ~, ~, stats] = regress(Ystim, Xstim);
    rr = 3 + (ss-1)*3 + 2; % stimsite1 = rows 4:6, stimsite2 = rows 7:9
    Bpreds.Coef(rr) = b(2); Bpreds{rr, {'Conf1', 'Conf2'}} = bint(2, :); Bpreds{rr, {'dfnum', 'dfden'}} = [1, length(Ystim) - 2]; Bpreds{rr, {'F', 'p'}} = stats(2:3);

    % Sham vs 0msISI (at the current stimsite). Pull out indices that are either sham or 0ms at current stimsite
    validIdxes0 = ~isnan(events.respTime) & (events.e2vType == 1 | events.e2vType == 2 & strcmpi(events.electrical_stimulation_site, stimsites{ss}));
    Ystim = events.respTime(validIdxes0);
    Xstim = events.e2vType(validIdxes0); Xstim(Xstim == 1) = 0; Xstim(Xstim == 2) = 1; % pull out e2vtype, but dummycode sham as 0 and stim as 1
    assert(length(unique(Xstim)) == 2, 'Error: Unaccounted e2vTypes');
    Xstim = [ones(sum(validIdxes0), 1), Xstim];
    [b, bint, ~, ~, stats] = regress(Ystim, Xstim);
    rr = 3 + (ss-1)*3 + 3; % stimsite1 = rows 4:6, stimsite2 = rows 7:9
    Bpreds.Coef(rr) = b(2); Bpreds{rr, {'Conf1', 'Conf2'}} = bint(2, :); Bpreds{rr, {'dfnum', 'dfden'}} = [1, length(Ystim) - 2]; Bpreds{rr, {'F', 'p'}} = stats(2:3);

end

% save output as table
if saveOutputs, writetable(Bpreds, fullfile(outdir, 'univariate_response_times.tsv'), 'FileType', 'text', 'Delimiter', '\t'); end

%% 3b) Calculate multivariate predictors of response time using stim vs sham, adjusting for run, onset, image

% T contains the response variable and all predictors
idxesValid = ~isnan(events.respTime) & ~isnan(events.e2vType); % had response and belongs to 1 of the 4 e2vtypes
T = events(idxesValid, {'respTime', 'run', 'onset', 'trial_type', 'electrical_stimulation_site', 'e2vType'});
T.stim_cond = arrayfun(@(ii) sprintf('%s_%d', T.electrical_stimulation_site{ii}, T.e2vType(ii)), 1:height(T), 'UniformOutput', false)'; % join stim site and e2v
T.stim_cond(strcmpi(T.stim_cond, '_1')) = {'1'}; % remove the underscore so it can come first (be coded as intercept)

% convert to categorical
T.trial_type = categorical(T.trial_type);
T.stim_cond = categorical(T.stim_cond);
T.run = T.run - 1; % so the first run is intercept

% run model
mdlMulti = fitlm(T, 'respTime ~ run + onset + trial_type + stim_cond');
BpredsMulti = mdlMulti.Coefficients;

dfe = mdlMulti.DFE;

% add confidence interval to coefs
BpredsMulti.conf1 = BpredsMulti{:, 'Estimate'} + tinv(0.025, dfe)*BpredsMulti{:, 'SE'};
BpredsMulti.conf2 = BpredsMulti{:, 'Estimate'} + tinv(0.975, dfe)*BpredsMulti{:, 'SE'};
BpredsMulti = BpredsMulti(:, {'Estimate', 'conf1', 'conf2', 'tStat', 'pValue'}); % rearrange cols

if saveOutputs, writetable(BpredsMulti, fullfile(outdir, sprintf('multivariate_response_times_dfe=%d.tsv', dfe)), 'FileType', 'text', 'Delimiter', '\t'); end


%% 4a) Logistic models to predict response accuracy

% pull out all events with button press and has a stimtype (to be consistent with stim analysis)
T = events(~isnan(events.respTime) & ~isnan(events.e2vType), :);

% Create the actual numeric coherence values for fitting.
T.coh = zeros(height(T), 1);
T.coh(contains(T.trial_type, 'low')) = 25;
T.coh(contains(T.trial_type, 'med')) = 50;
T.coh(contains(T.trial_type, 'high')) = 100;

% Change coherence to noise values (or not)
%T.coh = 100 - T.coh;

% Pull the image into a separate col
T.img = cellfun(@(s) s(1:4), T.trial_type_sep, 'UniformOutput', false);

% Set run number to start at 1
T.run = T.run - 1;

% Use a logistic regression model to determine if image affects accuracy. If not we can pool both images
% Recall: estimates refer to ln-odds-ratio. predict(mdl, [0s]) is same as exp(intercept) / (1 + exp(intercept))
mdlLog1 = fitglm(T, 'respCorrectOrig ~ coh', 'Distribution', 'Binomial'); % simple model using just coherence. Sanity check that coherence affects accuracy
mdlLog2 = fitglm(T, 'respCorrectOrig ~ coh + img', 'Distribution', 'Binomial'); % check if image influences accuracy
mdlLog3 = fitglm(T, 'respCorrectOrig ~ coh*img', 'Distribution', 'Binomial'); % check if image influences accuracy, interaction

% Calculate 95% conf intervals of model 3 fit and save. note that the "tStat" in the tables is actually a z-stat
mdlLog3Tbl = mdlLog3.Coefficients;
mdlLog3Tbl = renamevars(mdlLog3Tbl, 'tStat', 'zStat');
mdlLog3Tbl.conf1 = mdlLog3Tbl.Estimate + mdlLog3Tbl.SE * norminv(0.025);
mdlLog3Tbl.conf2 = mdlLog3Tbl.Estimate + mdlLog3Tbl.SE * norminv(0.975);
mdlLog3Tbl = mdlLog3Tbl(:, {'Estimate', 'conf1', 'conf2', 'zStat', 'pValue'}); % resort
if saveOutputs, writetable(mdlLog3Tbl, fullfile(outdir, 'mdl_logistic_coh_img'), 'FileType', 'text', 'Delimiter', '\t'); end

% Test if stim condition predicts accuracy with logistic function, also add in other regressors to adjust (run, onset)
T.stim_cond = arrayfun(@(ii) sprintf('%s_%d', T.electrical_stimulation_site{ii}, T.e2vType(ii)), 1:height(T), 'UniformOutput', false)'; % join stim site and e2v
T.stim_cond(strcmpi(T.stim_cond, '_1')) = {'1'}; % remove the underscore so it can come first (be coded as intercept)
T.stim_cond = categorical(T.stim_cond);
[~, ix] = sort(T.stim_cond);
Tsorted = T(ix, :); % sort so '1' comes first

mdlLog4 = fitglm(Tsorted, 'respCorrectOrig ~ coh + run + onset + stim_cond', 'Distribution', 'Binomial');
mdlLog4Tbl = mdlLog4.Coefficients;
mdlLog4Tbl = renamevars(mdlLog4Tbl, 'tStat', 'zStat');
mdlLog4Tbl.conf1 = mdlLog4Tbl.Estimate + mdlLog4Tbl.SE * norminv(0.025);
mdlLog4Tbl.conf2 = mdlLog4Tbl.Estimate + mdlLog4Tbl.SE * norminv(0.975);
mdlLog4Tbl = mdlLog4Tbl(:, {'Estimate', 'conf1', 'conf2', 'zStat', 'pValue'}); % resort
if saveOutputs, writetable(mdlLog4Tbl, fullfile(outdir, 'mdl_logistic_full'), 'FileType', 'text', 'Delimiter', '\t'); end

% mdlLog5: Merge mdlLog3 and mdlLog4 together (just add img to mdlLog4)
mdlLog5 = fitglm(Tsorted, 'respCorrectOrig ~ img + coh + run + onset + stim_cond', 'Distribution', 'Binomial');
mdlLog5Tbl = mdlLog5.Coefficients;
mdlLog5Tbl = renamevars(mdlLog5Tbl, 'tStat', 'zStat');
mdlLog5Tbl.conf1 = mdlLog5Tbl.Estimate + mdlLog5Tbl.SE * norminv(0.025);
mdlLog5Tbl.conf2 = mdlLog5Tbl.Estimate + mdlLog5Tbl.SE * norminv(0.975);
mdlLog5Tbl = mdlLog5Tbl(:, {'Estimate', 'conf1', 'conf2', 'zStat', 'pValue'}); % resort
if saveOutputs, writetable(mdlLog5Tbl, fullfile(outdir, 'mdl_logistic_full_img'), 'FileType', 'text', 'Delimiter', '\t'); end


%% 4a2 Chi-square test of independence on 25, 50, 100 coh cases: H: proportion correct ~ stim?
% performing across all 7 stim conds at the same time, then bonferroni correcting


%Tsite = T(T.e2vType == 1 | strcmpi(T.electrical_stimulation_site, stimsites{ss}), :);

% Chi-sq on all stim conditions together, for each coherence level
for cc = 2:4
    [tbl, chi2, p, labels] = crosstab(T.respCorrect(T.coh == cohs(cc)), T.stim_cond(T.coh == cohs(cc)));
    tbl = array2table(tbl, 'rownames', labels(1:2, 1), 'variablenames', labels(:, 2));
    disp(tbl);
    fprintf('coh = %d, chi-2 = %0.2f, p=%0.2e\n\n\n', cohs(cc), chi2, p);
end

%% 4b) Fit psychometric (Weibull) curves to accuracy data

% First fit Weibull function to all trials regardless of stim condition

% calculate accuracy for each coherence level from T
accs = zeros(4, 1);
for ii = 1:4
    accs(ii) = sum(T.respCorrectOrig == 1 & T.coh == cohs(ii)) / sum(T.coh == cohs(ii));
end

% fit
ft = fittype('0.5 + 0.5*(1 - exp(-(x/a)^b))', 'independent', 'x', 'coefficients', {'a', 'b'}); % 0.5 is lower asymptote
initialPts = [40, 2]; % initial guess points. First param = coherence at ~63% correct
lb = [0, 0];
fitOptions = fitoptions('Method', 'NonlinearLeastSquares', 'StartPoint', initialPts, 'Lower', lb);
[weibullFit, gof] = fit(cohs, accs, ft, fitOptions);

% get fit over range and confidence intervals
xvals = (0:100)';
[yci, ypred] = predint(weibullFit, xvals, 0.95, 'observation', 'off'); % observation bounds

figure('Position', [200, 200, 400, 400]); hold on
plot(xvals, ypred, 'k-', 'Linewidth', 1.5);
plot(cohs, accs, 'k.', 'MarkerSize', 15);
xlim([-10, 110]); ylim([0.4, 1.1]);
xlabel('Coherence'); ylabel('Accuracy');

% Now get fit for one stim site, iterated across ISI

% calculate accuracy by ISI for a chosen stim site
ss = 1; % choose stim site and pull out table subset
Tsite = T(T.e2vType == 1 | strcmpi(T.electrical_stimulation_site, stimsites{ss}), :);

accsByStim = zeros(4, 4); % 4 columns correspond to sham, 0ms, 100ms, and 200ms
for ee = 1:4
    for ii = 1:4
            accsByStim(ii, ee) = sum(Tsite.e2vType == ee & Tsite.respCorrectOrig == 1 & Tsite.coh == cohs(ii)) / ...
                sum(Tsite.e2vType == ee & Tsite.coh == cohs(ii));
    end
end

% do fits
weibullCoefs = zeros(2, 4); % a, and b from weibull fit, for each of the 4 conditions
ypreds = zeros(length(xvals), 4); % predicted yvalues
for ee = 1:4
    weibullFit = fit(cohs, accsByStim(:, ee), ft, fitOptions);
    weibullCoefs(1, ee) = weibullFit.a;
    weibullCoefs(2, ee) = weibullFit.b;
    ypreds(:, ee) = feval(weibullFit, xvals); % just evaluate the fit over the range
end

% Plot
figure('Position', [200, 200, 400, 400]); hold on
for ee = 1:4
    plot(cohs, accsByStim(:, ee), '.', 'Color', isiCols(ee, :), 'MarkerSize', 15);
    plot(xvals, ypreds(:, ee), 'Color', isiCols(ee, :), 'Linewidth', 1.5);
end
xlim([-10, 110]); ylim([0.2, 1.1]);
set(gca, 'xtick', cohs);
xlabel('Coherence'); ylabel('Accuracy');
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('weibull_byISI_%s', stimsites{ss})), 'png');
    saveas(gcf, fullfile(outdir, sprintf('weibull_byISI_%s', stimsites{ss})), 'svg');
end


%% 4c) significance testing that stim conditions yielded different weibull fits in 4b

weibullVars = var(weibullCoefs, [], 2)';

% Statistical testing on the variances of fitted alpha and beta
nperms = 10000;

if exist(fullfile(outdir, sprintf('weibull_null_%s_nperms-%d.mat', stimsites{ss}, nperms)), 'file')
    fprintf('Loading Weibull variance null distribution from file\n');
    load(fullfile(outdir, sprintf('weibull_null_%s_nperms-%d.mat', stimsites{ss}, nperms)));
else
    fprintf('CALCULATING Weibull variance null distribution for nperms=%d\n', nperms);
    rng('default');
    weibullVarsNull = zeros(nperms, 2); % col 1= var of a, col2 = var of b
    nn = 1; % counter, only increments if fitting works
    while nn <= nperms
        if ~mod(nn, 100), showdot(nn/100, 10); end
    
        % randperm the e2vType. Do we need to do separately for each coherence or no?
        TsitePerm = Tsite;
        TsitePerm.e2vType = TsitePerm.e2vType(randperm(length(TsitePerm.e2vType)));
    
        % calculate accuracies as above
        accsByStimPerm = zeros(4, 4); % 4 columns correspond to sham, 0ms, 100ms, and 200ms
        for ee = 1:4
            for ii = 1:4
                    accsByStimPerm(ii, ee) = sum(TsitePerm.e2vType == ee & TsitePerm.respCorrectOrig == 1 & TsitePerm.coh == cohs(ii)) / ...
                        sum(TsitePerm.e2vType == ee & TsitePerm.coh == cohs(ii));
            end
        end

        if any(isnan(accsByStimPerm(:)))
            fprintf('Encountered nan, retrying');
            continue
        end
    
        % Fit null to weibull
        weibullCoefsPerm = zeros(2, 4); % a, and b from weibull fit, for each of the 4 conditions
        for ee = 1:4
            weibullFit = fit(cohs, accsByStimPerm(:, ee), ft, fitOptions);
            weibullCoefsPerm(1, ee) = weibullFit.a;
            weibullCoefsPerm(2, ee) = weibullFit.b;
        end
    
        weibullVarsNull(nn, :) = var(weibullCoefsPerm, [], 2);
        nn = nn + 1; % iterate counter if successful
    end
    fprintf('\n');
    
    % Save null to file
    save(fullfile(outdir, sprintf('weibull_null_%s_nperms-%d.mat', stimsites{ss}, nperms)), 'weibullVarsNull');
end

% calculate p for alpha and beta separately, right tailed -- looking for probability of observing equal or higher variance in null
ps = [sum(weibullVarsNull(:, 1) >= weibullVars(1)) / length(weibullVarsNull), sum(weibullVarsNull(:, 2) >= weibullVars(2)) / length(weibullVarsNull)];
pjoint = sum(weibullVarsNull(:, 1) >= weibullVars(1) & weibullVarsNull(:, 2) >= weibullVars(2)) / length(weibullVarsNull);