%analyze and plot the decoding data %NOTE: this program takes a really long time because of large ANOVAs clear all; load ('R_2R.mat'); load ('RAW.mat'); load ('PoolDec_2R.mat'); load ('UnitDec_2R.mat'); sucrose=[1 0.6 0.1]; maltodextrin=[.9 0.3 .9]; water=[0.00 0.75 0.75]; total=[0.3 0.1 0.8]; NAc=[0.5 0.1 0.8]; VP=[0.3 0.7 0.1]; NAcP{5,1}=[1 0.7 1]; NAcP{4,1}=[0.85 0.3 0.9]; NAcP{3,1}=[0.6 0.1 0.8]; NAcP{2,1}=[0.3 0.05 0.5]; NAcP{1,1}=[0.2 0 0.3]; VPP{5,1}=[0.6 1 0.25]; VPP{4,1}=[0.4 0.9 0.15]; VPP{3,1}=[0.3 0.7 0.1]; VPP{2,1}=[0.1 0.4 0.05]; VPP{1,1}=[0.03 0.2 0]; NAcShuff=[0.3 0.05 0.5]; VPShuff=[0.05 0.4 0]; %load parameters BinDura=R_2R.Param.BinDura; bins=R_2R.Param.bins; binint=R_2R.Param.binint; binstart=R_2R.Param.binstart; NumNeurons=R_2R.Param.NumNeurons; repetitions=length(PoolDec{1,1}.True{1,1}(:,1)); testbins=R_2R.Param.SelectiveBins; xaxis=linspace(binstart+BinDura(2)/2,binstart+(bins-1)*binint+BinDura(2)/2,bins); %include all bins %% single unit decoding %divide neurons up by region NAneurons=strcmp(R_2R.Ninfo(:,3),'NA'); VPneurons=strcmp(R_2R.Ninfo(:,3),'VP'); %if you change this to e=1:3, you can also look at nonselective neurons for e=1:2 %different selections of neurons figure; %pick which set of neurons: all, reward-specific, or non-reward specific if e==1 selection=R_2R.SucN<2; end %all neurons if e==2 selection=R_2R.SucN | R_2R.MalN; end %reward-selective neurons if e==3 selection=(R_2R.SucN | R_2R.MalN) == 0; end %reward-nonselective neurons %get average accuracy for each bin and run stats comparing region for i = 1:bins AvgAccNAc(1,i)=nanmean(UnitDec.True(NAneurons&selection,i)); %average accuracy NAc SEMAccNAc(1,i)=nanste(UnitDec.True(NAneurons&selection,i),1); %SEM AvgAccNAcShuff(1,i)=nanmean(UnitDec.Shuff(NAneurons&selection,i)); %average accuracy NAcShuff SEMAccNAcShuff(1,i)=nanste(UnitDec.Shuff(NAneurons&selection,i),1); %SEM AvgAccVP(1,i)=nanmean(UnitDec.True(VPneurons&selection,i)); %average accuracy NAcShuff SEMAccVP(1,i)=nanste(UnitDec.True(VPneurons&selection,i),1); %SEM AvgAccVPShuff(1,i)=nanmean(UnitDec.Shuff(VPneurons&selection,i)); %average accuracy NAcShuff SEMAccVPShuff(1,i)=nanste(UnitDec.Shuff(VPneurons&selection,i),1); %SEM end %prepare shading upSEMNAc=AvgAccNAc+SEMAccNAc; downSEMNAc=AvgAccNAc-SEMAccNAc; upSEMVP=AvgAccVP+SEMAccVP; downSEMVP=AvgAccVP-SEMAccVP; upSEMNAcShuff=AvgAccNAcShuff+SEMAccNAcShuff; downSEMNAcShuff=AvgAccNAcShuff-SEMAccNAcShuff; upSEMVPShuff=AvgAccVPShuff+SEMAccVPShuff; downSEMVPShuff=AvgAccVPShuff-SEMAccVPShuff; %plotting decoder accuracies over time subplot(2,3,1); hold on; plot(xaxis,AvgAccNAc(1:bins),'Color', NAc,'linewidth',1.5); %accumbens plot(xaxis,AvgAccVP(1:bins),'Color', VP,'linewidth',1.5); %vp %shuffled plot(xaxis,AvgAccNAcShuff(1:bins),'Color', 'k','linewidth',1.5); plot(xaxis,AvgAccVPShuff(1:bins),'Color', 'k','linewidth',1.5); %error patch([xaxis,xaxis(end:-1:1)],[upSEMNAc,downSEMNAc(end:-1:1)],NAc,'EdgeColor','none');alpha(0.5); patch([xaxis,xaxis(end:-1:1)],[upSEMNAcShuff,downSEMNAcShuff(end:-1:1)],'k','EdgeColor','none');alpha(0.5); patch([xaxis,xaxis(end:-1:1)],[upSEMVP,downSEMVP(end:-1:1)],VP,'EdgeColor','none');alpha(0.5); patch([xaxis,xaxis(end:-1:1)],[upSEMVPShuff,downSEMVPShuff(end:-1:1)],'k','EdgeColor','none');alpha(0.5); xlabel('Seconds from reward delivery'); ylabel('Accuracy'); title('Unit decoding accuracy'); legend('NAc units','VP units','Shuffled','Location','northwest'); axis([xaxis(1) xaxis(end) 0.47 max(AvgAccVP)+0.04]); %find and plot bins exceeding confidence interval for i=1:bins %confidence interval for NAc x = UnitDec.Shuff(NAneurons&selection,i); % Create Data SEM = nanstd(x)/sqrt(length(x)); % Standard Error ts = tinv([0.005 0.995],length(x)-1); % T-Score CI = nanmean(x) + ts*SEM; % Confidence Intervals if nanmean(UnitDec.True(NAneurons&selection,i))>CI(2) NAcConf(1,i)=1; else NAcConf(1,i)=0; end %confidence interval for VP x = UnitDec.Shuff(VPneurons&selection,i); % Create Data SEM = nanstd(x)/sqrt(length(x)); % Standard Error ts = tinv([0.005 0.995],length(x)-1); % T-Score CI = nanmean(x) + ts*SEM; % Confidence Intervals if nanmean(UnitDec.True(VPneurons&selection,i))>CI(2) VPConf(1,i)=1; else VPConf(1,i)=0; end end %find consecutive bins R_2R.UnitNAcComp{e,1}=zeros(1,bins); R_2R.UnitVPComp{e,1}=zeros(1,bins); for i=2:bins-1 if NAcConf(1,i)==1 if NAcConf(1,i-1)==1 | NAcConf(1,i+1)==1 R_2R.UnitNAcComp{e,1}(1,i)=1; end end if VPConf(1,i)==1 if VPConf(1,i-1)==1 | VPConf(1,i+1)==1 R_2R.UnitVPComp{e,1}(1,i)=1; end end end %check to see if the first bin above shuffled data would be different %if using fewer VP neurons to match NAc for j=1:20 for i=1:bins matchedneurons=cat(1,ones(sum(NAneurons&selection),1),zeros(sum(VPneurons&selection)-sum(NAneurons&selection),1)); matchedneurons=(matchedneurons(randperm(length(matchedneurons)))==1); %confidence interval for VP VPselectionSh=UnitDec.Shuff(VPneurons&selection,i); VPselection=UnitDec.True(VPneurons&selection,i); x = VPselectionSh(matchedneurons,1); % Create Data SEM = nanstd(x)/sqrt(length(x)); % Standard Error ts = tinv([0.005 0.995],length(x)-1); % T-Score CI = nanmean(x) + ts*SEM; % Confidence Intervals if nanmean(VPselection(matchedneurons,1))>CI(2) R_2R.VPaboveshuff{e,1}(j,i)=1; else R_2R.VPaboveshuff{e,1}(j,i)=0; end end end %plot plot(xaxis,R_2R.UnitVPComp{e,1}*0.48,'*','Color',VP); plot(xaxis,R_2R.UnitNAcComp{e,1}*0.485,'*','Color',NAc); %multiple comparisons for NAc vs VP %make a matrix indicating which region each neuron-bin comes from nbinregion=[]; binname=[]; rat=[]; for i=1:bins %total number of bins being tested nbinregion=cat(2,nbinregion,NAneurons); binname=cat(2,binname,i*ones(length(NAneurons),1)); rat=cat(2,rat,R_2R.Ninfo(:,4)); end testtrue=UnitDec.True(selection,testbins(1)+1:testbins(2)+1); %becaue the bin 3 is actually the 4th entry, etc testshuff=UnitDec.Shuff(selection,testbins(1)+1:testbins(2)+1); testregion=nbinregion(selection,testbins(1)+1:testbins(2)+1); testbinname=binname(selection,testbins(1)+1:testbins(2)+1); testrat=rat(selection,testbins(1)+1:testbins(2)+1); %find effects of real vs shuffled, region, and bins on accuracy [~,R_2R.UnitDecStatsSubj{e,1},R_2R.UnitDecStatsSubj{e,2}]=anovan(cat(1,testtrue(:),testshuff(:)),{cat(1,zeros(length(testtrue(:)),1),ones(length(testshuff(:)),1)),cat(1,testregion(:),testregion(:)),cat(1,testbinname(:),testbinname(:)),cat(1,testrat(:),testrat(:))},'varnames',{'real vs shuffled','region','bins','subject'},'random',[4],'nested',[0 0 0 0;0 0 0 0;0 0 0 0;0 1 0 0],'display','off','model','full'); % %do post-hoc comparisons % %this is commented out because can't do post-hoc comparison with nested anova (which we're using for including random effect of subject) % %if you do want to look at this in matlab, need to do the anova without subject (as written out here) % % [~,R_2R.UnitDecStats{e,1},R_2R.UnitDecStats{e,2}]=anovan(cat(1,testtrue(:),testshuff(:)),{cat(1,zeros(length(testtrue(:)),1),ones(length(testshuff(:)),1)),cat(1,testregion(:),testregion(:)),cat(1,testbinname(:),testbinname(:))},'varnames',{'real vs shuffled','region','bins'},'display','off','model','full'); % [c,~,~,names]=multcompare(R_2R.UnitDecStats{e,2},'Dimension',[1 2 3],'CType','tukey-kramer','display','off'); % % %find post-hoc differences % R_2R.UnitNAcVPComp{e,1}=NaN(2,bins); % for i=1:1+testbins(2)-testbins(1) % %NAc vs VP % Sel=c(:,1)==4*(i-1)+1 & c(:,2)==4*(i-1)+3; % if c(Sel,6)<0.05 R_2R.UnitNAcVPComp{e,1}(1,i+testbins(1))=1; else R_2R.UnitNAcVPComp{e,1}(1,i+testbins(1))=0; end % R_2R.UnitNAcVPComp{e,1}(2,i+testbins(1))=c(Sel,6); % end % % %add it to plot % plot(xaxis,R_2R.UnitNAcVPComp{e,1}(1,:)*(max(AvgAccVP)+0.015),'*','Color','k'); %plotting CDF at peak bin subplot(2,3,4) hold on; %NAc [~,NAcbin]=max(AvgAccNAc); [cdfNAc,xNAc] = ecdf(UnitDec.True(NAneurons&selection,NAcbin)); [cdfNAcSh,xNAcSh] = ecdf(UnitDec.Shuff(NAneurons&selection,NAcbin)); plot(xNAc,cdfNAc,'Color',NAc,'linewidth',1.5); %VP [~,VPbin]=max(AvgAccVP); [cdfVP,xVP] = ecdf(UnitDec.True(VPneurons&selection,VPbin)); [cdfVPSh,xVPSh] = ecdf(UnitDec.Shuff(VPneurons&selection,VPbin)); plot(xVP,cdfVP,'Color',VP,'linewidth',1.5); %shuffled plot(xNAcSh,cdfNAcSh,'Color','k','linewidth',1.5); xlabel('Decoding accuracy') plot(xVPSh,cdfVPSh,'Color','k','linewidth',1.5); axis([0 1 0 1]); plot([0.5 0.5],[0 1],':','color','k','linewidth',0.75); title(['Accuracy at peak bin: ' num2str(((NAcbin-1)*binint)+(binstart+BinDura(2)/2)) 's NAc, ' num2str(((VPbin-1)*binint)+(binstart+BinDura(2)/2)) 's VP']); xlabel('Decoding accuracy'); ylabel('Cumulative fraction of population'); legend('NAc units','VP units','Shuffled','Location','northwest'); %stats comparing peak bins [~,R_2R.UnitDecPeakBinSubj{e,1},~]=anovan(cat(1,UnitDec.True(NAneurons&selection,NAcbin),UnitDec.Shuff(NAneurons&selection,NAcbin),UnitDec.True(VPneurons&selection,VPbin),UnitDec.Shuff(VPneurons&selection,VPbin)),... {cat(1,zeros(sum(NAneurons&selection),1),ones(sum(NAneurons&selection),1),zeros(sum(VPneurons&selection),1),ones(sum(VPneurons&selection),1)),... cat(1,zeros(sum(NAneurons&selection),1),zeros(sum(NAneurons&selection),1),ones(sum(VPneurons&selection),1),ones(sum(VPneurons&selection),1)),... cat(1,R_2R.Ninfo(NAneurons&selection,4),R_2R.Ninfo(NAneurons&selection,4),R_2R.Ninfo(VPneurons&selection,4),R_2R.Ninfo(VPneurons&selection,4))},'varnames',{'real vs shuffled','region','subject'},'random',[3],'nested',[0 0 0;0 0 0;0 1 0],'display','off','model','full'); %% pooled decoding %reset matrices for stats analysis A=[]; B=[]; C=[]; D=[]; %get average accuracy for each bin for v=1:length(PoolDec{e,1}.True) %each condition (number of neurons used) for i = 1:bins %NAc if length(PoolDec{e,1}.True{v,1})>1 %in case analysis wasn't run because not enough neurons PoolAccNAc{v,1}(1,i)=nanmean(PoolDec{e,1}.True{v,1}(:,i)); %average accuracy NAc PoolSEMNAc{v,1}(1,i)=nanste(PoolDec{e,1}.True{v,1}(:,i),1); %SEM PoolAccNAcShuff{v,1}(1,i)=nanmean(PoolDec{e,1}.Shuff{v,1}(:,i)); %average accuracy NAcShuff PoolSEMNAcShuff{v,1}(1,i)=nanste(PoolDec{e,1}.Shuff{v,1}(:,i),1); %SEM else PoolAccNAc{v,1}(1,i)=NaN; PoolSEMNAc{v,1}(1,i)=NaN; PoolAccNAcShuff{v,1}(1,i)=NaN; PoolSEMNAcShuff{v,1}(1,i)=NaN; end %VP if length(PoolDec{e,2}.True{v,1})>1 %in case analysis wasn't run because not enough neurons PoolAccVP{v,1}(1,i)=nanmean(PoolDec{e,2}.True{v,1}(:,i)); %average accuracy VP PoolSEMVP{v,1}(1,i)=nanste(PoolDec{e,2}.True{v,1}(:,i),1); %SEM PoolAccVPShuff{v,1}(1,i)=nanmean(PoolDec{e,2}.Shuff{v,1}(:,i)); %average accuracy VPshuff PoolSEMVPShuff{v,1}(1,i)=nanste(PoolDec{e,2}.Shuff{v,1}(:,i),1); %SEM else PoolAccVP{v,1}(1,i)=NaN; PoolSEMVP{v,1}(1,i)=NaN; PoolAccVPShuff{v,1}(1,i)=NaN; PoolSEMVPShuff{v,1}(1,i)=NaN; end end %find the time of the most accurate bin for each of the repetitions %(Following reward delivery) time0bin=round((((0-BinDura(2)/2)-binstart)/binint)); for j = 1:repetitions %NAc if length(PoolDec{e,1}.True{v,1})>1 %in case analysis wasn't run because not enough neurons [~,PeakBinsNAc{v,1}(j,1)]=max(PoolDec{e,1}.True{v,1}(j,time0bin+1:end)); PeakBinsNAc{v,1}(j,2)=(PeakBinsNAc{v,1}(j,1)-1)*binint; else PeakBinsNAc{v,1}(j,1)=NaN; PeakBinsNAc{v,1}(j,2)=NaN; end %VP if length(PoolDec{e,2}.True{v,1})>1 %in case analysis wasn't run because not enough neurons [~,PeakBinsVP{v,1}(j,1)]=max(PoolDec{e,2}.True{v,1}(j,time0bin+1:end)); PeakBinsVP{v,1}(j,2)=(PeakBinsVP{v,1}(j,1)-1)*binint; else PeakBinsVP{v,1}(j,1)=NaN; PeakBinsVP{v,1}(j,2)=NaN; end end %prepare shading PupSEMNAc{v,1}=PoolAccNAc{v,1}+PoolSEMNAc{v,1}; PdownSEMNAc{v,1}=PoolAccNAc{v,1}-PoolSEMNAc{v,1}; PupSEMVP{v,1}=PoolAccVP{v,1}+PoolSEMVP{v,1}; PdownSEMVP{v,1}=PoolAccVP{v,1}-PoolSEMVP{v,1}; PupSEMNAcShuff{v,1}=PoolAccNAcShuff{v,1}+PoolSEMNAcShuff{v,1}; PdownSEMNAcShuff{v,1}=PoolAccNAcShuff{v,1}-PoolSEMNAcShuff{v,1}; PupSEMVPShuff{v,1}=PoolAccVPShuff{v,1}+PoolSEMVPShuff{v,1}; PdownSEMVPShuff{v,1}=PoolAccVPShuff{v,1}-PoolSEMVPShuff{v,1}; %here I plot the main lines for each condition that will have an entry %in the legend (legend goes by order plotted) %plotting decoder accuracy over time subplot(2,3,2); %accumbens hold on; plot(xaxis,PoolAccNAc{v,1}(1:bins),'Color', NAcP{v,1},'linewidth',1); %plot(xaxis,PoolAccNAcShuff{v,1}(1:66),'Color', NAcShuff,'linewidth',3); %plot(xaxis,NAcSig{v,1}-0.53,'*','Color', 'k'); %patch([xaxis,xaxis(end:-1:1)],[PupSEMNAc{v,1},PdownSEMNAc{v,1}(end:-1:1)],NAc{v,1},'EdgeColor','none');alpha(0.5); %patch([xaxis,xaxis(end:-1:1)],[PupSEMNAcShuff{v,1},PdownSEMNAcShuff{v,1}(end:-1:1)],NAcShuff,'EdgeColor','none');alpha(0.5); xlabel('Seconds post reward delivery'); ylabel('Accuracy'); title('NAc decoding accuracy'); axis([xaxis(1) xaxis(end) 0.4 1]); subplot(2,3,3); %VP hold on; plot(xaxis,PoolAccVP{v,1}(1:bins),'Color', VPP{v,1},'linewidth',1); %plot(xaxis,PoolAccVPShuff{v,1}(1:66),'Color', VPShuff,'linewidth',3); %plot(xaxis,VPSig{v,1}-0.53,'*','Color','k'); %patch([xaxis,xaxis(end:-1:1)],[PupSEMVP{v,1},PdownSEMVP{v,1}(end:-1:1)],VP{v,1},'EdgeColor','none');alpha(0.5); %patch([xaxis,xaxis(end:-1:1)],[PupSEMVPShuff{v,1},PdownSEMVPShuff{v,1}(end:-1:1)],VPShuff,'EdgeColor','none');alpha(0.5); xlabel('Seconds post reward delivery'); ylabel('Accuracy'); title('VP decoding accuracy'); axis([xaxis(1) xaxis(end) 0.4 1]); %comparison at best bin subplot(2,3,5); hold on; [~,NAcbin]=max(PoolAccNAc{v,1}); [~,VPbin]=max(PoolAccVP{v,1}); if length(PoolDec{e,1}.True{v,1})>1 && length(PoolDec{e,2}.True{v,1})>1 %make matrices that will include best bin at all levels A=cat(1,A,PoolDec{e,1}.True{v,1}(:,NAcbin)); B=cat(1,B,PoolDec{e,1}.Shuff{v,1}(:,NAcbin)); C=cat(1,C,PoolDec{e,2}.True{v,1}(:,VPbin)); D=cat(1,D,PoolDec{e,2}.Shuff{v,1}(:,VPbin)); end errorbar(NumNeurons(v),PoolAccNAc{v,1}(NAcbin),PoolSEMNAc{v,1}(NAcbin),'*','Color', NAcP{v,1},'linewidth',1.5); errorbar(NumNeurons(v),PoolAccVP{v,1}(VPbin),PoolSEMVP{v,1}(VPbin),'*','Color', VPP{v,1},'linewidth',1.5); xlabel('Neurons used'); ylabel('Accuracy'); title('Accuracy for most predictive bin'); legend('NAc','VP','Location','southeast'); axis([0 155 0.5 1]); %plotting time of most accurate bin subplot(2,3,6); hold on; errorbar(NumNeurons(v),nanmean(PeakBinsNAc{v,1}(:,2)),nanste(PeakBinsNAc{v,1}(:,2),1),'*','Color', NAcP{v,1},'linewidth',1.5); errorbar(NumNeurons(v),nanmean(PeakBinsVP{v,1}(:,2)),nanste(PeakBinsVP{v,1}(:,2),1),'*','Color', VPP{v,1},'linewidth',1.5); xlabel('Neurons used'); ylabel('Seconds post reward delivery') title('Time of most accurate bin'); %legend('NAc','VP','Location','southeast'); axis([0 155 0 2.5]); legend('NAc','VP','location','southeast') end %conditions %now plot the shading, significance, and the shuffled separately for v=1:length(PoolDec{e,1}.True) %each condition (number of neurons used) %plotting decoder accuracy over time subplot(2,3,2); %accumbens hold on; plot(xaxis,PoolAccNAcShuff{v,1}(1:66),'Color', 'k','linewidth',1); patch([xaxis,xaxis(end:-1:1)],[PupSEMNAc{v,1},PdownSEMNAc{v,1}(end:-1:1)],NAcP{v,1},'EdgeColor','none');alpha(0.5); patch([xaxis,xaxis(end:-1:1)],[PupSEMNAcShuff{v,1},PdownSEMNAcShuff{v,1}(end:-1:1)],'k','EdgeColor','none');alpha(0.5); xlabel('Seconds post reward delivery'); title('NAc decoding accuracy'); axis([xaxis(1) xaxis(end) 0.4 1]); legend('10 units','25 units','50 units','100 units','150 units','Shuffled','location','northwest') subplot(2,3,3); %VP hold on; plot(xaxis,PoolAccVPShuff{v,1}(1:66),'Color', 'k','linewidth',1); patch([xaxis,xaxis(end:-1:1)],[PupSEMVP{v,1},PdownSEMVP{v,1}(end:-1:1)],VPP{v,1},'EdgeColor','none');alpha(0.5); patch([xaxis,xaxis(end:-1:1)],[PupSEMVPShuff{v,1},PdownSEMVPShuff{v,1}(end:-1:1)],'k','EdgeColor','none');alpha(0.5); xlabel('Seconds post reward delivery'); title('VP decoding accuracy'); axis([xaxis(1) xaxis(end) 0.4 1]); legend('10 units','25 units','50 units','100 units','150 units','Shuffled','location','northwest') end %stats! %ANOVA for effects of number of neurons and region on time of peak accuracy %time PeakTimes=[]; for i=1:length(NumNeurons) if length(PoolDec{e,1}.True{i,1})>1 && length(PoolDec{e,2}.True{i,1})>1 PeakTimes(1:repetitions,i)=PeakBinsNAc{i,1}(:,2); PeakTimes(1+repetitions:repetitions+repetitions,i)=PeakBinsVP{i,1}(:,2); end end [~,R_2R.PoolDecTimeStats{e,1},R_2R.PoolDecTimeStats{e,2}]=anova2(PeakTimes,repetitions,'off'); %columns is ensemble size, rows is region %ANOVA for effect of ensemble size, region, and shuff vs true on accuracy %of peak bin size=[]; for i=1:length(NumNeurons) if length(PoolDec{e,1}.True{i,1})>1 && length(PoolDec{e,2}.True{i,1})>1 size=cat(1,size,i*ones(repetitions,1)); end end %with shuff vs true % size2=cat(1,size,size,size,size); % shuffd=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1)); % region=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1)); % [~,R_2R.AccAnova,~]=anovan(cat(1,A,B,C,D),{shuffd,region,size2},'varnames',{'real vs shuffled','region','ens size'},'display','off','model','full'); %without shuff vs true size3=cat(1,size,size); region3=cat(1,ones(length(A),1),zeros(length(C),1)); [~,R_2R.PoolDecAccStats{e,1},R_2R.PoolDecAccStats{e,2}]=anovan(cat(1,A,C),{region3,size3},'varnames',{'region','ens size'},'display','off','model','full'); end %selections save('R_2R.mat','R_2R');