NAcVP2018/MatlabScripts/g_Fig3FigS6Decoding.m

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));

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');

for e=1:3 %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=[];
    for i=1:bins
        nbinregion=cat(2,nbinregion,NAneurons);
        binname=cat(2,binname,i*ones(length(NAneurons),1));
    end
    
    testtrue=UnitDec.True(selection,:);
    testshuff=UnitDec.Shuff(selection,:);
    testregion=nbinregion(selection,:);
    testbins=binname(selection,:);
    
    
    %find effects of real vs shuffled, region, and bins on accuracy
    [~,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,testbins(:),testbins(:))},'varnames',{'real vs shuffled','region','bins'},'display','off','model','full');

    %do post-hoc comparisons
    [c,~,~,names]=multcompare(R_2R.UnitDecStats{e,2},'Dimension',[1 2 3],'CType','tukey-kramer','display','off');

    %find post-hoc differences 
    for i=1:bins
        %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)=1; else R_2R.UnitNAcVPComp{e,1}(1,i)=0; end
        R_2R.UnitNAcVPComp{e,1}(2,i)=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.UnitDecPeakBin{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))},'varnames',{'real vs shuffled','region'},'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
        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,:));
                PeakBinsNAc{v,1}(j,2)=binstart-binint+PeakBinsNAc{v,1}(j,1)*binint+BinDura(2)/2;
            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,:));
                PeakBinsVP{v,1}(j,2)=binstart-binint+PeakBinsVP{v,1}(j,1)*binint+BinDura(2)/2;
            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:66),'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:66),'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');