From 54a195ad1b01b948d8540bd2c3a6f33be844cfc2 Mon Sep 17 00:00:00 2001
From: Tamas Borbath <tamas.borbath@tuebingen.mpg.de>
Date: Wed, 16 Dec 2020 19:24:37 +0100
Subject: [PATCH] Cleaned up a bit and adapted plot_all_test_data_profit code
to new Bland-Altman plots. However the code is still not clean
---
createTestData/plot_all_test_data_profit.m | 241 ++++++++++++++++-----
1 file changed, 183 insertions(+), 58 deletions(-)
diff --git a/createTestData/plot_all_test_data_profit.m b/createTestData/plot_all_test_data_profit.m
index ac80b74..c8f8f08 100644
--- a/createTestData/plot_all_test_data_profit.m
+++ b/createTestData/plot_all_test_data_profit.m
@@ -27,12 +27,8 @@ exportFolders = {'0 Water References',...
pathNameExport = 'ProFit test data';
pathBaseExportFiles = pathToDataFolder(pathNameExport, exportFolders);
-% pathBaseExportFiles = pathBaseExportFiles(1:end-1);
-% pathBaseExportFiles = [pathBaseExportFiles ' - Weighted Cost Fun R3 0.25scale m_15 trunc\'];
-% pathBaseExportFiles = [pathBaseExportFiles ' - Weighted Cost Fun R3 0.25scale m_15 sinebell\'];
-% pathBaseExportFiles = [pathBaseExportFiles ' - Weighted Cost Fun R3 0.25scale m_15 matched\'];
-% pathBaseExportFiles = [pathBaseExportFiles ' - R1 0.25scale m_15 trunc\'];
-% pathBaseExportFiles = [pathBaseExportFiles ' - R2 0.25scale m_15 trunc\'];
+pathBaseExportFiles = pathBaseExportFiles(1:end-1);
+pathBaseExportFiles = [pathBaseExportFiles ' - R3\'];
dataExportPathBaseRef = strcat(pathBaseExportFiles, exportFolders{2}, '\');
dataExportPathBase = strcat(pathBaseExportFiles, exportFolders{reconOption+1}, '\');
@@ -105,7 +101,7 @@ if fitLCModel
else
SoftwareName = 'ProFit';
end
-[metConcAllSubjLCModel, metConcAllSubjRefLCModel, metConcAllSubjNormLCModel] = ...
+[metConcAllSubjLCModel, metConcAllSubjRefLCModel, metConcAllSubjNormLCModel, MC_H2O_ConcAllSubjLCModel, MC_H2O_ConcAllSubjRefLCModel] = ...
getAllMetaboliteConc(dataExportPathBase, dataExportPathBaseRef, truncSuffix, subjects, TEs, ...
namingSuffixesAve, namingSuffixesCoil, fitLCModel, reconOption, ...
numberOfMet, metabolitesLCModel, metabolitesProFit, metabolitesLabels, indecesMainMet, ...
@@ -117,7 +113,7 @@ if fitLCModel
else
SoftwareName = 'ProFit';
end
-[metConcAllSubjProFit, metConcAllSubjRefProFit, metConcAllSubjNormProFit] = ...
+[metConcAllSubjProFit, metConcAllSubjRefProFit, metConcAllSubjNormProFit, MC_H2O_ConcAllSubjProFit, MC_H2O_ConcAllSubjRefProFit] = ...
getAllMetaboliteConc(dataExportPathBase, dataExportPathBaseRef, truncSuffix, subjects, TEs, ...
namingSuffixesAve, namingSuffixesCoil, fitLCModel, reconOption, ...
numberOfMet, metabolitesLCModel, metabolitesProFit, metabolitesLabels, indecesMainMet, ...
@@ -136,15 +132,17 @@ metConcAllSubjLCModelScaled = metConcAllSubjLCModel;
metConcAllSubjRefLCModelScaled = metConcAllSubjRefLCModel;
% BA plot paramters
-tit = 'Fitting Repeatability '; % figure title
-subsets = {'32 ave.', '64 ave.'};
-% subsets = {'32 ave.1','32 ave.2','32 ave.3', '64 ave.1', '64 ave.2', '64 ave.3'};
+tit = 'Repeatability '; % figure title
+subsets = {'32 ave.1','32 ave.2','32 ave.3', '64 ave.1', '64 ave.2', '64 ave.3'};
+subsets_refit = {'32 ave.', '64 ave.'};
softwareNames = {'ProFit', 'LCModel'};
gnames = {subsets, softwareNames}; % names of groups in data {dimension 1 and 2}
label = {'Ref. Fits','Subset Fits','ratio'}; % Names of data sets
+label_refit = {'c^{fits1}_{i,k}','c^{fits2}_{i,k}','arb.u.'}; % Names of data sets
label_SW_comp = {'LCModel','ProFit','arb.u.'}; % Names of data sets
corrinfo = {'n','SSE','r2','eq'}; % stats to display of correlation scatter plot
BAinfo = {'RPC','ks'}; % stats to display on Bland-ALtman plot
+BAinfo_refit = {''}; % stats to display on Bland-ALtman plot
limits = 'auto'; % how to set the axes limits
if 1 % colors for the data sets may be set as:
colors = 'br'; % character codes
@@ -152,61 +150,147 @@ else
colors = [0 0 1;... % or RGB triplets
1 0 0];
end
+
+paramsBA.data1TreatmentMode = 'Compare';
+paramsBA.diffValueMode = 'percent';
+
numOfSimulations = length(namingSuffixesAve)*length(namingSuffixesCoil);
-rpc_table = cell(numberOfMet+3,3);
+%there seems to be a fit error in the fitting of even the MC water data.
+% BlandAltman(MC_H2O_ConcAllSubjLCModel, MC_H2O_ConcAllSubjProFit*1e-6, label_SW_comp,[tit 'Water'],subsets,...
+% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'showFitCI',' on',...
+% 'baStatsMode','Gaussian','forceZeroIntercept','on', 'diffValueMode', 'percent');
+%
+% BlandAltman(MC_H2O_ConcAllSubjRefLCModel, MC_H2O_ConcAllSubjRefProFit*1e-6, label_SW_comp,[tit 'Water'],subsets,...
+% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'showFitCI',' on',...
+% 'baStatsMode','Gaussian','forceZeroIntercept','on', 'diffValueMode', 'percent');
+
+rpc_table = cell(numberOfMet+5,5);
rpc_table{1,2} = softwareNames{1};
-rpc_table{1,3} = softwareNames{2};
-rpc_table(2:end-2,1) = metabolitesToDisplay;
+rpc_table{1,4} = softwareNames{2};
+rpc_table{2,2} = subsets_refit{1};
+rpc_table{2,3} = subsets_refit{2};
+rpc_table{2,4} = subsets_refit{1};
+rpc_table{2,5} = subsets_refit{2};
+rpc_table(3:end-3,1) = metabolitesToDisplay;
+rpc_table{end-1,1} = 'Mean Main';
rpc_table{end,1} = 'Mean';
+
+%title('Reproducibility');
+subplotsX = 4;
+subplotsY = 4;
+totalSubplots = subplotsX * subplotsY;
+rpcSum = zeros(1,numberOfMet);
+diffMedian = zeros(1,numberOfMet);
for indexMet = 1: numberOfMet
if indexMet == metaboliteRef
continue;
end
+ if mod(indexMet,totalSubplots/2) == 1
+ figure;
+ end
+ runningIndexPlot = (indexMet-1)+1+floor((indexMet-1)/subplotsY)*subplotsY-floor((indexMet-1)/totalSubplots*2)*totalSubplots;
+
%ProFit Concentrations
- metConcAllSubjProFitPerCr_ = metConcAllSubjProFitScaled(:, :, indexMet);
- metConcAllSubjProFitPerCr_ = reshape(metConcAllSubjProFitPerCr_,11*3,2);%group same SNR
- metConcAllSubjRefProFitPerCr_ = metConcAllSubjRefProFitScaled(:, indexMet);
- metConcAllSubjRefProFitPerCr_ = repmat(metConcAllSubjRefProFitPerCr_,1,6);
- metConcAllSubjRefProFitPerCr_ = reshape(metConcAllSubjRefProFitPerCr_,11*3,2);%group same SNR
+ metConcAllSubjProFitScaled_ = metConcAllSubjProFitScaled(:, :, indexMet);
+ metConcAllSubjProFitScaled_ = reshape(metConcAllSubjProFitScaled_,11*3,2);%group same SNR
+ metConcAllSubjRefProFitScaled_ = metConcAllSubjRefProFitScaled(:, indexMet);
+ metConcAllSubjRefProFitScaled_ = repmat(metConcAllSubjRefProFitScaled_,1,6);
+ metConcAllSubjRefProFitScaled_ = reshape(metConcAllSubjRefProFitScaled_,11*3,2);%group same SNR
%LCModel Concentrations
- metConcAllSubjLCModelPerCr_ = metConcAllSubjLCModelScaled(:, :, indexMet);
- metConcAllSubjLCModelPerCr_ = reshape(metConcAllSubjLCModelPerCr_,11*3,2);%group same SNR
- metConcAllSubjRefLCModelPerCr_ = metConcAllSubjRefLCModelScaled(:, indexMet);
- metConcAllSubjRefLCModelPerCr_ = repmat(metConcAllSubjRefLCModelPerCr_,1,6);
- metConcAllSubjRefLCModelPerCr_ = reshape(metConcAllSubjRefLCModelPerCr_,11*3,2);%group same SNR
+ metConcAllSubjLCModelScaled_ = metConcAllSubjLCModelScaled(:, :, indexMet);
+ metConcAllSubjLCModelScaled_ = reshape(metConcAllSubjLCModelScaled_,11*3,2);%group same SNR
+ metConcAllSubjRefLCModelScaled_ = metConcAllSubjRefLCModelScaled(:, indexMet);
+ metConcAllSubjRefLCModelScaled_ = repmat(metConcAllSubjRefLCModelScaled_,1,6);
+ metConcAllSubjRefLCModelScaled_ = reshape(metConcAllSubjRefLCModelScaled_,11*3,2);%group same SNR
%ProFit&LCModel Concentrations
- metConcAllSubjPerCr_ = cat(3, metConcAllSubjProFitPerCr_, metConcAllSubjLCModelPerCr_);
- metConcAllSubjRefPerCr_ = cat(3, metConcAllSubjRefProFitPerCr_, metConcAllSubjRefLCModelPerCr_);
- %ref Fits vs subfits
-% BlandAltman(metConcAllSubjRefPerCr_, metConcAllSubjPerCr_, label,[tit metabolitesLCModel{indexMet}],gnames,...
+ metConcAllSubjScaled_ = cat(3, metConcAllSubjProFitScaled_, metConcAllSubjLCModelScaled_);
+ metConcAllSubjRefScaled_ = cat(3, metConcAllSubjRefProFitScaled_, metConcAllSubjRefLCModelScaled_);
+%% ref Fits vs subfits: Code not used anymore in this style
+% BlandAltman(metConcAllSubjRefScaled_, metConcAllSubjScaled_, label,[tit metabolitesLCModel{indexMet}],gnames,...
% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'colors',colors, 'showFitCI',' on',...
% 'baStatsMode','non-parametric','forceZeroIntercept','on');
- %LCModel vs ProFit
- BlandAltman(metConcAllSubjLCModelPerCr_, metConcAllSubjProFitPerCr_, label_SW_comp,[tit metabolitesLCModel{indexMet}],subsets,...
- 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'colors',colors, 'showFitCI',' on',...
- 'baStatsMode','non-parametric','forceZeroIntercept','on', 'diffValueMode', 'percent');
- %ProFit ref Fits vs subfits
-% rpcProFit = BlandAltman(metConcAllSubjRefProFitPerCr_, metConcAllSubjProFitPerCr_, label,[tit softwareNames{1} ' ' metabolitesLCModel{indexMet}],subsets,...
+% % Reference fits vs all per Software
+% %ProFit ref Fits vs subfits
+% rpcProFit = BlandAltman(metConcAllSubjRefProFitScaled_, metConcAllSubjProFitScaled_, label,[tit softwareNames{1} ' ' metabolitesLCModel{indexMet}],subsets,...
% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'colors',[0.75 0.5 0; 0.5 0.2 0], 'showFitCI',' on',...
% 'baStatsMode','Gaussian','forceZeroIntercept','on', 'diffValueMode', 'percent');
% %LCModel ref Fits vs subfits
-% rpcLCModel = BlandAltman(metConcAllSubjRefLCModelPerCr_, metConcAllSubjLCModelPerCr_, label,[tit softwareNames{2} ' ' metabolitesLCModel{indexMet}],subsets,...
+% rpcLCModel = BlandAltman(metConcAllSubjRefLCModelScaled_, metConcAllSubjLCModelScaled_, label,[tit softwareNames{2} ' ' metabolitesLCModel{indexMet}],subsets,...
% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'colors',[0, 0.5, 0.75; 0, 0.25 0.5], 'showFitCI',' on',...
% 'baStatsMode','Gaussian','forceZeroIntercept','on', 'diffValueMode', 'percent');
-% rpc_table{indexMet+1,2} = rpcProFit;
-% rpc_table{indexMet+1,3} = rpcLCModel;
+% %rpc_table{indexMet+1,2} = rpcProFit;
+% %rpc_table{indexMet+1,3} = rpcLCModel;
+%% LCModel vs ProFit
+% BlandAltman(metConcAllSubjLCModelScaled_, metConcAllSubjProFitScaled_, label_SW_comp,[tit metabolitesLCModel{indexMet}],subsets,...
+% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'colors',colors, 'showFitCI',' on',...
+% 'baStatsMode','non-parametric','forceZeroIntercept','on', 'diffValueMode', 'percent');
+% [rpc, ~, stat] = BlandAltman(metConcAllSubjRefLCModelScaled_, metConcAllSubjRefProFitScaled_, label_SW_comp,[tit metabolitesLCModel{indexMet}],subsets,...
+% 'corrInfo',corrinfo,'baInfo',BAinfo,'axesLimits',limits,'colors',colors, 'showFitCI',' on',...
+% 'baStatsMode','non-parametric','forceZeroIntercept','on', 'diffValueMode', 'percent');
+% rpcs(indexMet) = rpc;
+% diffMedian(indexMet) = stat.differenceMedian;
+%% Subfits vs subfits per Software
+ %ProFit ref Fits vs subfits
+ rpc_ProFit_32_ave = calculate_rpc_percent(metConcAllSubjProFitScaled_(1:11,1), metConcAllSubjProFitScaled_(12:22,1));
+ rpc_ProFit_64_ave = calculate_rpc_percent(metConcAllSubjProFitScaled_(1:11,2), metConcAllSubjProFitScaled_(12:22,2));
+ legendText = {['{\color{blue}\bf32-ave: RPC=' mynum2str(rpc_ProFit_32_ave,2) '%}'];['{\color{red}\bf64-ave: RPC=' mynum2str(rpc_ProFit_64_ave,2) '%}']};
+
+ sProFit = subplot(subplotsX,subplotsY,runningIndexPlot);
+ [rpcProFit, fig, stats] = BlandAltmanOnly(sProFit, metConcAllSubjProFitScaled_(1:11,:), metConcAllSubjProFitScaled_(12:22,:),...
+ label_refit,[metabolitesToDisplay{indexMet}],legendText,...
+ 'corrInfo',corrinfo,'baInfo',BAinfo_refit,'axesLimits',limits,'colors','br', 'showFitCI',' on',...
+ 'baStatsMode','Gaussian','forceZeroIntercept','on', 'diffValueMode', 'percent');
+
+ %LCModel ref Fits vs subfits
+ rpc_LCModel_32_ave = calculate_rpc_percent(metConcAllSubjLCModelScaled_(1:11,1), metConcAllSubjLCModelScaled_(12:22,1));
+ rpc_LCModel_64_ave = calculate_rpc_percent(metConcAllSubjLCModelScaled_(1:11,2), metConcAllSubjLCModelScaled_(12:22,2));
+ legendText = {['{\color{blue}\bf32-ave: RPC=' mynum2str(rpc_LCModel_32_ave,2) '%}'];['{\color{red}\bf64-ave: RPC=' mynum2str(rpc_LCModel_64_ave,2) '%}']};
+
+ sLCModel = subplot(subplotsX,subplotsY,runningIndexPlot+subplotsY);
+ [rpcLCModel, fig, stats] = BlandAltmanOnly(sLCModel, metConcAllSubjRefLCModelScaled_(1:11,:), metConcAllSubjLCModelScaled_(12:22,:), ...
+ label_refit,[metabolitesToDisplay{indexMet}], legendText,...
+ 'corrInfo',corrinfo,'baInfo',BAinfo_refit,'axesLimits',limits,'colors','br', 'showFitCI',' on',...
+ 'baStatsMode','Gaussian','forceZeroIntercept','on', 'diffValueMode', 'percent');
+
+ if mod(indexMet,subplotsY)==1
+ pos1 = sProFit.Position;
+ annotation('textarrow', [pos1(1), 0], [pos1(2)+pos1(4)/2, 0],'String',softwareNames{1}, ...
+ 'HeadStyle','none','LineStyle', 'none', 'TextRotation',90,'units','normalized', 'FontSize', 16, ...
+ 'FontWeight','bold','HorizontalAlignment','center', 'VerticalAlignment','middle');
+ pos2 = sLCModel.Position;
+ annotation('textarrow', [pos1(1), 0], [pos2(2)+pos2(4)/2, 0],'String',softwareNames{2}, ...
+ 'HeadStyle','none','LineStyle', 'none', 'TextRotation',90,'units','normalized', 'FontSize', 16, ...
+ 'FontWeight','bold','HorizontalAlignment','center', 'VerticalAlignment','middle');
+ end
+ rpc_table{indexMet+2,2} = rpc_ProFit_32_ave;
+ rpc_table{indexMet+2,3} = rpc_ProFit_64_ave;
+ rpc_table{indexMet+2,4} = rpc_LCModel_32_ave;
+ rpc_table{indexMet+2,5} = rpc_LCModel_64_ave;
end
-rpc_meanProFit = mean(cell2mat(rpc_table(2:end-2,2)));
-rpc_meanLCModel = mean(cell2mat(rpc_table(2:end-2,3)));
-rpc_table{end,2} = rpc_meanProFit;
-rpc_table{end,3} = rpc_meanLCModel;
+rpc_values = cell2mat(rpc_table(3:end-2,2:5));
+rpc_meanProFit_32ave_main = mean(rpc_values(indecesMainMet,1));
+rpc_meanProFit_64ave_main = mean(rpc_values(indecesMainMet,2));
+rpc_meanLCModel_32ave_main = mean(rpc_values(indecesMainMet,3));
+rpc_meanLCModel_64ave_main = mean(rpc_values(indecesMainMet,4));
+rpc_table{end-1,2} = rpc_meanProFit_32ave_main;
+rpc_table{end-1,3} = rpc_meanProFit_64ave_main;
+rpc_table{end-1,4} = rpc_meanLCModel_32ave_main;
+rpc_table{end-1,5} = rpc_meanLCModel_64ave_main;
+rpc_meanProFit_32ave = mean(rpc_values(:,1));
+rpc_meanProFit_64ave = mean(rpc_values(:,2));
+rpc_meanLCModel_32ave = mean(rpc_values(:,3));
+rpc_meanLCModel_64ave = mean(rpc_values(:,4));
+rpc_table{end,2} = rpc_meanProFit_32ave;
+rpc_table{end,3} = rpc_meanProFit_64ave;
+rpc_table{end,4} = rpc_meanLCModel_32ave;
+rpc_table{end,5} = rpc_meanLCModel_64ave;
xlswrite([pathBaseExportFiles, 'InVivoResults_RPC.xlsx'], rpc_table)
-
+%% code will break here. This used to be the old style plotting/results without Bland-Altman plots.
meanMetDeviation = mean(abs(metaboliteConcentrationsAllSubjNorm), 'omitnan');
stdMetDeviation = std(abs(metaboliteConcentrationsAllSubjNorm), 'omitnan');
@@ -265,7 +349,8 @@ end
% metaboliteConcentrations(1) = metaboliteConcentrations(1) / 1e6; % adjust the MMB to be within range
end
-function [metaboliteConcentrationsAllSubj, metaboliteConcentrationsAllSubjRef, metaboliteConcentrationsAllSubjNorm] = ...
+function [metaboliteConcentrationsAllSubj, metaboliteConcentrationsAllSubjRef, metaboliteConcentrationsAllSubjNorm,...
+ MC_H2O_ConcentrationsAllSubj, MC_H2O_ConcentrationsAllSubjRef] = ...
getAllMetaboliteConc(dataExportPathBase, dataExportPathBaseRef, truncSuffix, subjects, TEs, ...
namingSuffixesAve, namingSuffixesCoil, fitLCModel, reconOption, ...
numberOfMet, metabolitesLCModel, metabolitesProFit, metabolitesLabels, indecesMainMet, ...
@@ -308,25 +393,34 @@ meanFQN = 0;
meanFQN2 = 0;
numOfSimulations = length(namingSuffixesAve)*length(namingSuffixesCoil);
metaboliteConcentrationsAllSubjNorm = zeros(numberOfSubjects * numOfSimulations, numberOfMet);
+MC_H2O_ConcentrationsAllSubj = zeros(numberOfSubjects, numOfSimulations);
metaboliteConcentrationsAllSubj = zeros(numberOfSubjects, numOfSimulations, numberOfMet);
metaboliteConcentrationsAllSubjRef = zeros(numberOfSubjects, numberOfMet);
+MC_H2O_ConcentrationsAllSubjRef = zeros(numberOfSubjects, 1);
%% do the actual plotting
for indexSubj = 1:numberOfSubjects
if fitLCModel
load([dataExportPathRef{indexSubj}, 'concentrations' truncSuffix '.mat'], 'concentrations');
+ load([dataExportPathRef{indexSubj}, 'concentrations', truncSuffix '_MC_H2O.mat'], 'concentrations_MC_H2O');
concentrationsRef = concentrations;
+ concentrations_MC_H2ORef = concentrations_MC_H2O;
load([dataExportPath{indexSubj}, 'concentrations' truncSuffix '.mat'], 'concentrations');
+ load([dataExportPath{indexSubj}, 'concentrations', truncSuffix '_MC_H2O.mat'], 'concentrations_MC_H2O');
else
load([dataExportPathRef{indexSubj}, 'concentrations' truncSuffix '_profit.mat'], 'concentrations', 'activeMetabolites', 'FQNs', 'FQNs2');
load([dataExportPathRef{indexSubj}, 'concentrations', truncSuffix '_H2O_profit.mat'], 'concentrations_H2O');
+ load([dataExportPathRef{indexSubj}, 'concentrations', truncSuffix '_MC_H2O_profit.mat'], 'concentrations_MC_H2O');
concentrationsRef = concentrations;
concentrations_H2ORef = cell2mat(concentrations_H2O);
+ concentrations_MC_H2ORef = cell2mat(concentrations_MC_H2O);
activeMetabolitesRef = activeMetabolites;
FQNsRef = FQNs;
FQNsRef2 = FQNs2;
load([dataExportPath{indexSubj}, 'concentrations' truncSuffix '_profit.mat'], 'concentrations', 'activeMetabolites', 'FQNs', 'FQNs2');
load([dataExportPath{indexSubj}, 'concentrations', truncSuffix '_H2O_profit.mat'], 'concentrations_H2O');
+ load([dataExportPath{indexSubj}, 'concentrations', truncSuffix '_MC_H2O_profit.mat'], 'concentrations_MC_H2O');
concentrations_H2O = cell2mat(concentrations_H2O);
+ concentrations_MC_H2O = cell2mat(concentrations_MC_H2O);
end
% load stats (SNR, FWHM, etc)
load([dataExportPathRef{indexSubj}, 'stats' truncSuffix '.mat'], 'stats');
@@ -355,12 +449,16 @@ for indexSubj = 1:numberOfSubjects
currentConcentrationRef = concentrationsRef{indexTE,1, 1};
if fitLCModel
+% MC_H2O_ConcentrationsRef = concentrations_MC_H2ORef{indexTE,1, 1};
+% MC_H2O_ConcentrationsAllSubjRef(indexSubj) = MC_H2O_ConcentrationsRef{2,2};
metaboliteConcentrationsRef = matchFittedMetabolites(fitLCModel, currentConcentrationRef, metabolitesLCModel, [], []);
metaboliteConcentrationsRef(end) = metaboliteConcentrationsRef(end) .* 1e-8;
else
metaboliteConcentrationsRef_notScaled = matchFittedMetabolites(fitLCModel, currentConcentrationRef, metabolitesProFit, proFitIteration, activeMetabolitesRef);
%scale the concentations similarly as LCModel
metaboliteConcentrationsRef = metaboliteConcentrationsRef_notScaled * 40873 / concentrations_H2ORef * 2;
+% MC_H2O_ConcentrationsRef = concentrations_MC_H2ORef * 40873 / concentrations_H2ORef * 2;
+% MC_H2O_ConcentrationsAllSubjRef(indexSubj) = MC_H2O_ConcentrationsRef;
end
if ~fitLCModel
currentFQNRef = FQNsRef{indexTE,1, 1};
@@ -390,17 +488,23 @@ for indexSubj = 1:numberOfSubjects
figure(indexFigure);
hold on
if fitLCModel
+% MC_H2O_Concentrations = concentrations_MC_H2O{indexTE,indexAverageDelete, indexCoilDelete};
+% MC_H2O_ConcentrationsAllSubj(indexSubj, runningIndex) = MC_H2O_Concentrations{2,2};
metaboliteConcentrations = matchFittedMetabolites(fitLCModel, currentConcentration, metabolitesLCModel, [], []);
metaboliteConcentrations(end) = metaboliteConcentrations(end) .* 1e-8;
else
metaboliteConcentrations_notScaled = matchFittedMetabolites(fitLCModel, currentConcentration, metabolitesProFit, proFitIteration, activeMetabolites);
%scale the concentations similarly as LCModel
metaboliteConcentrations = metaboliteConcentrations_notScaled * 40873 / concentrations_H2O(runningIndex) * 2;
+% MC_H2O_Concentrations = concentrations_MC_H2O(runningIndex) * 40873 / concentrations_H2O(runningIndex) * 2;
+% MC_H2O_ConcentrationsAllSubj(indexSubj, runningIndex) = MC_H2O_Concentrations;
end
if runningIndex <= 3 %scale with coresponding averages
metaboliteConcentrations = metaboliteConcentrations * 3;
+ MC_H2O_ConcentrationsAllSubj(indexSubj, runningIndex) = MC_H2O_ConcentrationsAllSubj(indexSubj, runningIndex) * 3;
else
metaboliteConcentrations = metaboliteConcentrations *3/2;
+ MC_H2O_ConcentrationsAllSubj(indexSubj, runningIndex) = MC_H2O_ConcentrationsAllSubj(indexSubj, runningIndex) *3/2;
end
if ~fitLCModel
currentFQN = FQNs{indexTE,indexAverageDelete, indexCoilDelete};
@@ -448,22 +552,43 @@ for indexSubj = 1:numberOfSubjects
mainMetIterator = iterator(indecesMainMet);
title(SoftwareName)
- % if fitLCModel
- % title([SoftwareName ... ' Subject: ', subjects{indexSubj}, ' TE ', num2str(TEs(indexTE)) ' ms'...
- % ' - mean change (main mets): ', num2str(mean(meanPercentualChange./iterator,'omitnan'),3) '%; '...
- % ' (', num2str(mean(mainMetMeanChange./mainMetIterator,'omitnan'),3) '%); '...
- % num2str(sum(nansIterator)), ' (', num2str(sum(mainMetNaNs)), ') NaNs'...
- % ]);
- % else
- % title([SoftwareName ... ' Subject: ', subjects{indexSubj}, ' TE ', num2str(TEs(indexTE)) ' ms'...
- % ' - mean change (main mets): ', num2str(mean(meanPercentualChange./iterator,'omitnan'),3) '%; '...
- % ' (', num2str(mean(mainMetMeanChange./mainMetIterator,'omitnan'),3) '%); '...
- % num2str(sum(nansIterator)), ' (', num2str(sum(mainMetNaNs)), ') NaNs'...
- % ', FQN_{Ref}:' num2str(FQN_Ref/numberOfSubjects,3), ', FQN_{mean}:' num2str(meanFQN/max(iterator),3)...
- % ', FQN_{Ref}:' num2str(FQN_Ref2/numberOfSubjects,3), ', FQN_{mean}:' num2str(meanFQN2/max(iterator),3)...
- % ]);
- % end
+ if fitLCModel
+ title([SoftwareName ... ' Subject: ', subjects{indexSubj}, ' TE ', num2str(TEs(indexTE)) ' ms'...
+ ' - mean change (main mets): ', num2str(mean(meanPercentualChange./iterator,'omitnan'),3) '%; '...
+ ' (', num2str(mean(mainMetMeanChange./mainMetIterator,'omitnan'),3) '%); '...
+ num2str(sum(nansIterator)), ' (', num2str(sum(mainMetNaNs)), ') NaNs'...
+ ]);
+ else
+ title([SoftwareName ... ' Subject: ', subjects{indexSubj}, ' TE ', num2str(TEs(indexTE)) ' ms'...
+ ' - mean change (main mets): ', num2str(mean(meanPercentualChange./iterator,'omitnan'),3) '%; '...
+ ' (', num2str(mean(mainMetMeanChange./mainMetIterator,'omitnan'),3) '%); '...
+ num2str(sum(nansIterator)), ' (', num2str(sum(mainMetNaNs)), ') NaNs'...
+ ', FQN_{Ref}:' num2str(FQN_Ref/numberOfSubjects,3), ', FQN_{mean}:' num2str(meanFQN/max(iterator),3)...
+ ', FQN_{Ref}:' num2str(FQN_Ref2/numberOfSubjects,3), ', FQN_{mean}:' num2str(meanFQN2/max(iterator),3)...
+ ]);
+ end
end
end
+end
+
+function rpc = calculate_rpc_percent(dataS1, dataS2)
+ data.set1 = dataS1;
+ data.set2 = dataS2;
+ s = size(data.set1);
+ if ~isequal(s,size(data.set2))
+ error('data1 and data2 must have the same size');
+ end
+
+ data.set1 = reshape(data.set1, [numel(data.set1),1]);
+ data.set2 = reshape(data.set2, [numel(data.set2),1]);
+ data.mask = isfinite(data.set1) & isnumeric(data.set1) & isfinite(data.set2) & isnumeric(data.set2);
+ data.maskedSet1 = data.set1(data.mask);
+ data.maskedSet2 = data.set2(data.mask);
+
+ data.maskedBaRefData = mean([data.maskedSet1,data.maskedSet2],2);
+ data.maskedDifferences = (data.maskedSet2-data.maskedSet1) ./ data.maskedBaRefData*100;
+
+ differenceSTD = std(data.maskedDifferences, 'omitnan');
+ rpc = 1.96*differenceSTD;
end
\ No newline at end of file
--
GitLab