%% stats for linear models comparison %% DATA =================================================================== clear stats ecc=[]; sz=[]; % eccentricity and size signal = []; % categorical 1=MRI, 2=MUA, 3-7=LFP bands area=[]; % categorical 1=V1,4=V4 model=[]; % categorical 1=lin, 2=lin_ng, 3=css, 4=dog % v1 ecc = [ecc; ... mri1(m).ECC; ... mua1(m).ECC; ... ]; sz = [ sz; ... mri1(m).S; ... mua1(m).S; ... ]; signal = [signal;... 1*ones(size(mri1(m).S)); ... 2*ones(size(mua1(m).S)); ... ]; area = [area; ... 1*ones(size(mri1(m).S)); ... 1*ones(size(mua1(m).S)); ... ]; for i=1:5 ecc = [ecc; lfp1(i,m).ECC]; sz = [ sz; lfp1(i,m).S]; signal = [signal; (2+i)*ones(size(lfp1(i,m).S))]; area = [area; 1*ones(size(lfp1(i,m).S))]; end % v4 ecc = [ecc; ... mri4(m).ECC; ... mua4(m).ECC; ... ]; sz = [ sz; ... mri4(m).S; ... mua4(m).S; ... ]; signal = [signal;... 1*ones(size(mri4(m).S)); ... 2*ones(size(mua4(m).S)); ... ]; area = [area; ... 4*ones(size(mri4(m).S)); ... 4*ones(size(mua4(m).S)); ... ]; for i=1:5 ecc = [ecc; lfp4(i,m).ECC]; sz = [ sz; lfp4(i,m).S]; signal = [signal; (2+i)*ones(size(lfp4(i,m).S))]; area = [area; 4*ones(size(lfp4(i,m).S))]; end % sigtype = {1,2,3,4,5,6,7;'mri' ,'mua','alpha','beta','theta','hgam','lgam'}; signal_num=signal; for i=1:length(signal_num) sn{i,1}=sigtype{2,signal_num(i)}; end signal=categorical(signal); area=categorical(area); areaname = {' V1', 'V4'}; DT=table(signal,area, ecc, sz); %% PLOT SCATTER AND FITS ================================================== f_xmod = figure; set(f_xmod,'Position',[100 100 2000 800]); CLR=[... 0 0 0;... 0.85,0.33,0.10;... 0.93,0.69,0.13;... 0.47,0.67,0.19;... 0.49,0.18,0.56;... 0.30,0.75,0.93;... 0.00,0.45,0.74;... ]; % do the modelfits and plot scatter and fitlines % v1 ----- DT1=DT(DT.area==categorical(1) & DT.ecc<=MaxECC(1),:); mdl_v1 = fitlm(DT1,'sz ~ 1 + signal*ecc'); slidx = 1+ length(mdl_v1.Coefficients.Estimate)/2; v1idx=[]; for j=1:length(mdl_v1.CoefficientNames) if mdl_v1.CoefficientNames{j}(1:3) == 'sig' & ... mdl_v1.CoefficientNames{j}(end-2:end) ~= 'ecc' v1idx = [v1idx str2double(mdl_v1.CoefficientNames{j}(8:end))]; end end msz = 40; subplot(2,2,1); ll={}; US=[sigtype{1,:}]; hold on; ii=1; for i=1:length(US) sel=DT1.signal==categorical(i); if sum(sel)>0 scatter(DT1.ecc(sel),DT1.sz(sel),msz,... 'Marker','o',... 'MarkerEdgeColor',CLR(i,:),'MarkerFaceColor',CLR(i,:),... 'MarkerEdgeAlpha',0,'MarkerFaceAlpha',0.2); ll{ii}=sigtype{2,i}; ii=ii+1; end end legend(ll,'Location','NorthEastOutside'); title('V1');xlabel('Ecc');ylabel('Sz'); set(gca,'xlim',[0 12.5],'ylim',[0 7],'TickDir','out'); subplot(2,2,2); ll={}; US=[sigtype{1,:}]; hold on; ii=1; for i=1:length(US) sel=DT1.signal==categorical(i); if sum(sel)>0 x=[0 1.5*MaxECC(1)]; if i==1 y=[... mdl_v1.Coefficients.Estimate(1) + ... mdl_v1.Coefficients.Estimate(slidx).*x ... ]; plot(x,y, 'LineWidth',5, 'Color',CLR(i,:)); else j = find(v1idx==i,1,'first'); y=[... mdl_v1.Coefficients.Estimate(1) + ... mdl_v1.Coefficients.Estimate(1+j) + ... (mdl_v1.Coefficients.Estimate(slidx) + ... mdl_v1.Coefficients.Estimate(slidx+j)).*x ... ]; plot(x,y, 'LineWidth',3, 'Color',CLR(i,:)); end ll{ii}=sigtype{2,i}; ii=ii+1; end end legend(ll,'Location','NorthEastOutside'); title('V1 FIT');xlabel('Ecc');ylabel('Sz'); set(gca,'xlim',[0 12.5],'ylim',[0 7],'TickDir','out'); % v4 DT4=DT(DT.area==categorical(4) & DT.ecc<=MaxECC(2),:); mdl_v4 = fitlm(DT4,'sz ~ 1 + signal*ecc'); v4idx=[]; for j=1:length(mdl_v4.CoefficientNames) if mdl_v4.CoefficientNames{j}(1:3) == 'sig' & ... mdl_v4.CoefficientNames{j}(end-2:end) ~= 'ecc' v4idx = [v4idx str2double(mdl_v4.CoefficientNames{j}(8:end))]; end end subplot(2,2,3); ll={}; US=[sigtype{1,:}]; hold on; ii=1; for i=1:length(US) sel=DT4.signal==categorical(i); if sum(sel)>0 %scatter(DT4.ecc(sel),DT4.sz(sel),'filled','CData',CLR(i,:)); scatter(DT4.ecc(sel),DT4.sz(sel),msz,... 'Marker','o',... 'MarkerEdgeColor',CLR(i,:),'MarkerFaceColor',CLR(i,:),... 'MarkerEdgeAlpha',0,'MarkerFaceAlpha',0.25); ll{ii}=sigtype{2,i}; ii=ii+1; end end legend(ll,'Location','NorthEastOutside'); title('V4');xlabel('Ecc');ylabel('Sz'); set(gca,'xlim',[0 12.5],'ylim',[0 7],'TickDir','out'); subplot(2,2,4); ll={}; US=[sigtype{1,:}]; hold on; ii=1; for i=1:length(US) sel=DT4.signal==categorical(i); if sum(sel)>0 x=[0 1.5*MaxECC(2)]; if i==1 y=[... mdl_v4.Coefficients.Estimate(1) + ... mdl_v4.Coefficients.Estimate(slidx).*x ... ]; plot(x,y, 'LineWidth',5, 'Color',CLR(i,:)); else j = find(v4idx==i,1,'first'); y=[... mdl_v4.Coefficients.Estimate(1) + ... mdl_v4.Coefficients.Estimate(1+j) + ... (mdl_v4.Coefficients.Estimate(slidx) + ... mdl_v4.Coefficients.Estimate(slidx+j)).*x ... ]; plot(x,y, 'LineWidth',3, 'Color',CLR(i,:)); end ll{ii}=sigtype{2,i}; ii=ii+1; end end legend(ll,'Location','NorthEastOutside'); title('V4 FIT');xlabel('Ecc');ylabel('Sz'); set(gca,'xlim',[0 12.5],'ylim',[0 7],'TickDir','out'); sgtitle(['MODEL: ' MMS{m} ', Rth_mri: ' num2str(Rth_mri) ... ', Rth_ephys:' num2str(Rth_ephys)],'interpreter','none') if SaveFigs saveas(f_xmod,fullfile(pngfld, ['XMOD_REGR_' MMS{m} '.png'])); saveas(f_xmod,fullfile(svgfld, ['XMOD_REGR_' MMS{m} '.svg'])); end if CloseFigs; close(f_xmod); end %% WHICH SIGNAL TYPES HAVE SIGNIFICANT ECC-SZ RELATION ==================== areas = [1 4]; st=[sigtype{1,:}];sn=sigtype(2,:); for a = 1:length(areas) % loop over areas stats.area(a).signalswithslope = []; for s = st % loop over signals sel = ... DT.area==categorical(areas(a)) & ... DT.signal==categorical(s) & ... DT.ecc<=MaxECC(a); if sum(sel)>1 tbl = DT(sel,:); mdl = fitlm(tbl,'sz ~ 1 + ecc'); stats.area(a).signal(s).name = sn{s}; stats.area(a).signal(s).mdl = mdl; stats.area(a).signal(s).anova = anova(mdl); stats.area(a).signal(s).ic = mdl.Coefficients(1,:); stats.area(a).signal(s).sl = mdl.Coefficients(2,:); CI = coefCI(mdl); stats.area(a).signal(s).icCI = CI(1,:); stats.area(a).signal(s).slCI = CI(2,:); if stats.area(a).signal(s).sl.pValue < 0.05 stats.area(a).signalswithslope = [... stats.area(a).signalswithslope s]; end fprintf(['Area ' num2str(areas(a)) ', Signal: ' num2str(sn{s}) ... ', slope: ' num2str(stats.area(a).signal(s).sl.Estimate) ... ', n = ' num2str(stats.area(a).signal(s).anova.DF(2)+1) ... ', t = ' num2str(stats.area(a).signal(s).sl.tStat) ... ', p = ' num2str(stats.area(a).signal(s).sl.pValue) '\n']); end end end %% TEST SIGNALS WITH SIGNIFICANT SLOPE TOGETHER IN ONE LINEAR MODEL ======= for a = 1:length(areas) % loop over areas sel = ... DT.area==categorical(areas(a)) & ... ismember(DT.signal,categorical(stats.area(a).signalswithslope)) & ... DT.ecc<=MaxECC(a); % recreate a table to avoid empty signal categories signal2 = zeros(size(signal_num)); for i=1:length(stats.area(a).signalswithslope) signal2(signal_num==stats.area(a).signalswithslope(i))=i; end signal2 = signal2(sel); area2 = area(sel); ecc2=ecc(sel); sz2=sz(sel); signal2=categorical(signal2); tbl2=table(signal2,area2, ecc2, sz2); mdl2 = fitlm(tbl2,'sz2 ~ signal2*ecc2'); stats.area(a).full.mdl = mdl2; stats.area(a).full.anova = anova(mdl2); stats.area(a).full.signals = sn(stats.area(a).signalswithslope); fprintf('========================================\n') fprintf(['All significant signal types in one model V' num2str(areas(a)) '\n']); fprintf('========================================\n') fprintf('Eccentricity\n') fprintf(['F = ' num2str(stats.area(a).full.anova.F(2)) ... ', df = ' num2str(stats.area(a).full.anova.DF(2)) ... ', p = ' num2str(stats.area(a).full.anova.pValue(2)) '\n']); fprintf('Signal\n') fprintf(['F = ' num2str(stats.area(a).full.anova.F(1)) ... ', df = ' num2str(stats.area(a).full.anova.DF(1)) ... ', p = ' num2str(stats.area(a).full.anova.pValue(1)) '\n']); fprintf('Ecc*Signal\n') fprintf(['F = ' num2str(stats.area(a).full.anova.F(3)) ... ', df = ' num2str(stats.area(a).full.anova.DF(3)) ... ', p = ' num2str(stats.area(a).full.anova.pValue(3)) '\n']); end clear signal2; %% TEST SIGNALS WITH SIGNIFICANT SLOPE AGAINST MRI ======================== for a = 1:length(areas) % loop over areas ss = stats.area(a).signalswithslope; ss(ss==1)=[]; sidx=1; fprintf('========================================\n') fprintf(['MRI vs EPHYS V' num2str(areas(a)) '\n']); fprintf('========================================\n') for s = ss sel = ... DT.area == categorical(areas(a)) & ... ismember(DT.signal,categorical([1 s])) & ... DT.ecc <= MaxECC(a) ; % recreate a table to avoid empty signal categories signal2 = zeros(size(signal_num)); for i=1:length(stats.area(a).signalswithslope) signal2(signal_num==stats.area(a).signalswithslope(i))=i; end signal2 = signal2(sel); area2 = area(sel); ecc2=ecc(sel); sz2=sz(sel); signal2=categorical(signal2); tbl2=table(signal2,area2, ecc2, sz2); mdl2 = fitlm(tbl2,'sz2 ~ signal2*ecc2'); stats.area(a).sig_vsMRI(sidx).name = sn{s}; stats.area(a).sig_vsMRI(sidx).mdl = mdl2; stats.area(a).sig_vsMRI(sidx).anova = anova(mdl2); fprintf(['MRI vs ' stats.area(a).sig_vsMRI(sidx).name '\n']); fprintf('Eccentricity\n') fprintf(['F = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.F(2)) ... ', df = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.DF(2)) ... ', p = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.pValue(2)) '\n']); fprintf('Signal\n') fprintf(['F = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.F(1)) ... ', df = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.DF(1)) ... ', p = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.pValue(1)) '\n']); fprintf('Ecc*Signal\n') fprintf(['F = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.F(3)) ... ', df = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.DF(3)) ... ', p = ' num2str(stats.area(a).sig_vsMRI(sidx).anova.pValue(3)) '\n']); sidx=sidx+1; end end %% PLOT SLOPES WITH CI AND INDICATE SIGN DIFF FROM MRI ==================== f_slopes = figure; set(f_slopes,'Position',[100 100 1200 400]); % V1 ----- subplot(1,2,1); hold on; mri_sl = stats.area(1).signal(1).sl.Estimate; mri_sl_bot = stats.area(1).signal(1).slCI(1); mri_sl_top = stats.area(1).signal(1).slCI(2); x=[0 10]; xarea=[x fliplr(x)]; yarea=[mri_sl_bot mri_sl_bot mri_sl_top mri_sl_top]; fill(xarea,yarea,'r','EdgeColor','none','FaceAlpha',0.25) plot(x,[mri_sl mri_sl],'r','Linewidth',3); for s=2:7 NoData=false; if isempty(stats.area(1).signal(s).sl) NoData=true; end if~NoData Y=[stats.area(1).signal(s).sl.Estimate ... stats.area(1).signal(s).slCI(1) ... stats.area(1).signal(s).slCI(2)]; % check if this is significantly different from MRI isSim=false; for ss = 1:length(stats.area(1).sig_vsMRI) if strcmp(stats.area(1).signal(s).name,... stats.area(1).sig_vsMRI(ss).name) & ... stats.area(1).sig_vsMRI(ss).anova.pValue(3) > 0.05 isSim=true; end end if isSim bar(s,Y(1),'FaceColor',[0.4660, 0.6740, 0.1880]); errorbar(s,Y(1),Y(1)-Y(2),Y(3)-Y(1),'k','Linestyle','none'); elseif NoData % nothing to plot else bar(s,Y(1),'FaceColor',[.25 .25 .25]); errorbar(s,Y(1),Y(1)-Y(2),Y(3)-Y(1),'k','Linestyle','none'); end end end set(gca,'xlim',[1.5 7.5],'ylim',[-0.1 0.5]) set(gca,'xtick',1:7,'xticklabels',sigtype(2,:)) ylabel('Ecc-Sz SLOPE'); title('V1','interpreter','none'); % V4 ----- subplot(1,2,2); hold on; mri_sl = stats.area(2).signal(1).sl.Estimate; mri_sl_bot = stats.area(2).signal(1).slCI(1); mri_sl_top = stats.area(2).signal(1).slCI(2); x=[0 10]; xarea=[x fliplr(x)]; yarea=[mri_sl_bot mri_sl_bot mri_sl_top mri_sl_top]; fill(xarea,yarea,'r','EdgeColor','none','FaceAlpha',0.25) plot(x,[mri_sl mri_sl],'r','Linewidth',3); for s=2:7 NoData=false; if isempty(stats.area(2).signal(s).sl) NoData=true; end if~NoData Y=[stats.area(2).signal(s).sl.Estimate ... stats.area(2).signal(s).slCI(1) ... stats.area(2).signal(s).slCI(2)]; % check if this is significantly different from MRI isSim=false; for ss = 1:length(stats.area(2).sig_vsMRI) if strcmp(stats.area(2).signal(s).name,... stats.area(2).sig_vsMRI(ss).name) & ... stats.area(2).sig_vsMRI(ss).anova.pValue(3) > 0.05 isSim=true; end end if isSim bar(s,Y(1),'FaceColor',[0.4660, 0.6740, 0.1880]); errorbar(s,Y(1),Y(1)-Y(2),Y(3)-Y(1),'k','Linestyle','none'); elseif NoData % nothing to plot else bar(s,Y(1),'FaceColor',[.25 .25 .25]); errorbar(s,Y(1),Y(1)-Y(2),Y(3)-Y(1),'k','Linestyle','none'); end end end set(gca,'xlim',[1.5 7.5],'ylim',[-0.10 0.5]) set(gca,'xtick',1:7,'xticklabels',sigtype(2,:)) ylabel('Ecc-Sz SLOPE'); title('V4','interpreter','none'); sgtitle(['MODEL: ' MMS{m} ', Rth_mri: ' num2str(Rth_mri) ... ', Rth_ephys:' num2str(Rth_ephys)],'interpreter','none') if SaveFigs saveas(f_slopes,fullfile(pngfld, ['XMOD_SlopeComparison_' MMS{m} '.png'])); saveas(f_slopes,fullfile(svgfld, ['XMOD_SlopeComparison_' MMS{m} '.svg'])); end if CloseFigs; close(f_slopes); end