covariates_analysis.m

function [R2_losses] = covariates_analysis(config_file)
%% Covariates Analysis

%{
Created by Luca La Fisca
ISIA Lab, Faculty of Engineering University of Mons, Mons (Belgium)
luca.lafisca@umons.ac.be
Source: La Fisca et. al., "Biases in BCI Experiments: Do we really need to 
balance stimulus properties across categories?", Frontiers in Computational
Neuroscience, Volume 16, Pages 153-172, 22 November 2022
Copyright (C) 2022 - UMons
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
This library cannot be used for commercial use without the agreement of the
author (Luca La Fisca).
%}

%% (DATA SPECIFIC) Set paths and variable names
%read config file
if ~exist('config_file', 'var')
    config_file = 'config.json';
end
fid = fopen(config_file);
raw = fread(fid,inf);
str = char(raw');
fclose(fid);
config = jsondecode(str);

PATH_TO_ROOT                = config.PATH_TO_ROOT;
PATH_TO_SOURCE              = config.PATH_TO_SOURCE;
PATH_TO_FIELDTRIP           = config.PATH_TO_FIELDTRIP;
PATH_TO_LIMO                = config.PATH_TO_LIMO;
PATH_TO_FT2LIMO             = config.PATH_TO_FT2LIMO;
PATH_TO_COV_ANALYSIS        = config.PATH_TO_COV_ANALYSIS;

ext = ['.' config.datatype];
config.ext = ext;
save_choice = config.save_choice;


% output folder (derivatives)
PATH_TO_DERIV               = fullfile(PATH_TO_ROOT, 'derivatives');
config.PATH_TO_DERIV        = PATH_TO_DERIV;

% add toolboxes to path
cd(PATH_TO_FIELDTRIP)
ft_defaults
addpath(PATH_TO_LIMO)
addpath(genpath(fullfile(PATH_TO_LIMO,'external')))
addpath(genpath(fullfile(PATH_TO_LIMO,'limo_cluster_functions')))
addpath(genpath(PATH_TO_FT2LIMO))
addpath(genpath(PATH_TO_COV_ANALYSIS))

BIDS_FOLDER = PATH_TO_ROOT;
config.BIDS_FOLDER = BIDS_FOLDER;

cd(PATH_TO_ROOT)

%% Selection of variables/confounders/covariates (cf. paper section 2.4)
[correlation_val] = var_selection(config);
disp('correlation values:')
disp(correlation_val)

%% Linear Modeling (cf. paper section 2.5)

%% Create the models
% model_names = {'model_cat', 'model_psycho', 'model_image', 'model_psycho_image'};
model_names = split(config.model_names,',');

selected_regressors = {[],4:13,14:26,4:26}; % the regressors corresponding to each model_name
trialinfo_filename = 'trialinfo_psycho_image'; % give the trialinfo corresponding to the complete model

trial_start = config.trial_start; %starting time of the trial in ms
trial_end = config.trial_end; %ending time of the trial in ms

% Information about categorical variables
contrast.mat = [1 -1 0];
regress_cat = {1:2 ,1;
                0  ,0};
            
for i = 1:length(model_names)
    model_name = model_names{i};
    regressors = selected_regressors{i};
    model = create_model(PATH_TO_DERIV,PATH_TO_SOURCE,SOURCE_ANALYSIS,task_name,...
                trialinfo_filename,trial_start,trial_end,regressors,regress_cat);
    if isempty(model.cont_files{1})
        model.cont_files = {};
    end
    if save_choice
        save(fullfile(BIDS_FOLDER,'derivatives',[model_name,'.mat']),'model')
    end
    
    % Remove previous GLM folder (automatically generated by LIMO)
    dinfo = dir(fullfile(BIDS_FOLDER,'sub-*'));
    subj = {dinfo.name};
    for subj_name = subj
        subj_name = subj_name{1};
        del_folder = fullfile(BIDS_FOLDER,'derivatives',subj_name,'eeg','GLM_OLS_Time_Channels');
        [root,name,~] = fileparts(del_folder);
        cd(root)
        if exist(del_folder,'dir')
            rmdir(name,'s')
        end
    end

    % Run limo_batch on complete model
    cd(PATH_TO_ROOT)
    option = 'both'; % or 'model specification', 'contrast only' or 'both'
    [LIMO_files, procstatus] = limo_batch(option,model,contrast);

    % Rename GLM folder
    dinfo = dir(fullfile(BIDS_FOLDER,'sub-*'));
    subj = {dinfo.name};
    for subj_name = subj
        subj_name = subj_name{1};
        del_folder = fullfile(BIDS_FOLDER,'derivatives',subj_name,'eeg','GLM_OLS_Time_Channels');
        [root,name,~] = fileparts(del_folder);
        cd(root)
        if exist(del_folder,'dir')
            new_name = [model_name '_GLM_OLS_Time_Channels'];
            if exist(fullfile(root,new_name),'dir')
                rmdir(new_name,'s')
            end
            movefile(name,new_name);
            cd(fullfile(BIDS_FOLDER,'derivatives',subj_name,'eeg',new_name))
            load('LIMO.mat')
            LIMO.dir = pwd;
            save('LIMO.mat','LIMO')
        end
    end
end

% check the design matrix
tmp = LIMO.design.X;
tmp(tmp(:,1:2)==0) = -1;
figure;imagesc(tmp);colormap(gray);caxis([-1,1]);colorbar;
ylabel('trials / subjects')
title('Check Design Matrix')

% Define the corresponding xtickslabels:

% xticks(1:3);
% xticklabels({'Manufactured\newlineitem', 'Natural\newlineitem', 'Error'});

% xticks(1:11);
% xticklabels({'Manufactured\newlineitem', 'Natural\newlineitem', '#Phoneme\newlinePrimer','#Phoneme\newlineTarget', ...
%     'Movie\newlineFrequency\newlinePrimer', 'Movie\newlineFrequency\newlineTarget', 'AoA\newlinePrimer','AoA\newlineTarget',...
%     'Familiarity\newlinePrimer','Familiarity\newlineTarget', 'Error'});

% xticks(1:16);
% xticklabels({'Manufactured\newlineitem', 'Natural\newlineitem', 'Correlation\newlinePrimer', 'Correlation\newlineTarget',...
%     'Compactness\newlinePrimer', 'Compactness\newlineTarget', 'Ratio\newlinePrimer','Ratio\newlineTarget',...
%     'FreqEnergy\newlinePrimer', 'FreqEnergy\newlineTarget', 'Contrast\newlinePrimer','Contrast\newlineTarget',...
%     'frequency\newlinePrimer', 'frequency\newlineTarget', 'Visual\newlineSimilarity', 'Error'});

xticks(1:24);
xticklabels({'Man.\newlineitem', 'Nat.\newlineitem', '#Phon.\newlinePrimer','#Phon.\newlineTarget', ...
    'MovFreq.\newlinePrimer', 'MovFreq.\newlineTarget', 'AoA\newlinePrimer','AoA\newlineTarget',...
    'Fam.\newlinePrimer','Fam.\newlineTarget',...
    'Corr.\newlinePrimer', 'Corr.\newlineTarget',...
    'Compact.\newlinePrimer', 'Compact.\newlineTarget', 'Ratio\newlinePrimer','Ratio\newlineTarget',...
    'FreqEn.\newlinePrimer', 'FreqEn.\newlineTarget', 'Contr.\newlinePrimer','Contr.\newlineTarget',...
    'freq.\newlinePrimer', 'freq.\newlineTarget', 'Visual\newlineSim.', 'Error'});

%% Create the naive models (it is a long process)
% model_names = {'model_psycho', 'model_image', 'model_psycho_image'}; %Note: a naive categorical model is useless
model_names = split(config.model_names_naive,',');
selected_regressors = {4:13,14:26,4:26}; % the regressors corresponding to each model_name
trialinfo_filename = 'trialinfo_psycho_image'; % give the trialinfo corresponding to the complete model

for k = 1:length(model_names)
    model_name = model_names{k};
    regressors = selected_regressors{k};

    model = load(fullfile(BIDS_FOLDER,'derivatives',[model_name,'.mat']));
    model = model.(cell2mat(fieldnames(model)));
    numberOfCov = size(model.cont_files{1},2);

    cd(PATH_TO_ROOT)
    dinfo = dir(fullfile(PATH_TO_DERIV,'sub-*'));
    subj = {dinfo.name};
    rng('shuffle')
    option = 'both'; % or 'model specification', 'contrast only' or 'both'

    % Initiate the parameters on which the mean will be computed
    tmp_betas = cell(length(subj),1);
    tmp_con = cell(length(subj),1);
    tmp_r2 = cell(length(subj),1);
    betas_mean = cell(length(subj),1);
    con_mean = cell(length(subj),1);
    r2_mean = cell(length(subj),1);

    % repeat 30 times to get mean naive (/!\ this operation takes a lot of time!)
    for j=1:30
        i = 1;
        for subj_name = subj
            subj_name = subj_name{1};
            load(fullfile(PATH_TO_DERIV,subj_name, 'eeg', [model_name '_GLM_OLS_Time_Channels'], 'LIMO.mat'))
            idx = ~isnan(model.cont_files{i});
            model_tmp = model.cont_files{i};
            delete(fullfile(BIDS_FOLDER,'limo_batch_report','GLM_OLS_Time_Channels',['subject' num2str(i)],'PIPE.lock'))
            naive_path = fullfile(PATH_TO_DERIV,subj_name, 'eeg', [model_name '_naive_GLM_OLS_Time_Channels']);
            glm_path = fullfile(PATH_TO_DERIV,subj_name, 'eeg', 'GLM_OLS_Time_Channels');
            if ~exist(glm_path,'dir')
                mkdir(glm_path)
            end
            LIMO.dir = glm_path;
            if isfield(LIMO,'contrast')
                LIMO = rmfield(LIMO,'contrast');
                save(fullfile(glm_path,'LIMO.mat'),'LIMO')
            end
            null_mat = mvnrnd(zeros(1,numberOfCov),cov(LIMO.design.X(:,3:numberOfCov+2)),size(LIMO.design.X,1));
            model_tmp(sum(idx,2)==numberOfCov,:) = null_mat;
            model.cont_files{i} = model_tmp;

            i = i + 1;
        end

        [LIMO_files, procstatus] = limo_batch(option,model,contrast);

        i=1;
        for subj_name = subj
            subj_name = subj_name{1};
            tmp = load(LIMO_files.mat{i});
            load(fullfile(tmp.LIMO.dir,'Betas.mat'));
            load(fullfile(tmp.LIMO.dir,'con_1.mat'));
            load(fullfile(tmp.LIMO.dir,'R2.mat'));
            if j==1
                betas_mean{i} = Betas;
                con_mean{i} = con;
                r2_mean{i} = R2;
            else
                betas_mean{i} = ((j-1)*betas_mean{i} + Betas)./j;
                con_mean{i} = ((j-1)*con_mean{i} + con)./j;
                r2_mean{i} = ((j-1)*r2_mean{i} + R2)./j;
            end
            LIMO = tmp.LIMO;
            save(fullfile(LIMO.dir,'LIMO.mat'),'LIMO')

            i = i + 1;
        end	
        save(fullfile(PATH_TO_DERIV,[model_name '_betas_mean.mat']), 'betas_mean')
        save(fullfile(PATH_TO_DERIV,[model_name '_con_mean.mat']), 'con_mean')
        save(fullfile(PATH_TO_DERIV,[model_name '_r2_mean.mat']), 'r2_mean')
    end

    i=1;
    for subj_name = subj
        subj_name = subj_name{1};
    %     load(fullfile(PATH_TO_DERIV,subj_name, 'eeg', 'GLM_OLS_Time_Channels', 'LIMO.mat'))

    %     Betas = betas_mean{i};
    %     save(fullfile(LIMO.dir,'Betas.mat'),'Betas')
    %     con = con_mean{i};
    %     save(fullfile(LIMO.dir,'con_1.mat'),'con')

        Betas = betas_mean{i};
        save(fullfile(PATH_TO_DERIV,subj_name, 'eeg', 'GLM_OLS_Time_Channels','Betas.mat'),'Betas')
        con = con_mean{i};
        save(fullfile(PATH_TO_DERIV,subj_name, 'eeg', 'GLM_OLS_Time_Channels','con_1.mat'),'con') 
        R2 = r2_mean{i};
        save(fullfile(PATH_TO_DERIV,subj_name, 'eeg', 'GLM_OLS_Time_Channels','R2.mat'),'R2') 

        i = i + 1;
    end

    if save_choice
        save(fullfile(PATH_TO_DERIV, [model_name,'_naive.mat']), 'model')
    end

    % Rename naive GLM folder

    dinfo = dir(fullfile(BIDS_FOLDER,'sub-*'));
    subj = {dinfo.name};
    for subj_name = subj
        subj_name = subj_name{1};
        del_folder = fullfile(BIDS_FOLDER,'derivatives',subj_name,'eeg','GLM_OLS_Time_Channels');
    %     del_folder = fullfile(BIDS_FOLDER,'derivatives',subj_name,'eeg','new_nat_man_naive_GLM_OLS_Time_Channels');
        [root,name,~] = fileparts(del_folder);
        cd(root)
        if exist(del_folder,'dir')
            new_name = [model_name '_naive_GLM_OLS_Time_Channels'];
    %         new_name = 'GLM_OLS_Time_Channels';
            if exist(fullfile(root,new_name),'dir')
                rmdir(new_name,'s')
            end
            movefile(name,new_name);
    %         copyfile(name,new_name);
            cd(fullfile(BIDS_FOLDER,'derivatives',subj_name,'eeg',new_name))
            if exist('LIMO.mat','file')
                load('LIMO.mat')
                LIMO.dir = pwd;
                save('LIMO.mat','LIMO')
            end
        end
    end
end

%% Compute corrected R2 (cf. paper figure 7)
dinfo = dir(fullfile(PATH_TO_DERIV,'sub-*'));
% model_names = {'model_psycho', 'model_image', 'model_psycho_image'}; %Note: a corrected r2 for categorical model is useless
model_names = split(config.model_names,',');
for k = 1:length(model_names)
    model_name = model_names{k};
    absolute_R2 = [];
    naive_R2 = [];
    corrected_R2 = [];
    for i = numel( dinfo ):-1:1
        if i >= 10
            subfolder = ['sub-0' num2str(i)];
        else
            subfolder = ['sub-00' num2str(i)];
        end
        load(fullfile(PATH_TO_DERIV,subfolder,'eeg',[model_name '_GLM_OLS_Time_Channels'],'R2.mat'))
        absolute_R2(:,:,i) = R2(:,:,1);
        if k > 1 %Note: categorical model does not have a naive version
            load(fullfile(PATH_TO_DERIV,subfolder,'eeg',[model_name '_naive_GLM_OLS_Time_Channels'],'R2.mat'))
            naive_R2(:,:,i) = R2(:,:,1);
            corrected_R2(:,:,i) = absolute_R2(:,:,i) - naive_R2(:,:,i);
        end
    end
    if save_choice
        save(fullfile(PATH_TO_DERIV,[model_name '_absolute_R2.mat']),'absolute_R2')
        save(fullfile(PATH_TO_DERIV,[model_name '_naive_R2.mat']),'naive_R2')
        save(fullfile(PATH_TO_DERIV,[model_name '_corrected_R2.mat']),'corrected_R2')
    end
end

%% Find clusters of significant contrast and R2 (2nd level analysis)

% modify Beta_files_GLM_OLS_Time_Channels.txt and con_1_files_GLM_OLS_Time_Channels files to access the desired models
filePath = fullfile(PATH_TO_ROOT,'Beta_files_GLM_OLS_Time_Channels.txt');
fileContentsBetas = fileread(filePath);
filePath = fullfile(PATH_TO_ROOT,'con_1_files_GLM_OLS_Time_Channels.txt');
fileContentsCon= fileread(filePath);

model_names = split(config.model_names,',');

for k = 1:length(model_names)
    model_name = model_names{k};
    % Replace the target string with the new string (Beta files)
    newContents = strrep(fileContentsBetas, 'GLM_OLS_Time_Channels', sprintf("%s_GLM_OLS_Time_Channels",model_name));
    newFilePath = fullfile(PATH_TO_ROOT,['Beta_files_GLM_OLS_Time_Channels_' model_name '.txt']);
    % Write the new contents back to the file
    fid = fopen(newFilePath, 'w');
    if fid == -1
        error('Cannot open file for writing: %s', newFilePath);
    end
%     fwrite(fid, newContents);
    fprintf(fid, '%s', newContents);
    fclose(fid);
    
    % same for con_1
    newContents = strrep(fileContentsCon, 'GLM_OLS_Time_Channels', sprintf("%s_GLM_OLS_Time_Channels",model_name));
    newFilePath = fullfile(PATH_TO_ROOT,['con_1_files_GLM_OLS_Time_Channels_' model_name '.txt']);
    % Write the new contents back to the file
    fid = fopen(newFilePath, 'w');
    if fid == -1
        error('Cannot open file for writing: %s', newFilePath);
    end
    fwrite(fid, newContents);
    fclose(fid);

    cd(PATH_TO_ROOT)
    load(fullfile(PATH_TO_DERIV,[model_name '.mat']))
    expected_chanlocs = limo_avg_expected_chanlocs(PATH_TO_DERIV, model.defaults, model_name);

    % Categorical clusters
    my_param = 'con_1';
    LIMOfiles = fullfile(PATH_TO_ROOT,sprintf('%s_files_GLM_OLS_Time_Channels_%s.txt',my_param,model_name));
    if ~exist(fullfile(PATH_TO_ROOT,[model_name '_' my_param]),'dir')
        mkdir(fullfile(PATH_TO_ROOT,[model_name '_' my_param]))
    end
    cd(fullfile(PATH_TO_ROOT,[model_name '_' my_param]))
    LIMOPath = limo_random_select('one sample t-test',expected_chanlocs,'LIMOfiles',... 
        LIMOfiles,'analysis_type','Full scalp analysis',...
        'type','Channels','nboot',100,'tfce',1,'skip design check','yes');

    p = 0.05;
    MCC = 3; % TFCE
    load('LIMO.mat')
    [~, mask, ~] = limo_stat_values('one_sample_ttest_parameter_1.mat',p,MCC,LIMO);
    if save_choice
        save(fullfile(PATH_TO_DERIV,['mask_' model_name '.mat']), 'mask')
    end
    
    % R2 (explained variance) clusters
    my_param = 'R2';
    if ~exist(fullfile(PATH_TO_ROOT,[model_name '_' my_param]),'dir')
        mkdir(fullfile(PATH_TO_ROOT,[model_name '_' my_param]))
    end
    cd(fullfile(PATH_TO_ROOT,[model_name '_' my_param]))
    LIMO.dir = pwd;
    LIMO.data.data_dir = pwd;
    LIMO.design.method = 'Trimmed Mean';
    LIMO.design.tfce = 0;
    LIMO.design.bootstrap = 100;
    
    load(fullfile(PATH_TO_DERIV,[model_name '_absolute_R2.mat'])) % load the corresponding R2 values
    if k > 1 % no naive model for categorical model
        load(fullfile(PATH_TO_DERIV,[model_name '_naive_R2.mat'])) % load the corresponding R2 values
        limo_random_robust(3,absolute_R2,naive_R2,1,LIMO) % perform the t-test
    else
        limo_random_robust(1,absolute_R2,1,LIMO) % perform the t-test
    end

    p = 0.05;
    MCC = 2; % MCC
    load('LIMO.mat')
    if k > 1 % no naive model for categorical model
        [~, mask, ~] = limo_stat_values('paired_samples_ttest_parameter_1.mat',p,MCC,LIMO);
    else
        [~, mask, ~] = limo_stat_values('one_sample_ttest_parameter_1.mat',p,MCC,LIMO);
    end
    if save_choice
        save(fullfile(PATH_TO_DERIV,['mask_R2_' model_name '.mat']), 'mask')
    end
   
end

%% Boxplot the R2 of each model within specific clusters
% Note: here is an example for the clusters of significant contrast found within the psycho-image
% model. Modify it following your needs.

% /!\ IMPORTANT /!\ "psycho-image" - "image" = "psycho" (cf. paper figure 9)

model_names = split(config.model_names,','); %Note: a corrected r2 for categorical model is useless

% Load the targeted mask
my_param = 'con_1';
model_name = model_names{end};
load(fullfile(PATH_TO_DERIV,['mask_' model_name '.mat'])) % load the mask

% TmR2 = cell(length(model_names),1);
TmR2 = cell(0);
load(fullfile(PATH_TO_DERIV,[model_names{end} '_absolute_R2.mat'])) % load R2 values of the complete model
R2_all_cov = absolute_R2;
for k = 1:length(model_names)-1
    load(fullfile(PATH_TO_DERIV,[model_names{k} '_absolute_R2.mat'])) % load the corresponding R2 values
    if k > 1 % no need to subtract for categorical model
        R2 = R2_all_cov - absolute_R2;
        TmR2{end+1} = limo_trimmed_mean(R2,20,0.05);
    else
        TmR2{end+1} = limo_trimmed_mean(absolute_R2,20,0.05);
    end
    tmp = TmR2{end}(:,:,2);
    tmp(~mask) = nan;
    TmR2{end} = tmp(:);
    TmR2{end} = TmR2{end}(~isnan(TmR2{end}));
end

load(fullfile(PATH_TO_DERIV,[model_names{end} '_naive_R2.mat'])) % load R2 values of the complete model
R2_all_cov = naive_R2;
for k = length(model_names)+2:2*length(model_names)-1
    model_name = model_names{k-length(model_names)};
    load(fullfile(PATH_TO_DERIV,[model_name '_naive_R2.mat'])) % load the corresponding R2 values
    R2 = R2_all_cov - naive_R2;
    TmR2{end+1} = limo_trimmed_mean(R2,20,0.05);
    tmp = TmR2{end}(:,:,2);
    tmp(~mask) = nan;
    TmR2{end} = tmp(:);
    TmR2{end} = TmR2{end}(~isnan(TmR2{end}));
end

target_group = []; % Initialize an empty array
for i = 1:length(TmR2)
    target_group = [target_group, TmR2{i}]; % Concatenate each element
end
isout = isoutlier(target_group,'quartiles');
xClean = target_group;
xClean(isout) = NaN;
figure;boxplot(xClean)
xticks([1,2,3,4,5])
xticklabels({'categorial', 'image', 'psycho', 'naive-image', 'naive-psycho'})
fontSize = 10;
title(sprintf("Explained variance by model\n(R^2 values)"),'FontSize',fontSize)
[est,HDI]=data_plot(xClean,'estimator','trimmed mean'); % test with estimator 
xticks([1,2.25,3.5,4.75,6])
xticklabels({'categorial', 'image', 'psycho', 'naive-image', 'naive-psycho'})
fontSize = 10;
title(sprintf("Explained variance by model\n(R^2 values)"),'FontSize',fontSize)

%% Quantify the separability (cf. paper figure 9)
R2_losses = cell(0);
% Here is the example of the R2 loss between categories and psycho
% covariates. Adapt it to your data.

% 1) "computed" psycho effect
load(fullfile(PATH_TO_DERIV,[model_names{end} '_absolute_R2.mat'])) % load R2 values of the psycho-image model
R2_all_cov = absolute_R2;
load(fullfile(PATH_TO_DERIV,[model_names{3} '_absolute_R2.mat'])) % load R2 values of the image model
R2_computed_psycho = R2_all_cov - absolute_R2;

% 2) "computed" categorical effect
load(fullfile(PATH_TO_DERIV,[model_names{2} '_absolute_R2.mat'])) % load R2 values of the psycho model
R2_computed_cat = absolute_R2 - R2_computed_psycho;

% 3) R2 loss
load(fullfile(PATH_TO_DERIV,[model_names{1} '_absolute_R2.mat'])) % load R2 values of the categorical model
R2_loss = absolute_R2 - R2_computed_cat;

R2_losses{end+1} = R2_loss;

% Same for the R2 loss between categories and image covariates
% 1) "computed" image effect
load(fullfile(PATH_TO_DERIV,[model_names{end} '_absolute_R2.mat'])) % load R2 values of the psycho-image model
R2_all_cov = absolute_R2;
load(fullfile(PATH_TO_DERIV,[model_names{2} '_absolute_R2.mat'])) % load R2 values of the psycho model
R2_computed_image = R2_all_cov - absolute_R2;

% 2) "computed" categorical effect
load(fullfile(PATH_TO_DERIV,[model_names{3} '_absolute_R2.mat'])) % load R2 values of the image model
R2_computed_cat = absolute_R2 - R2_computed_image;

% 3) R2 loss
load(fullfile(PATH_TO_DERIV,[model_names{1} '_absolute_R2.mat'])) % load R2 values of the categorical model
R2_loss = absolute_R2 - R2_computed_cat;

R2_losses{end+1} = R2_loss;

end