diff --git a/createTestData/plot_all_test_data_profit.m b/createTestData/plot_all_test_data_profit.m index ac80b7430f4fcc8c610c89e8c1aab8eca12dbbc1..c8f8f0844e52f708224140071c0c8f0cb40e2d96 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