NMM_QCfit.m

%% CODE TO PERFORM ADDITIONAL QC MRS ANALYSIS REQUIRED IN 1ST ROUND REV

% This code extracts QC NMM fit metrics (R2 and correlations between
% observed and modeled data)

% It needs
% a "your_paths.txt" file with all the info on your paths. An example is
% provided on GIT + the DCM estimated.
%-----------------------------------------------------------------------
% Written by I. Paparella
% Cyclotron Research Centre, University of Liege, Belgium
% 07th November 2024

%% Path info
paths_txt = fileread('your_paths.txt') ;
paths_txt = char(strsplit(paths_txt, '\n')) ;

pathtodata = (strtrim(paths_txt(2,:)));

%%  Select the subject to analyse (all except 001, 007 and 011 which are males)
cd(pathtodata);
[subdir] = spm_select(Inf,'dir','Select subject_ID directory');
IDs = str2num(subdir(:,end-2:end));

%% Initalize variables
%MRS_QC_metrics = zeros(size(subdir,1),(length(Metrics)*length(Metabs)));
R2 = [];
r = [];
modes= [1 2 3];
%% Loop across subjects to analyse and get subjective parameters
for isub = 1:size(subdir,1)
    
    % Extract path info per subject
    subject_ID = sprintf('sub-%03d', IDs(isub));
    subject_DCMest = fullfile (pathtodata, subject_ID, 'TMS_EEG');
    
    % Load DCM file
    DCMfile = spm_select('List',subject_DCMest,'DCM_DefaultCodes_.*.mat$');
    SubDCM=load(fullfile(subject_DCMest, '/', DCMfile(1,:)));
    
% Compute R2    
    for m = 1:length(modes)
        % compute and store R2 for ERP
        MSE_1 = sum(SubDCM.DCM.R{1}(:,m).^2); % for condition 1
        SS_1  = sum((SubDCM.DCM.H{1}(:,m)+SubDCM.DCM.R{1}(:,m)).^2);
        R2(isub,m,1) = 1-(MSE_1/SS_1);
        clear MSE_1 SS_1
    end
    
% Compute correlations between observed and modeled data    
    for m = 1:length(modes)
        % correlation
        y1 = SubDCM.DCM.H{1}(:,m)+SubDCM.DCM.R{1}(:,m);        
        yhat1 = SubDCM.DCM.H{1}(:,m);
        r(isub,m,1) = corr(y1(:,1),yhat1(:,1),'type','Pearson');
    end    
end

% Convert to table and save
r=array2table(r);
R2=array2table(R2);
IDs_table = table(IDs);

QC_metrics_full = [IDs_table r R2];

QC_metrics_full.Properties.VariableNames = {'Sbj_number'; 'Mode1_r'; 'Mode2_r'; 'Mode3_r';
    'Mode1_R2'; 'Mode2_R2'; 'Mode3_R2'};
writetable(QC_metrics_full,fullfile(pathtodata,'DCM_QC_1stRev_AllFemale.csv'),'WriteRowNames',true);