function analysis_bcfs_value

% Analysis of the raw data set of the article "Associations with 
% monetary values do not influence access to awareness for faces"

% This script requires the Matlab Statistics Toolbox and the Bayes Factor
% Package (https://github.com/klabhub/bayesFactor).

%% Define general variables

% cutoff for anticipatory responses (in ms)
aresp_cutoff = 200;

% number of time segments for time-dependent analysis of post-conditioning
% data
n_timeseg = 2;

% number of quantiles for hierarchical shift functions (see Rousselet &
% Wilcox, 2019, for further details)
quant=10;


%% Read in data

% select file
[file,path] = uigetfile({'*.txt';'*.tsv'},'Select raw data file');

data = readmatrix(fullfile(path,file),'NumHeaderLines',1,'FileType','text');

% read header
file_id = fopen(fullfile(path,file));
header = textscan(file_id,'%s',size(data,2),'Delimiter','\t');
fclose(file_id);

% get column numbers for the relevant variables
col_subject = find(cellfun(@(x) contains(x,'SubjectID'),header{1}));
col_session = find(cellfun(@(x) contains(x,'BCFS_session'),header{1}));
col_value = find(cellfun(@(x) contains(x,'Face value'),header{1}));
col_loc = find(cellfun(@(x) contains(x,'Target location'),header{1}));
col_loc_resp = find(cellfun(@(x) contains(x,'Response: location'),header{1}));
col_rt = find(cellfun(@(x) contains(x,'Response: time'),header{1}));
col_trial = find(cellfun(@(x) contains(x,'Trial_number'),header{1}));
col_faceid = find(cellfun(@(x) contains(x,'Face_exemplar'),header{1}));
col_cond = find(cellfun(@(x) contains(x,'QueryTrialsAccuracy'),header{1}));

%% Trial exclusion

% Two types of trials are excluded:
%   - incorrect (or missing) responses
%   - anticipatory (i.e. too fast) responses

% add another column to the data matrix that indicates whether a trial will
% be included or excluded
trial_incl = ones(size(data,1),1);
trial_incl(find(data(:,col_rt) <= aresp_cutoff)) = 0;
trial_incl(find(data(:,col_loc) ~= data(:,col_loc_resp))) = 0;
trial_incl(find(isnan(data(:,col_rt)))) = 0;

data = [data trial_incl];
header{1} = [header{1}; {'inclusion'}];

%% Divide data matrix into sessions

% division into cell array such that there is one cell per session
cell_sizes = find(abs(diff([data(:,col_session); 1])));
cell_sizes = [cell_sizes(1); diff(cell_sizes)];

data_session = mat2cell(data,cell_sizes,size(data,2));

%% Compute median RT and number of missing responses for each stimulus
%% value per session

% the variable "median_rts" contains the subject number (1st column),
% session number (2nd column), and the median for the different face values
% in the following order: [-2, -.1, +.1, +2]
median_rts = cell2mat(cellfun(@(x) [x(1,col_subject) ...
    x(1,col_session) median_rt_value(x,col_value,col_rt)], ...
    data_session,'UniformOutput',false));

%% Compute the number of missing and incorrect responses

% the variable "missing_incorr_resp" contains the subject number (1st column),
% session number (2nd column), the number of missing responses for the
% different face values in the following order: [-2, -.1, +.1, +2], and the
% number of incorrect responses for the different values in the same order
missing_incorr_resp = cell2mat(cellfun(@(x) [x(1,col_subject) ...
    x(1,col_session) missing_resp_value(x,col_value,col_loc_resp,col_loc)], ...
    data_session,'UniformOutput',false));

% number and proportion of missing responses across values
missing_resp = cell2mat(cellfun(@(x,y) [sum(x(3:6)) sum(x(3:6))/size(y,1)],...
    num2cell(missing_incorr_resp,2),data_session,'UniformOutput',false));

% number and proportion of incorrect responses across values
incorr_resp = cell2mat(cellfun(@(x,y) [sum(x(7:end)) sum(x(7:end))/size(y,1)],...
    num2cell(missing_incorr_resp,2),data_session,'UniformOutput',false));

% compute mean and standard deviation of missing responses across participants
missing_resp_subj = cellfun(@(x) mean(x(:,2)),mat2cell(missing_resp,...
    2*ones(size(missing_resp,1)/2,1),size(missing_resp,2)));
missing_resp_mean_sd = [mean(missing_resp_subj) std(missing_resp_subj)];

% compute mean and standard deviation of incorrect responses across participants
incorr_resp_subj = cellfun(@(x) mean(x(:,2)),mat2cell(incorr_resp,...
    2*ones(size(incorr_resp,1)/2,1),size(incorr_resp,2)));
incorr_resp_mean_sd = [mean(incorr_resp_subj) std(incorr_resp_subj)];




%% Rearrange data matrix to have one cell per subject

median_rts = mat2cell(median_rts,2*ones(size(median_rts,1)/2,1),...
    size(median_rts,2));

%% Get baseline median RTs

median_rts_baseline = cell2mat(cellfun(@(x) x(1,:),median_rts,...
    'UniformOutput',false));

%% Compute inference statistics for baseline RTs

% compute mean for high and low reward
stats_reward.baseline_mean.high_reward = mean(median_rts_baseline(:,6));
stats_reward.baseline_mean.low_reward = mean(median_rts_baseline(:,5));

% compute standard errors of the mean for high and low reward
stats_reward.baseline_sem.high_reward = ...
    stats_reward.baseline_mean.high_reward/sqrt(size(median_rts_baseline,1));
stats_reward.baseline_sem.low_reward = ...
    stats_reward.baseline_mean.low_reward/sqrt(size(median_rts_baseline,1));

% compute p-value and Bayes Factor for reward-related faces in baseline
% phase
[~,stats_reward_baseline.p,stats_reward_baseline.ci,...
    stats_reward_baseline.stats] = ttest(median_rts_baseline(:,5),...
    median_rts_baseline(:,6));
[stats_reward_baseline.bf10,~,~,~] = ...
    bf.ttest(median_rts_baseline(:,5)-median_rts_baseline(:,6));
stats_reward_baseline.bf01 = 1/stats_reward_baseline.bf10;


% compute mean change for high and low punishment
stats_pun.baseline_mean.high_pun = mean(median_rts_baseline(:,3));
stats_pun.baseline_mean.low_pun = mean(median_rts_baseline(:,4));

% compute standard errors of the mean for high and low punishment
stats_pun.baseline_sem.high_pun = stats_pun.baseline_mean.high_pun/...
    sqrt(size(median_rts_baseline,1));
stats_pun.baseline_sem.low_pun = stats_pun.baseline_mean.low_pun/...
    sqrt(size(median_rts_baseline,1));

% compute p-value and Bayes Factor for punishment-related faces in baseline
% phase
[~,stats_pun_baseline.p,stats_pun_baseline.ci,stats_pun_baseline.stats] = ...
    ttest(median_rts_baseline(:,3),median_rts_baseline(:,4));
[stats_pun_baseline.bf10,~,~,~] = ...
    bf.ttest(median_rts_baseline(:,3)-median_rts_baseline(:,4));
stats_pun_baseline.bf01 = 1/stats_pun_baseline.bf10;





%% Compute RT differences between sessions and perform inference statistics

median_rts_sessiondiff = cell2mat(cellfun(@(x) x(2,:)-x(1,:),median_rts,...
    'UniformOutput',false));

% compute mean change for high and low reward
stats_reward.mean_change.high_reward = mean(median_rts_sessiondiff(:,6));
stats_reward.mean_change.low_reward = mean(median_rts_sessiondiff(:,5));

% compute standard errors of the mean for the change for high and low
% reward
stats_reward.sem_change.high_reward = ...
    stats_reward.mean_change.high_reward/sqrt(size(median_rts_sessiondiff,1));
stats_reward.sem_change.low_reward = ...
    stats_reward.mean_change.low_reward/sqrt(size(median_rts_sessiondiff,1));

% compute p-value and Bayes Factor for reward-related faces
[~,stats_reward.p,stats_reward.ci,stats_reward.stats] = ...
    ttest(median_rts_sessiondiff(:,5),median_rts_sessiondiff(:,6));
[stats_reward.bf10,~,~,~] = ...
    bf.ttest(median_rts_sessiondiff(:,5)-median_rts_sessiondiff(:,6));
stats_reward.bf01 = 1/stats_reward.bf10;

% compute mean change for high and low punishment
stats_pun.mean_change.high_pun = mean(median_rts_sessiondiff(:,3));
stats_pun.mean_change.low_pun = mean(median_rts_sessiondiff(:,4));

% compute standard errors of the mean for the change for high and low
% punishment
stats_pun.sem_change.high_pun = ...
    stats_pun.mean_change.high_pun/sqrt(size(median_rts_sessiondiff,1));
stats_pun.sem_change.low_pun = ...
    stats_pun.mean_change.low_pun/sqrt(size(median_rts_sessiondiff,1));

% compute p-value and Bayes Factor for punishment-related faces
[~,stats_pun.p,stats_pun.ci,stats_pun.stats] = ...
    ttest(median_rts_sessiondiff(:,3),median_rts_sessiondiff(:,4));
[stats_pun.bf10,~,~,~] = ...
    bf.ttest(median_rts_sessiondiff(:,3)-median_rts_sessiondiff(:,4));
stats_pun.bf01 = 1/stats_pun.bf10;


%% Produce figure for reward-related RT differences

figure;

% define range of axes
axis_min = min(min(median_rts_sessiondiff(:,3:end)));
axis_min = floor(axis_min/250)*250;
axis_max = max(max(median_rts_sessiondiff(:,3:end)));
axis_max = ceil(axis_max/250)*250;

% plot data
patch([axis_min; axis_max; axis_min],[axis_min; axis_max;axis_max],[.85 .85 .85]); hold on
plot(median_rts_sessiondiff(:,6),median_rts_sessiondiff(:,5),'o',...
    'MarkerSize',6,'MarkerFaceColor',[.55 .55 .55],'MarkerEdgeColor','w');

% label axes
xlabel('\Delta RT (post - pre) high reward (ms)','FontWeight','bold',...
    'Color',[.8 .1 .7]);
ylabel('\Delta RT (post - pre) low reward (ms)','FontWeight','bold',...
    'Color',[.1 .7 .8]);
set(gca,'FontSize',8.5);
set(gca,'LineWidth',1.5);

% include diagonal line
line([axis_min axis_max], [axis_min axis_max],'Color','black');

% include means for each condition
line([mean(median_rts_sessiondiff(:,6)) mean(median_rts_sessiondiff(:,6))],...
    [axis_min mean(median_rts_sessiondiff(:,6))],...
    'LineStyle','--','Color',[.8 .1 .7]);

line([axis_min mean(median_rts_sessiondiff(:,5))], ...
    [mean(median_rts_sessiondiff(:,5)) mean(median_rts_sessiondiff(:,5))],...
    'LineStyle','--','Color',[.1 .7 .8]);


xlim([axis_min axis_max]); ylim([axis_min axis_max]);

% modify axes
set(gca,'XTick',axis_min:500:axis_max,'TickDir','out');
set(gca,'YTick',axis_min:500:axis_max,'TickDir','out');

box off

% open dialog box to save figure
[fig_filename,fig_type,user_canceled] = imputfile;

if ~user_canceled   % save figure
    
    [filepath,filename,~] = fileparts(fig_filename);
    
    set(gcf,'PaperPositionMode','auto');
    set(gcf,'PaperUnits','centimeters');
    set(gcf,'PaperPosition',[1 1 8.5 8]);
    
    if strcmp(fig_type,'tif')
        fig_type = 'tiff';
    end
    
    print(gcf,['-d' fig_type],'-r600',fullfile(filepath,filename));
end

close(gcf);


%% Produce figure for punishment-related RT differences

figure;

% plot data (range of axes is taken from reward data)
patch([axis_min; axis_max; axis_min],[axis_min; axis_max;axis_max],[.85 .85 .85]); hold on
plot(median_rts_sessiondiff(:,3),median_rts_sessiondiff(:,4),'o',...
    'MarkerSize',6,'MarkerFaceColor',[.55 .55 .55],'MarkerEdgeColor','w');

% label axes
xlabel('\Delta RT (post - pre) high punishment (ms)','FontWeight','bold',...
    'Color',[.8 .1 .7]);
ylabel('\Delta RT (post - pre) low punishment (ms)','FontWeight','bold',...
    'Color',[.1 .7 .8]);
set(gca,'FontSize',8.5);
set(gca,'LineWidth',1.5);

% include diagonal line
line([axis_min axis_max], [axis_min axis_max],'Color','black');

% include means for each condition
line([mean(median_rts_sessiondiff(:,3)) mean(median_rts_sessiondiff(:,3))],...
    [axis_min mean(median_rts_sessiondiff(:,3))],...
    'LineStyle','--','Color',[.8 .1 .7]);

line([axis_min mean(median_rts_sessiondiff(:,4))], ...
    [mean(median_rts_sessiondiff(:,4)) mean(median_rts_sessiondiff(:,4))],...
    'LineStyle','--','Color',[.1 .7 .8]);


xlim([axis_min axis_max]); ylim([axis_min axis_max]);

% modify axes
set(gca,'XTick',axis_min:500:axis_max,'TickDir','out');
set(gca,'YTick',axis_min:500:axis_max,'TickDir','out');

box off

% open dialog box to save figure
[fig_filename,fig_type,user_canceled] = imputfile;

if ~user_canceled   % save figure
    
    [filepath,filename,~] = fileparts(fig_filename);
    
    set(gcf,'PaperPositionMode','auto');
    set(gcf,'PaperUnits','centimeters');
    set(gcf,'PaperPosition',[1 1 8.5 8]);
    
    if strcmp(fig_type,'tif')
        fig_type = 'tiff';
    end
    
    print(gcf,['-d' fig_type],'-r600',fullfile(filepath,filename));
end

close(gcf);

%% Time-dependent analysis of post-conditioning RT data

% The post-conditioning RT are divided into two halves to examine whether 
% potential effects of learned values might have decayed over time.

% compute median RTs for the different time segments
% different rows for the post-conditioning RTs represent different time
% segments; the columns are: subject ID (first column), session (second
% column), face values in the following order: [-2 -0.1 +0.1 +2] (third to
% sixth column)
median_rts_time = cell2mat(cellfun(@(x)  ...
    median_rt_value_time(x,col_value,col_rt,col_session,...
    col_trial,col_subject,n_timeseg),data_session,'UniformOutput',false));

% rearrange matrix to have one cell per subject
median_rts_time = mat2cell(median_rts_time,...
    (n_timeseg+1)*ones(size(median_rts_time,1)/(n_timeseg+1),1),...
    size(median_rts_time,2));

% subtract baseline from all time segments of post-conditioning data
% (the last column contains the time segment index)
median_rts_time_sessiondiff = cell2mat(cellfun(@(x) [x(2:end,1:2) ...
    x(2:end,3:end) - x(1,3:end) (1:n_timeseg)'], median_rts_time, ...
    'UniformOutput',false));

% prepare data for time x value ANOVA 
% ANOVA factors are: subject (random factor), value (fixed factor), and
% time segment (fixed factor)
anova_factors = {repmat(median_rts_time_sessiondiff(:,1),2,1), ...
    [repmat(1,size(median_rts_time_sessiondiff,1),1); ...
    repmat(2,size(median_rts_time_sessiondiff,1),1)], ...
    repmat(median_rts_time_sessiondiff(:,7),2,1)};
anova_factornames = {'subject','value','time'};
anova_reward = [median_rts_time_sessiondiff(:,5); ...
    median_rts_time_sessiondiff(:,6)];
anova_pun = [median_rts_time_sessiondiff(:,3); ...
    median_rts_time_sessiondiff(:,4)];

% compute ANOVA for reward
[anova_time.reward.p, anova_time.reward.tbl, anova_time.reward.stats] = ...
    anovan(anova_reward, anova_factors, 'model', 'interaction', 'random', 1, ...
    'varnames', anova_factornames, 'display', 'on');

% compute Bayes Factor from F-value
df = cell2mat(anova_time.reward.tbl(find(cellfun(@(x) strcmp(x,'value*time'),...
    anova_time.reward.tbl(:,1))),3));
dfe = cell2mat(anova_time.reward.tbl(find(cellfun(@(x) strcmp(x,'Error'),...
    anova_time.reward.tbl(:,1))),3));
F_reward = cell2mat(anova_time.reward.tbl(find(cellfun(@(x) strcmp(x,'value*time'),...
    anova_time.reward.tbl(:,1))),6));

anova_time.reward.bf10 = 1/(sqrt((length(median_rts_time)^df) * ...
    (1+F_reward*df/dfe)^-(length(median_rts_time))));
anova_time.reward.bf01 = 1/anova_time.reward.bf10;

% compute ANOVA for punishment
[anova_time.pun.p, anova_time.pun.tbl, anova_time.pun.stats] = ...
    anovan(anova_pun, anova_factors, 'model', 'interaction', 'random', 1, ...
    'varnames', anova_factornames, 'display', 'on');

% compute Bayes Factor from F-value
F_pun = cell2mat(anova_time.pun.tbl(find(cellfun(@(x) strcmp(x,'value*time'),...
    anova_time.pun.tbl(:,1))),6));

anova_time.pun.bf10 = 1/(sqrt((length(median_rts_time)^df) * ...
    (1+F_pun*df/dfe)^-(length(median_rts_time))));
anova_time.pun.bf01 = 1/anova_time.pun.bf10;






%% Compute correlation between accuracy in conditioning phase and RTs
%% during bCFS

% The accuracies in the conditioning phase are listed in four columns of the
% data matrix. Get the first line from each participant's data to read
% their accuracies.

subj_ind = unique(data(:,col_subject));

cond_acc = cell2mat(cellfun(@(x) data(find(data(:,col_subject) == x,1),...
    col_cond),num2cell(subj_ind),'UniformOutput',false));

% compute mean accuracy across the four conditioning blocks for each
% participant
mean_cond_acc = mean(cond_acc,2);

% get the difference in the RT change between high and low reward (negative
% values indicate stronger RT reduction for high vs. low reward)
RT_change_diff_reward = median_rts_sessiondiff(:,6) - ...
    median_rts_sessiondiff(:,5);

% get the difference in the RT change between high and low punishment
% (negative values indicate stronger RT reduction for high vs. low punishment)
RT_change_diff_pun = median_rts_sessiondiff(:,3) - ...
    median_rts_sessiondiff(:,4);

% Compute correlation between accuracy in conditioning phase and RT change.
% Since the accuracy data are heavily skewed (only few participants with
% low accuracy) the Spearman correlation coefficient is computed.

% Define ranks for the three variables
ranks = [tiedrank(mean_cond_acc) tiedrank(RT_change_diff_reward) ...
    tiedrank(RT_change_diff_pun)];

% compute Pearson correlation between ranks
[corr_cond_RT.reward.bf10 corr_cond_RT.reward.rho corr_cond_RT.reward.p] = ...
    bf.corr(ranks(:,1),ranks(:,2));
corr_cond_RT.reward.bf01 = 1/corr_cond_RT.reward.bf10;
[corr_cond_RT.pun.bf10 corr_cond_RT.pun.rho corr_cond_RT.pun.p] = ...
    bf.corr(ranks(:,1),ranks(:,3));
corr_cond_RT.pun.bf01 = 1/corr_cond_RT.pun.bf10;

%% Compare change in RTs between reward and punishment (across high and low values)

% the variable "median_rts_val" contains the subject number (1st column),
% session number (2nd column), and the median for the different face
% valences in the following order: [punishment reward]
median_rts_val = cell2mat(cellfun(@(x) [x(1,col_subject) x(1,col_session) ...
    median_rt_valence(x,col_value,col_rt)],data_session,...
    'UniformOutput',false));

% Rearrange data matrix to have one cell per subject
median_rts_val = mat2cell(median_rts_val,2*ones(size(median_rts_val,1)/2,1),...
    size(median_rts_val,2));

% compute difference between sessions
median_rts_val_sessiondiff = cell2mat(cellfun(@(x) x(2,:)-x(1,:),...
    median_rts_val,'UniformOutput',false));


% compute mean change for reward and punishment
stats_valence.mean_change.pun = mean(median_rts_val_sessiondiff(:,3));
stats_valence.mean_change.reward = mean(median_rts_val_sessiondiff(:,4));

% compute standard errors of the mean for the RT change for reward and
% punishment
stats_valence.sem_change.pun = ...
    stats_valence.mean_change.pun/sqrt(size(median_rts_val_sessiondiff,1));
stats_valence.sem_change.reward = ...
    stats_valence.mean_change.reward/sqrt(size(median_rts_val_sessiondiff,1));

% compute p-value and Bayes Factor for the difference between reward and
% punishment
[~,stats_valence.p,stats_valence.ci,stats_valence.stats] = ...
    ttest(median_rts_val_sessiondiff(:,3),median_rts_val_sessiondiff(:,4));
[stats_valence.bf10,~,~,~] = ...
    bf.ttest(median_rts_val_sessiondiff(:,3)-median_rts_val_sessiondiff(:,4));
stats_valence.bf01 = 1/stats_valence.bf10;

% compute mean change for the comparison between high reward and high
% punishment
stats_valence_highval.mean_change.high_pun = mean(median_rts_sessiondiff(:,3));
stats_valence_highval.mean_change.high_reward = mean(median_rts_sessiondiff(:,6));

% compute standard errors of the mean for the RT change for high reward and
% high punishment
stats_valence_highval.sem_change.high_pun = ...
    stats_valence_highval.mean_change.high_pun/sqrt(size(median_rts_sessiondiff,1));
stats_valence_highval.sem_change.high_reward = ...
    stats_valence_highval.mean_change.high_reward/sqrt(size(median_rts_sessiondiff,1));

% compute p-value and Bayes Factor for the difference between high reward and
% high punishment
[~,stats_valence_highval.p,stats_valence_highval.ci,stats_valence_highval.stats] = ...
    ttest(median_rts_sessiondiff(:,3),median_rts_sessiondiff(:,6));
[stats_valence_highval.bf10,~,~,~] = ...
    bf.ttest(median_rts_sessiondiff(:,3)-median_rts_sessiondiff(:,6));
stats_valence_highval.bf01 = 1/stats_valence_highval.bf10;


%% Compare change in RTs between high and low value (across reward and punishment)

% the variable "median_rts_sal" contains the subject number (1st column),
% session number (2nd column), and the median for the different face
% saliences in the following order: [0.1 (i.e. low value) 2 (i.e. high value)]
median_rts_sal = cell2mat(cellfun(@(x) [x(1,col_subject) x(1,col_session) ...
    median_rt_salience(x,col_value,col_rt)],data_session,...
    'UniformOutput',false));

% Rearrange data matrix to have one cell per subject
median_rts_sal = mat2cell(median_rts_sal,2*ones(size(median_rts_sal,1)/2,1),...
    size(median_rts_sal,2));

% compute difference between sessions
median_rts_sal_sessiondiff = cell2mat(cellfun(@(x) x(2,:)-x(1,:),...
    median_rts_sal,'UniformOutput',false));


% compute mean change for high and low values
stats_salience.mean_change.high = mean(median_rts_sal_sessiondiff(:,4));
stats_salience.mean_change.low = mean(median_rts_sal_sessiondiff(:,3));

% compute standard errors of the mean for the RT change for high and low
% values
stats_salience.sem_change.high = ...
    stats_salience.mean_change.high/sqrt(size(median_rts_sal_sessiondiff,1));
stats_salience.sem_change.low = ...
    stats_salience.mean_change.low/sqrt(size(median_rts_sal_sessiondiff,1));

% compute p-value and Bayes Factor for the difference between high and low
% values
[~,stats_salience.p,stats_salience.ci,stats_salience.stats] = ...
    ttest(median_rts_sal_sessiondiff(:,3),median_rts_sal_sessiondiff(:,4));
[stats_salience.bf10,~,~,~] = ...
    bf.ttest(median_rts_sal_sessiondiff(:,3)-median_rts_sal_sessiondiff(:,4));
stats_salience.bf01 = 1/stats_salience.bf10;

%% Hierarchical shift function for the analysis of the RT distributions

% define quantiles as proportion of the whole distribution
quant_perc = linspace(0,1,quant+1);
quant_perc(1) = []; quant_perc(end) = [];

% compute the quantiles for each session and stimulus value
% The first two columns contain the subject and the session number, columns
% 3 to 6 contain the quantiles for the different values in the following
% order: [-2 -0.1 +0.1 +2], the last column contains the number of the
% quantiles.
% If the RTs are computed for four different quantiles, for example, each 
% cell has four rows. The number in the third column of the first row
% contains the first quantile of the RT distribution for the value of -2
% ...
quantile_data = cell2mat(cellfun(@(x) [x(1,col_subject)*...
    ones(length(quant_perc),1) x(1,col_session)*ones(length(quant_perc),1) ...
    compute_quantiles(x,col_value,quant_perc,col_rt)], ...
    data_session, 'UniformOutput',false));

% rearrange data to have one cell per subject
quantile_data = mat2cell(quantile_data,(length(quant_perc)*2*...
    ones(size(quantile_data,1)/(length(quant_perc)*2),1)),...
    size(quantile_data,2));

% compute differences between the two sessions
% (1st column: subject number, 2nd to 5th column: RTs for values in order:
% [-2 -0.1 +0.1 +2], 6th column: quantile number)
quantiles_sessiondiff = cellfun(@(x) [x(1,col_subject)*...
    ones(length(quant_perc),1) x((size(x,1)/2)+1:end,3:6)-...
    x(1:size(x,1)/2,3:6) (1:length(quant_perc))'],quantile_data,...
    'UniformOutput',false);

% get quantile RTs for the different values (rows = subjects, columns =
% quantiles)
quantiles_highrew = cell2mat(cellfun(@(x) (x(:,5))', quantiles_sessiondiff,...
    'UniformOutput',false));
quantiles_highrew = num2cell(quantiles_highrew,1);

quantiles_lowrew = cell2mat(cellfun(@(x) (x(:,4))', quantiles_sessiondiff,...
    'UniformOutput',false));
quantiles_lowrew = num2cell(quantiles_lowrew,1);

quantiles_highpun = cell2mat(cellfun(@(x) (x(:,2))', quantiles_sessiondiff,...
    'UniformOutput',false));
quantiles_highpun = num2cell(quantiles_highpun,1);

quantiles_lowpun = cell2mat(cellfun(@(x) (x(:,3))', quantiles_sessiondiff,...
    'UniformOutput',false));
quantiles_lowpun = num2cell(quantiles_lowpun,1);

% compute mean and SEM for each quantile
quantile_analysis.high_reward.mean = cellfun(@(x) mean(x), ...
    quantiles_highrew);
quantile_analysis.high_reward.sem = cellfun(@(x) std(x)/sqrt(length(x)),...
    quantiles_highrew);

quantile_analysis.low_reward.mean = cellfun(@(x) mean(x), ...
    quantiles_lowrew);
quantile_analysis.low_reward.sem = cellfun(@(x) std(x)/sqrt(length(x)),...
    quantiles_lowrew);

quantile_analysis.high_pun.mean = cellfun(@(x) mean(x), ...
    quantiles_highpun);
quantile_analysis.high_pun.sem = cellfun(@(x) std(x)/sqrt(length(x)),...
    quantiles_highpun);

quantile_analysis.low_pun.mean = cellfun(@(x) mean(x), ...
    quantiles_lowpun);
quantile_analysis.low_pun.sem = cellfun(@(x) std(x)/sqrt(length(x)),...
    quantiles_lowpun);

% perform t-test between high and low reward for each quantile
[~, quantile_analysis.ttest.reward.p,quantile_analysis.ttest.reward.ci,...
    quantile_analysis.ttest.reward.stats] = ...
    cellfun(@(x,y) ttest(x,y), quantiles_highrew, quantiles_lowrew,...
    'UniformOutput',false);
quantile_analysis.ttest.reward.bf10 = ...
    cellfun(@(x,y) bf.ttest(x-y),quantiles_highrew, quantiles_lowrew,...
    'UniformOutput',false);
quantile_analysis.ttest.reward.bf01 = cellfun(@(x) 1/x,...
    quantile_analysis.ttest.reward.bf10, 'UniformOutput',false);

% perform t-test between high and low punishment for each quantile
[~, quantile_analysis.ttest.pun.p,quantile_analysis.ttest.pun.ci,...
    quantile_analysis.ttest.pun.stats] = ...
    cellfun(@(x,y) ttest(x,y), quantiles_highpun, quantiles_lowpun,...
    'UniformOutput',false);
quantile_analysis.ttest.pun.bf10 = ...
    cellfun(@(x,y) bf.ttest(x-y),quantiles_highpun, quantiles_lowpun,...
    'UniformOutput',false);
quantile_analysis.ttest.pun.bf01 = cellfun(@(x) 1/x,...
    quantile_analysis.ttest.pun.bf10, 'UniformOutput',false);

% plot hierarchical shift function for reward
figure;
plot((1:length(quant_perc))-.1,quantile_analysis.high_reward.mean,...
    'LineStyle','-','LineWidth',1.25,'Marker','o','MarkerFaceColor',...
    [0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',4,'Color','k');
hold on
plot((1:length(quant_perc))+.1,quantile_analysis.low_reward.mean,...
    'LineStyle','--','LineWidth',1.25,'Marker','d','MarkerFaceColor',...
    [0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',4,'Color','k');

% add SEMs
cellfun(@(x,y,z) line([x x],[z-y z+y],'Color','k','LineStyle','-',...
    'LineWidth',.5),num2cell((1:length(quant_perc))-.1),...
    num2cell(quantile_analysis.high_reward.sem),...
    num2cell(quantile_analysis.high_reward.mean));
cellfun(@(x,y,z) line([x x],[z-y z+y],'Color','k','LineStyle','-',...
    'LineWidth',.5),num2cell((1:length(quant_perc))+.1),...
    num2cell(quantile_analysis.low_reward.sem),...
    num2cell(quantile_analysis.low_reward.mean));

% axes
set(gca,'xlim',[.75 9.25]);
set(gca,'ylim',[-2000 0]);
set(gca,'XTick',1:length(quant_perc));
xlabel('Deciles','FontWeight','bold','FontSize',10);
ylabel('\Delta RT (post - pre) (ms)','FontWeight','bold','FontSize',10);
set(gca,'FontSize',9,'LineWidth',1.5);
box off;

% add legend 
legend('High reward','Low reward','Location','southwest','FontSize',8);

% add p-values and Bayes Factors
text(-.65,100,['p      ';'BF_{10}'],'FontWeight','bold','FontSize',7);
cellfun(@(x,y,p,b) text(x,y,sprintf('%0.2f\n%0.2f',p,b),...
    'HorizontalAlignment','center','FontSize',7),...
    num2cell(1:length(quant_perc)),...
    num2cell(100*ones(1,length(quant_perc))),...
    quantile_analysis.ttest.reward.p,quantile_analysis.ttest.reward.bf10);


% open dialog box to save figure
[fig_filename,fig_type,user_canceled] = imputfile;

if ~user_canceled   % save figure
    
    [filepath,filename,~] = fileparts(fig_filename);
    
    set(gcf,'PaperPositionMode','auto');
    set(gcf,'PaperUnits','centimeters');
    set(gcf,'PaperPosition',[1 1 8.5 8]);
    
    if strcmp(fig_type,'tif')
        fig_type = 'tiff';
    end
    
    print(gcf,['-d' fig_type],'-r600',fullfile(filepath,filename));
end

close(gcf);

% plot hierarchical shift function for punishment
figure;
plot((1:length(quant_perc))-.1,quantile_analysis.high_pun.mean,...
    'LineStyle','-','LineWidth',1.25,'Marker','o','MarkerFaceColor',...
    [0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',4,'Color','k');
hold on
plot((1:length(quant_perc))+.1,quantile_analysis.low_pun.mean,...
    'LineStyle','--','LineWidth',1.25,'Marker','d','MarkerFaceColor',...
    [0 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',4,'Color','k');

% add SEMs
cellfun(@(x,y,z) line([x x],[z-y z+y],'Color','k','LineStyle','-',...
    'LineWidth',.5),num2cell((1:length(quant_perc))-.1),...
    num2cell(quantile_analysis.high_pun.sem),...
    num2cell(quantile_analysis.high_pun.mean));
cellfun(@(x,y,z) line([x x],[z-y z+y],'Color','k','LineStyle','-',...
    'LineWidth',.5),num2cell((1:length(quant_perc))+.1),...
    num2cell(quantile_analysis.low_pun.sem),...
    num2cell(quantile_analysis.low_pun.mean));

% axes
set(gca,'xlim',[.75 9.25]);
set(gca,'ylim',[-2000 0]);
set(gca,'XTick',1:length(quant_perc));
xlabel('Deciles','FontWeight','bold','FontSize',10);
ylabel('\Delta RT (post - pre) (ms)','FontWeight','bold','FontSize',10);
set(gca,'FontSize',9,'LineWidth',1.5);
box off;

% add legend 
legend('High punishment','Low punishment','Location','southwest','FontSize',8);

% add p-values and Bayes Factors
text(-.65,100,['p      ';'BF_{10}'],'FontWeight','bold','FontSize',7);
cellfun(@(x,y,p,b) text(x,y,sprintf('%0.2f\n%0.2f',p,b),...
    'HorizontalAlignment','center','FontSize',7),...
    num2cell(1:length(quant_perc)),...
    num2cell(100*ones(1,length(quant_perc))),...
    quantile_analysis.ttest.pun.p,quantile_analysis.ttest.pun.bf10);


% open dialog box to save figure
[fig_filename,fig_type,user_canceled] = imputfile;

if ~user_canceled   % save figure
    
    [filepath,filename,~] = fileparts(fig_filename);
    
    set(gcf,'PaperPositionMode','auto');
    set(gcf,'PaperUnits','centimeters');
    set(gcf,'PaperPosition',[1 1 8.5 8]);
    
    if strcmp(fig_type,'tif')
        fig_type = 'tiff';
    end
    
    print(gcf,['-d' fig_type],'-r600',fullfile(filepath,filename));
end

close(gcf);




%% Perform multilevel modelling

% get variables for linear and generalized linear mixed-effects model
data_mem = data;

% Generalized linear mixed-effects models investigating the effect of
%   trial number and face identity on reaction times
data_glme = data_mem;
data_glme(~data_mem(:,end),:) = [];

% get data from all trials for each valence
glme_valence = [1 -1];
for v = 1:length(glme_valence)
    
    % identify trials for this valence
    trial_ind = find(sign(data_glme(:,col_value)) == glme_valence(v));
    
    % get all other data
    mem_subject = data_glme(trial_ind,col_subject);
    mem_session = data_glme(trial_ind,col_session);
    mem_value = data_glme(trial_ind,col_value);
    mem_loc = data_glme(trial_ind,col_loc);
    mem_trial = data_glme(trial_ind,col_trial);
    mem_faceid = data_glme(trial_ind,col_faceid);
    mem_rt = data_glme(trial_ind,col_rt);
    
    % create table
    glme_table = table(mem_rt,mem_faceid,mem_trial,mem_loc,mem_value,...
        mem_session,mem_subject);
    
    % fit model
    glme_rt_formula = (['mem_rt ~ mem_value*mem_session*mem_trial*mem_faceid + ' ...
        '(mem_value*mem_session*mem_trial*mem_faceid | mem_subject)']);
    glme_rt{v} = fitglme(glme_table,glme_rt_formula,...
        'Distribution','Gamma','DummyVarCoding','effects');

    
end
    

end

%% function median_rt_value
function median_rt = median_rt_value(session_data,col_value,col_rt)

% function to compute the median RT for each face value in each session
values = unique(session_data(:,col_value));
median_rt = nan(1,length(values));

for v = 1:length(values)
    
    % trial indices for included trials with this value
    
    trial_ind = find(session_data(:,col_value) == values(v) & ...
        session_data(:,end) == 1);
    
    % compute the median for these trials
    
    median_rt(v) = median(session_data(trial_ind,col_rt));
    
end

end

%% function missing_resp_value
function missing_incorr_resp = ...
    missing_resp_value(session_data,col_value,col_loc_resp,col_loc)

% function to compute the number of missing and incorrect responses for
% each face value in each session
values = unique(session_data(:,col_value));
missing_resp = zeros(1,length(values));
incorr_resp = zeros(1,length(values));

for v = 1:length(values)
    
    % number of missing responses (i.e. trials in which the response
    % location equals 0)
    missing_resp(v) = length(find(~session_data(:,col_loc_resp) & ...
        session_data(:,col_value) == values(v)));
    
    % number of incorrect responses (i.e. trials in which the location
    % indicated by the participant is not zero but does not match the
    % actual location)
    incorr_resp(v) = length(find((session_data(:,col_loc_resp) ~= ...
        session_data(:,col_loc)) & session_data(:,col_loc_resp)>0 & ...
        session_data(:,col_value) == values(v)));
    
end

% combine matrices with number of missing and incorrect responses
missing_incorr_resp = [missing_resp incorr_resp];

end


%% function median_rt_value_time
function median_rt_time = ...
    median_rt_value_time(session_data,col_value,col_rt,col_session,...
    col_trial,col_subject,nr_timeseg)

% function to compute the median RT for each face value in each session,
% separately for each time segment
values = unique(session_data(:,col_value));

% different time segments will only be analyzed for the post-conditioning
% session
if session_data(1,col_session) == 1
    median_rt_time = nan(1,length(values)+2);
    trials_timeseg = {1:size(session_data,1)};
elseif session_data(1,col_session) == 2
    median_rt_time = nan(nr_timeseg,length(values)+2);
    all_trial_ind = 1:size(session_data,1);
    trials_timeseg = mat2cell(all_trial_ind,1,...
        (length(all_trial_ind)/nr_timeseg)*ones(1,nr_timeseg));
end

% include subject number and session number in first two columns 
median_rt_time(:,1) = session_data(1,col_subject);
median_rt_time(:,2) = session_data(1,col_session);

for v = 1:length(values)
    
    for ts = 1:length(trials_timeseg)
    
    % trial indices for included trials with this value
    
    trial_ind = find(session_data(:,col_value) == values(v) & ...
        session_data(:,end) == 1 & ...
        ismember(session_data(:,col_trial),trials_timeseg{ts}));
    
    % compute the median for these trials
    
    median_rt_time(ts,v+2) = median(session_data(trial_ind,col_rt));
    
    end
    
end

end



%% function median_rt_valence

function median_rt = median_rt_valence(session_data,col_value,col_rt)

% function to compute the median RT for each face valence (i.e. reward vs.
% punishment) in each session
valence = unique(sign(session_data(:,col_value)));
median_rt = nan(1,length(valence));

for v = 1:length(valence)
    
    % trial indices for included trials with this value
    
    trial_ind = find(sign(session_data(:,col_value)) == valence(v) & ...
        session_data(:,end) == 1);
    
    % compute the median for these trials
    
    median_rt(v) = median(session_data(trial_ind,col_rt));
    
end

end

%% function median_rt_salience
function median_rt = median_rt_salience(session_data,col_value,col_rt)

% function to compute the median RT for each face salience (i.e. high vs.
% low value)in each session
salience = unique(abs(session_data(:,col_value)));
median_rt = nan(1,length(salience));

for v = 1:length(salience)
    
    % trial indices for included trials with this value
    
    trial_ind = find(session_data(:,col_value) == salience(v) & ...
        session_data(:,end) == 1);
    
    % compute the median for these trials
    
    median_rt(v) = median(session_data(trial_ind,col_rt));
    
end

end

%% function compute_quantiles
function quantile_data = compute_quantiles(data_session,col_value,...
    quant_perc,col_rt)

% compute quantiles for hierarchical shift function

values = unique(data_session(:,col_value));
quantile_data = nan(length(quant_perc),length(values));

for v = 1:length(values)
    
    % trial indices for valid trials with this value
    trial_ind = find(data_session(:,col_value) == values(v) & ...
        data_session(:,end) == 1);
    
    % compute quantiles
    quantile_data(:,v) = quantile(data_session(trial_ind,col_rt),quant_perc);
    
end

% add quantile number to data matrix
quantile_data = [quantile_data (1:length(quant_perc))'];

end