SubcorticalDevianceDetectionHumans/Deviance Detection_Humans/DD_PlottingConSCurve.m

%% define some variables

% general parameters
plotRespI = 0; % plot responses?
plotVarI = 1; % plot variables?
saveI = 0; % should the figures be saved? (0: no, 1: yes)
fixedY = 0; % should the y-axis be fixed?
y_lim = [-0.25,0.35]; % limits for fixed y-axis

% analysis parameters
modelTy = 1; % model type to be plotted (1: linear regression model, 2: flexible model)
filt = 2; % filter (1: none, 2: 60 Hz, 3: 70 Hz)
stPaPlot = [1,4]; % statistical parameters to be considered for plotting

% cosmetics
figWid = 1000; % width of figure window
figHei = 600; % height of figure window
winCol = [0.9290,0.6940,0.1250;0.4940,0.1840,0.5560]; % colors for analysis windows

%% load and concatenate data

% load and rename av data
load("DD_avData_Oddball_25ms_SbD.mat")
dataAv_25 = conSAv;
dataGrAv_25 = conSGrAv;
dataSe_25 = conSSe;
load("DD_avData_Oddball_30ms_SbD.mat")
dataAv_30 = conSAv;
dataGrAv_30 = conSGrAv;
dataSe_30 = conSSe;
load("DD_avData_Oddball_50ms_SbD.mat")
dataAv_50 = conSAv;
dataGrAv_50 = conSGrAv;
dataSe_50 = conSSe;

% concatenate data
dataAv = cat(1,dataAv_25,dataAv_30,dataAv_50);
dataGrAv = cat(1,dataGrAv_25,dataGrAv_30,dataGrAv_50);
dataSe = cat(1,dataSe_25,dataSe_30,dataSe_50);

% load model data
load('DD_ConSData.mat')

%% do some additional calculations

combNameS = split(combName,","); % split comb names at comma
nStPa = numel(stPaPlot);

nStims = 1; % 1: plots only low-freq. responses (remove line if both responses should be plotted)
nStPa = 1; % 1: plots only P2P responses (remove line if both responses should be plotted)

%% plot responses

if plotRespI==1
    for r = 1:nRepRate % for each rep rate
        pSmpl = size(dataGrAv{r},1); % number of points per standard
        timeTr = (0:nConS*pSmpl-1)/fs(r); % time trace for the whole standard cluster
        for s = 1:nStims % for each stimulus (A and B)
            switch s
                case 1
                    stimName = 'A';
                case 2
                    stimName = 'B';
            end
            % concatenate standard responses
            dataPlot = zeros(pSmpl*nConS,1);
            for cs = 1:nConS
                dataPlot(1+pSmpl*(cs-1):pSmpl*cs) = dataGrAv{r}(:,s,cs,filt);
            end
            figure('NumberTitle','off','Name',['conSResp_',stimName,'_RepRate',num2str(r),'_Filt',num2str(filt)],'Position',[0,0,figWid,figHei])
            title([convertStringsToChars(combNameS(s)),'-Freq. Chirp'])
            hold on
            plotConS = plot(timeTr,dataPlot,'k','Linewidth',2);
            for cn = 1:nConS
                x_min = timeTr(winAll{r,s}(:,1)+round(pSmpl*(cn-1)+1));
                x_max = timeTr(winAll{r,s}(:,2)+round(pSmpl*(cn-1)+1));
                x = [x_min;x_max;x_max;x_min];
                switch fixedY
                    case 0
                        y_limAuto = ylim;
                        y_min = y_limAuto(1);
                        y_max = y_limAuto(2);
                    case 1
                        y_min = y_lim(1);
                        y_max = y_lim(2);
                end
                y_temp = [y_min;y_min;y_max;y_max];
                y = repmat(y_temp,1,nWin);
                for w = 1:nWin
                    sigWin = patch('XData',x(:,w),'YData',y(:,w),'EdgeColor',winCol(w,:),'EdgeAlpha',1,'FaceColor','none','Linewidth',2,'LineStyle','--');
                    uistack(sigWin,'bottom')
                end
            end
            if fixedY==1
                ylim(y_lim)
            end
            set(gca,'FontSize',30,'Linewidth',2,'FontName','Arial')
            if saveI==1
                saveas(gcf,['conSResp_',stimName,'_RepRate',num2str(r),'_Filt',num2str(filt),'.jpg'])
                saveas(gcf,['conSResp_',stimName,'_RepRate',num2str(r),'_Filt',num2str(filt),'.svg'])
            end
        end
    end
end

%% plot response variables (with fitted curve)

if plotVarI==1
    for r = 1:nRepRate % for each rep rate
        pSmpl = size(dataGrAv{r},1); % number of points per standard
        x = (0:nConS*pSmpl-1)/fs(r); % time trace for the whole standard cluster
        for s = 1:nStims % for each stimulus (A and B)
            switch s
                case 1
                    stimName = 'A';
                case 2
                    stimName = 'B';
            end
            for pa = 1:nStPa % for each stat param
                switch pa
                    case 1 % P2P
                        yLabelNameL = 'ABR p-p Amp. [µV]';
                        yLabelNameR = 'Rel. ABR p-p Amp.';
                        win = 1; % analysis window
                    case 2 % latency
                        yLabelNameL = 'Peak V Latency [ms]';
                        yLabelNameR = 'Rel. Peak V Latency';
                        win = 2; % analysis window
                end
                stPaAv = reshape(mean(stPa(:,r,s,:,filt,win,stPaPlot(pa)),1),nConS,1);
                varMax = max(stPaAv);
                figure('NumberTitle','off','Name',['conSVar_',stimName,'_RepRate',num2str(r),'_Filt',num2str(filt),'_Param',num2str(pa)],'Position',[0,0,figWid,figHei])
                hold on
                box on
                title([convertStringsToChars(combNameS(s)),'-Freq. Chirp'])
                % plot data points and create absolute and relative y axis
                yyaxis left % absolute
                set(gca,'YColor','k')
                plotVar = plot(stPaAv,'k','LineStyle','none','Marker','.','MarkerSize',30,'Linewidth',2); % plot datapoints                
                % plot model
                switch modelTy
                    case 1 % linear regression model
                        plotMdl = plot(linFit{r,s,filt,win,stPaPlot(pa)}); % plot linear regression model
                        % model adjustments
                        plotMdl(1).Color = [0,0,0];
                        plotMdl(2).Color = [0,0,0];
                        plotMdl(3).Color = [0,0,0];
                        plotMdl(4).Color = [0,0,0];
                        % linFit{r,s,filt,win,stPaPlot(pa)} % print model variables
                    case 2 % flexible model
                        plotMdl = plot(flexFit{r,s,filt,win,stPaPlot(pa)},'k');  % plot flexible model
                        % model adjustments
                        plotMdl(1).Color = [0,0,0];
                        plotMdl.LineStyle = "-";
                        % flexFit{r,s,filt,win,stPaPlot(pa)}
                        % flexGof{r,s,filt,win,stPaPlot(pa)}
                end
                % adjust ylim
                yLimL = ylim;
                yMinL = yLimL(1)-diff(yLimL)*0.1;
                yMaxL = yLimL(2)+diff(yLimL)*0.1;
                ylim([yMinL,yMaxL]);
                yLimL = ylim;
                ylabel(yLabelNameL,'FontWeight','bold','FontSize',30)
                yyaxis right % relative
                set(gca,'YColor','k')
                ylim(yLimL/stPaAv(1)) % relative to first standard of cluster
                yLimR = ylim;
                ylabel(yLabelNameR,'FontWeight','bold','FontSize',30)
                % cosmetics
                set(plotMdl,'Linewidth',2);
                xticks(1:1:nConS)
                xticklabels(1:1:nConS)
                xlim([0,nConS+1])
                xLim = xlim;
                % create text
                switch modelTy
                    case 1
                        p = round(linFit{r,s,filt,win,stPaPlot(pa)}.Coefficients.pValue(2),3);
                        if p<0.001
                            ast = '***';
                        elseif p<0.01
                            ast = '**';
                        elseif p<0.05
                            ast = '*';
                        else
                            ast = '';
                        end
                        Rs = round(linFit{r,s,filt,win,stPaPlot(pa)}.Rsquared.Ordinary,3);
                        txt = ['\itp',ast,'\rm = ',num2str(p),', \itR^2\rm = ',num2str(Rs)];
                        text((xLim(2)*0.4),yLimR(2)-diff(yLimR)*0.1,txt,'FontSize',30)
                        legend([plotVar,plotMdl(2),plotMdl(3)],'Data','Lin. Reg. Model','95 % Conf. Int.','Location','northeast')
                    case 2

                end
                xlabel("Number of Consecutive Std. after Dev.",'FontWeight','bold','FontSize',30)
                title([convertStringsToChars(combNameS(s)),'-Freq. Chirp'])
                set(gca,'FontSize',30,'Linewidth',2,'FontName','Arial')
                if saveI==1
                    saveas(gcf,['conSVar_',stimName,'_RepRate',num2str(r),'_Filt',num2str(filt),'.jpg'])
                    saveas(gcf,['conSVar_',stimName,'_RepRate',num2str(r),'_Filt',num2str(filt),'.svg'])
                end
            end
        end
    end
end
bla = 1;