VentralPallidumAlcohol2019/AnalysisAndFigures/i_Figure4_DecodingSucroseVersusAlcohol.m

%Plotting decoding analysis
clear all;
load('UnitMdlDSSucrose.mat')
load('UnitMdlDSShuffSucrose.mat')
load('UnitMdlPavSucrose.mat')
load('UnitMdlPavShuffSucrose.mat')
if ~exist('RPPSall'),load('RPPSall.mat'); end
if ~exist('RPSall'),load('RPSall.mat'); end
load('UnitMdlDSAlcohol.mat')
load('UnitMdlDSShuffAlcohol.mat')
load('UnitMdlPavAlcohol.mat')
load('UnitMdlPavShuffAlcohol.mat')
if ~exist('RvppaALL'),load('RvppaALL.mat'); end
if ~exist('RvpppaALL'),load('RvpppaALL.mat'); end

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];
DS=[0    0.4470    0.7410];
Pav=[0.8500    0.3250    0.0980];
DSAlcohol=[0.4940    0.1840    0.5560];
PavAlcohol=[0.9290    0.6940    0.1250];
DSShuff=[0 0 0];
PavShuff=[0 0 0];

xaxis=linspace(-0.75,3.15,13);

selection1=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;
selection2=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;
selection3=RvpppaALL.Structure==10 & RvpppaALL.CSPlusRatio>=.5 & (RvpppaALL.CSPlusRatio./(RvpppaALL.CSPlusRatio+RvpppaALL.CSMinusRatio))>=.7;
selection4=RvppaALL.Structure==10 & RvppaALL.CSPlusRatio>=.5 & (RvppaALL.CSPlusRatio./(RvppaALL.CSPlusRatio+RvppaALL.CSMinusRatio))>=.7;

%% Figure code

figure;

%get average accuracy for each bin
for i = 1:13
    AvgAccDSSucrose(1,i)=nanmean(UnitMdlAccDSSucrose(selection1,i)); %average accuracy DS
    SEMAccDSSucrose(1,i)=nanste(UnitMdlAccDSSucrose(selection1,i),1); %SEM
    AvgAccDSShuffSucrose(1,i)=nanmean(UnitMdlAccDSShuffSucrose(selection1,i)); %average accuracy DSShuff
    SEMAccDSShuffSucrose(1,i)=nanste(UnitMdlAccDSShuffSucrose(selection1,i),1); %SEM
    AvgAccPavSucrose(1,i)=nanmean(UnitMdlAccPavSucrose(selection2,i)); %average accuracy DSShuff
    SEMAccPavSucrose(1,i)=nanste(UnitMdlAccPavSucrose(selection2,i),1); %SEM
    AvgAccPavShuffSucrose(1,i)=nanmean(UnitMdlAccPavShuffSucrose(selection2,i)); %average accuracy DSShuff
    SEMAccPavShuffSucrose(1,i)=nanste(UnitMdlAccPavShuffSucrose(selection2,i),1); %SEM
    
    DSSigSucrose(1,i) = ttest2(UnitMdlAccDSSucrose(:,i),UnitMdlAccDSShuffSucrose(selection1,i),'Vartype','unequal'); %welch's t-test
    PavSigSucrose(1,i) = ttest2(UnitMdlAccPavSucrose(:,i),UnitMdlAccPavShuffSucrose(selection2,i),'Vartype','unequal'); %welch's t-test
end

%get normalized accuracy increase over shuffled for each neuron, then take
%average and find SEM
for i = 1:13
    for j = 1:length(UnitMdlAccDSSucrose)
        NormAccDSSucrose(j,i) = UnitMdlAccDSSucrose(j,i)-nanmean(nanmean(UnitMdlAccDSShuffSucrose(:,i)));
    end
    for k = 1:length(UnitMdlAccPavSucrose)
        NormAccPavSucrose(k,i) = UnitMdlAccPavSucrose(k,i)-nanmean(nanmean(UnitMdlAccPavShuffSucrose(:,i)));
    end
    
    AvgNormDSSucrose(1,i)=nanmean(NormAccDSSucrose(selection1,i)); %average norm accuracy DS
    SEMNormDSSucrose(1,i)=nanste(NormAccDSSucrose(selection1,i),1); %SEM
    AvgNormPavSucrose(1,i)=nanmean(NormAccPavSucrose(selection2,i)); %average norm accuracy DS
    SEMNormPavSucrose(1,i)=nanste(NormAccPavSucrose(selection2,i),1); %SEM
    
    %NormSig(1,i) = ttest2(NormAccDS(:,i),NormAccPav(:,i),'Vartype','unequal'); %welch's t-test
    [~,NormSig(1,i)] = ranksum(NormAccDSSucrose(selection1,i),NormAccPavSucrose(selection2,i)); %wilcoxon's ranksum test
end

%prepare shading
upSEMDS=AvgAccDSSucrose+SEMAccDSSucrose;
downSEMDS=AvgAccDSSucrose-SEMAccDSSucrose;
upSEMPav=AvgAccPavSucrose+SEMAccPavSucrose;
downSEMPav=AvgAccPavSucrose-SEMAccPavSucrose;
upSEMDSShuff=AvgAccDSShuffSucrose+SEMAccDSShuffSucrose;
downSEMDSShuff=AvgAccDSShuffSucrose-SEMAccDSShuffSucrose;
upSEMPavShuff=AvgAccPavShuffSucrose+SEMAccPavShuffSucrose;
downSEMPavShuff=AvgAccPavShuffSucrose-SEMAccPavShuffSucrose;
upSEMNormDS=AvgNormDSSucrose+SEMNormDSSucrose;
downSEMNormDS=AvgNormDSSucrose-SEMNormDSSucrose;
upSEMNormPav=AvgNormPavSucrose+SEMNormPavSucrose;
downSEMNormPav=AvgNormPavSucrose-SEMNormPavSucrose;

%plotting decoder accuracy over time
subplot(3,4,1); %accumbens
hold on;
plot(xaxis,AvgAccDSSucrose(1:13),'Color', DS,'linewidth',3);
plot(xaxis,AvgAccDSShuffSucrose(1:13),'Color', DSShuff,'linewidth',3);
%plot(xaxis,DSSigSucrose-0.51,'*','Color', 'k');
patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DS,'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
xlabel('Seconds post cue');
legend('Sucrose DS','Shuffled','Location','best');
%title('Sucrose DS single unit decoding accuracy');
ylabel('Single Unit Decoding Accuracy');
axis([-.75 3.15 0.45 0.75]);

subplot(3,4,2); %Pav
hold on;
plot(xaxis,AvgAccPavSucrose(1:13),'Color', Pav,'linewidth',3);
plot(xaxis,AvgAccPavShuffSucrose(1:13),'Color', PavShuff,'linewidth',3);
%plot(xaxis,PavSig-0.51,'*','Color', 'k');
patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],Pav,'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
xlabel('Seconds post cue');
%title('Sucrose CS+ single unit decoding accuracy');
legend('Sucrose CS+','Shuffled','Location','northeast');
axis([-.75 3.15 0.45 0.75]);

%get average accuracy for each bin
for i = 1:13
    AvgAccDSAlcohol(1,i)=nanmean(UnitMdlAccDSAlcohol(selection3,i)); %average accuracy DS
    SEMAccDSAlcohol(1,i)=nanste(UnitMdlAccDSAlcohol(selection3,i),1); %SEM
    AvgAccDSShuffAlcohol(1,i)=nanmean(UnitMdlAccDSShuffAlcohol(selection3,i)); %average accuracy DSShuff
    SEMAccDSShuffAlcohol(1,i)=nanste(UnitMdlAccDSShuffAlcohol(selection3,i),1); %SEM
    AvgAccPavAlcohol(1,i)=nanmean(UnitMdlAccPavAlcohol(selection4,i)); %average accuracy DSShuff
    SEMAccPavAlcohol(1,i)=nanste(UnitMdlAccPavAlcohol(selection4,i),1); %SEM
    AvgAccPavShuffAlcohol(1,i)=nanmean(UnitMdlAccPavShuffAlcohol(selection4,i)); %average accuracy DSShuff
    SEMAccPavShuffAlcohol(1,i)=nanste(UnitMdlAccPavShuffAlcohol(selection4,i),1); %SEM
    
    DSSigAlcohol(1,i) = ttest2(UnitMdlAccDSAlcohol(:,i),UnitMdlAccDSShuffAlcohol(selection3,i),'Vartype','unequal'); %welch's t-test
    PavSigAlcohol(1,i) = ttest2(UnitMdlAccPavAlcohol(:,i),UnitMdlAccPavShuffAlcohol(selection4,i),'Vartype','unequal'); %welch's t-test
end

%get normalized accuracy increase over shuffled for each neuron, then take
%average and find SEM
for i = 1:13
    for j = 1:length(UnitMdlAccDSAlcohol)
        NormAccDSAlcohol(j,i) = UnitMdlAccDSAlcohol(j,i)-nanmean(nanmean(UnitMdlAccDSShuffAlcohol(:,i)));
    end
    for k = 1:length(UnitMdlAccPavAlcohol)
        NormAccPavAlcohol(k,i) = UnitMdlAccPavAlcohol(k,i)-nanmean(nanmean(UnitMdlAccPavShuffAlcohol(:,i)));
    end
    
    AvgNormDSAlcohol(1,i)=nanmean(NormAccDSAlcohol(selection3,i)); %average norm accuracy DS
    SEMNormDSAlcohol(1,i)=nanste(NormAccDSAlcohol(selection3,i),1); %SEM
    AvgNormPavAlcohol(1,i)=nanmean(NormAccPavAlcohol(selection4,i)); %average norm accuracy DS
    SEMNormPavAlcohol(1,i)=nanste(NormAccPavAlcohol(selection4,i),1); %SEM
    
    %NormSig(1,i) = ttest2(NormAccDS(:,i),NormAccPav(:,i),'Vartype','unequal'); %welch's t-test
    [~,NormSig(1,i)] = ranksum(NormAccDSAlcohol(selection3,i),NormAccPavAlcohol(selection4,i)); %wilcoxon's ranksum test
end

%prepare shading
upSEMDS=AvgAccDSAlcohol+SEMAccDSAlcohol;
downSEMDS=AvgAccDSAlcohol-SEMAccDSAlcohol;
upSEMPav=AvgAccPavAlcohol+SEMAccPavAlcohol;
downSEMPav=AvgAccPavAlcohol-SEMAccPavAlcohol;
upSEMDSShuff=AvgAccDSShuffAlcohol+SEMAccDSShuffAlcohol;
downSEMDSShuff=AvgAccDSShuffAlcohol-SEMAccDSShuffAlcohol;
upSEMPavShuff=AvgAccPavShuffAlcohol+SEMAccPavShuffAlcohol;
downSEMPavShuff=AvgAccPavShuffAlcohol-SEMAccPavShuffAlcohol;
upSEMNormDS=AvgNormDSAlcohol+SEMNormDSAlcohol;
downSEMNormDS=AvgNormDSAlcohol-SEMNormDSAlcohol;
upSEMNormPav=AvgNormPavAlcohol+SEMNormPavAlcohol;
downSEMNormPav=AvgNormPavAlcohol-SEMNormPavAlcohol;

%plotting decoder accuracy over time
subplot(3,4,5); %accumbens
hold on;
plot(xaxis,AvgAccDSAlcohol(1:13),'Color', DSAlcohol,'linewidth',3);
plot(xaxis,AvgAccDSShuffAlcohol(1:13),'Color', DSShuff,'linewidth',3);
%plot(xaxis,DSSigAlcohol-0.51,'*','Color', 'k');
patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DSAlcohol,'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
xlabel('Seconds post cue');
legend('Alcohol DS','Shuffled','Location','northeast');
%title('Alcohol DS single unit decoding accuracy');
axis([-.75 3.15 0.45 0.75]);
ylabel('Single Unit Decoding Accuracy');

subplot(3,4,6); %Pav
hold on;
plot(xaxis,AvgAccPavAlcohol(1:13),'Color', PavAlcohol,'linewidth',3);
plot(xaxis,AvgAccPavShuffAlcohol(1:13),'Color', PavShuff,'linewidth',3);
%plot(xaxis,PavSig-0.51,'*','Color', 'k');
patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],PavAlcohol,'EdgeColor','none');alpha(0.5);
patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
xlabel('Seconds post cue');
%title('Alcohol CS+ single unit decoding accuracy');
legend('Alcohol CS+','Shuffled','Location','northeast');
axis([-.75 3.15 0.45 0.75]);

%plotting cdf plots for DS at time of max predictive power
subplot(3,4,[3 4]);
DSbin=4;
[cdfDS,xDS] = ecdf(UnitMdlAccDSAlcohol(selection3,DSbin));
[cdfPav,xPav] = ecdf(UnitMdlAccPavAlcohol(selection4,DSbin));
[cdfDSsuc,xDSsuc] = ecdf(UnitMdlAccDSSucrose(selection1,DSbin));
[cdfPavsuc,xPavsuc] = ecdf(UnitMdlAccPavSucrose(selection2,DSbin));

%subplot(4,3,[3 6]);
hold on;
plot(xDSsuc,cdfDSsuc,'Color',DS,'linewidth',2);
plot(xPavsuc,cdfPavsuc,'Color',Pav,'linewidth',2);
plot(xDS,cdfDS,'Color',DSAlcohol,'linewidth',2);
plot(xPav,cdfPav,'Color',PavAlcohol,'linewidth',2);
plot([.5 .5], [0 1],'LineStyle',':','Color','k');
plot([AvgAccDSSucrose(1,4) AvgAccDSSucrose(1,4)], [0 1],'Color',DS);
plot([AvgAccDSAlcohol(1,4) AvgAccDSAlcohol(1,4)], [0 1],'Color',DSAlcohol);
plot([AvgAccPavSucrose(1,4) AvgAccPavSucrose(1,4)], [0 1],'Color',Pav);
plot([AvgAccPavAlcohol(1,4) AvgAccPavAlcohol(1,4)], [0 1],'Color',PavAlcohol);
axis([0.25 1 0 1]);
title(['Cumulative Distribution 0 to .3s']);
xlabel('Decoding Accuracy')
ylabel('Proportion of Population')
legend('Sucrose DS','Sucrose CS+','Alcohol DS','Alcohol CS+','Location','best');


%Plotting ensemble decoding

clear all;
DS=[0    0.4470    0.7410];
Pav=[0.8500    0.3250    0.0980];
DSAlcohol=[0.4940    0.1840    0.5560];
PavAlcohol=[0.9290    0.6940    0.1250];
NumNeurons=[1 5 10 25 50 100];

load('PoolMdlAccDSSucrose.mat','PoolMdlAccDSSucrose')
load('PoolMdlAccDSSucroseShuff.mat','PoolMdlAccDSSucroseShuff')
load('PoolMdlAccPavSucrose.mat','PoolMdlAccPavSucrose')
load('PoolMdlAccPavSucroseShuff.mat','PoolMdlAccPavSucroseShuff')
load('PoolMdlAccDSAlcohol.mat','PoolMdlAccDSAlcohol')
load('PoolMdlAccDSAlcoholShuff.mat','PoolMdlAccDSAlcoholShuff')
load('PoolMdlAccPavAlcohol.mat','PoolMdlAccPavAlcohol')
load('PoolMdlAccPavAlcoholShuff.mat','PoolMdlAccPavAlcoholShuff')


for i=1:6    
AvgAccDSSucrose(i)=nanmean(PoolMdlAccDSSucrose{i,1}); %average accuracy
SEMAccDSSucrose(i)=nanste(PoolMdlAccDSSucrose{i,1},1); %SEM

AvgAccPavSucrose(i)=nanmean(PoolMdlAccPavSucrose{i,1}); %average accuracy
SEMAccPavSucrose(i)=nanste(PoolMdlAccPavSucrose{i,1},1); %SEM

AvgAccDSAlcohol(i)=nanmean(PoolMdlAccDSAlcohol{i,1}); %average accuracy
SEMAccDSAlcohol(i)=nanste(PoolMdlAccDSAlcohol{i,1},1); %SEM

AvgAccPavAlcohol(i)=nanmean(PoolMdlAccPavAlcohol{i,1}); %average accuracy
SEMAccPavAlcohol(i)=nanste(PoolMdlAccPavAlcohol{i,1},1); %SEM

AvgAccDSSucroseShuff(i)=nanmean(PoolMdlAccDSSucroseShuff{i,1}); %average accuracy
SEMAccDSSucroseShuff(i)=nanste(PoolMdlAccDSSucroseShuff{i,1},1); %SEM

AvgAccPavSucroseShuff(i)=nanmean(PoolMdlAccPavSucroseShuff{i,1}); %average accuracy
SEMAccPavSucroseShuff(i)=nanste(PoolMdlAccPavSucroseShuff{i,1},1); %SEM

AvgAccDSAlcoholShuff(i)=nanmean(PoolMdlAccDSAlcoholShuff{i,1}); %average accuracy
SEMAccDSAlcoholShuff(i)=nanste(PoolMdlAccDSAlcoholShuff{i,1},1); %SEM

AvgAccPavAlcoholShuff(i)=nanmean(PoolMdlAccPavAlcoholShuff{i,1}); %average accuracy
SEMAccPavAlcoholShuff(i)=nanste(PoolMdlAccPavAlcoholShuff{i,1},1); %SEM
end



subplot(3,4,[7 8]);
errorbar(NumNeurons, AvgAccDSSucrose,SEMAccDSSucrose,'Color',DS,'linewidth',2);
hold on;
errorbar(NumNeurons, AvgAccPavSucrose,SEMAccPavSucrose, 'Color',Pav,'linewidth',2);
errorbar(NumNeurons, AvgAccDSAlcohol,SEMAccDSAlcohol, 'Color',DSAlcohol,'linewidth',2);
errorbar(NumNeurons, AvgAccPavAlcohol,SEMAccPavAlcohol, 'Color',PavAlcohol,'linewidth',2);
% errorbar(NumNeurons, AvgAccDSSucroseShuff,SEMAccDSSucroseShuff, 'Color','k');
% errorbar(NumNeurons, AvgAccPavSucroseShuff,SEMAccPavSucroseShuff, 'Color','k');
% errorbar(NumNeurons, AvgAccDSAlcoholShuff,SEMAccDSAlcoholShuff, 'Color','k');
% errorbar(NumNeurons, AvgAccPavAlcoholShuff,SEMAccPavAlcoholShuff, 'Color','k');
legend('Sucrose DS','Sucrose CS+','Alcohol DS','Alcohol CS+','Location','best');
plot([-5 105], [.5 .5],'LineStyle',':','Color','k');
xlabel('Neurons used');
ylabel('Mean Accuracy (+/- SEM)');
xticks([1 5 10 25 50 100])
%xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'})
axis([-5 105 0.4 1]); 

%% STATS

%Plotting single All decoding
bins = 1; %number of bins
binint = 0.3; %spacing of bins
binstart = 0; %start time of first bin relative to event
Dura=[0 0.3]; %size of bin used for decoding:
NumNeurons=[1 5 10 25 50 100]; %matrix of how many neurons to use on each iteration
repetitions=20; %how many times to run the analysis

%make matrices that will include best bin at all levels
A=[];
B=[];
C=[];
D=[];
A1=[];
B1=[];
C1=[];
D1=[];


for v=1:length(PoolMdlAccDSSucrose)
A=cat(1,A,PoolMdlAccDSSucrose{v,1});
B=cat(1,B,PoolMdlAccPavSucrose{v,1});
C=cat(1,C,PoolMdlAccDSAlcohol{v,1});
D=cat(1,D,PoolMdlAccPavAlcohol{v,1});
A1=cat(1,A1,PoolMdlAccDSSucroseShuff{v,1});
B1=cat(1,B1,PoolMdlAccPavSucroseShuff{v,1});
C1=cat(1,C1,PoolMdlAccDSAlcoholShuff{v,1});
D1=cat(1,D1,PoolMdlAccPavAlcoholShuff{v,1});
end


%ANOVA for effect of ensemble size, task and reward on accuracy 0-.3 sec
size=[];
for i=1:length(NumNeurons)
    size=cat(1,size,i*ones(repetitions,1));
end
size2=cat(1,size,size,size,size);
%shuffd=cat(1,zeros(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
reward=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1));
task=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{reward,task,size2},'varnames',{'reward','task','ens size'},'model','full');


rewardtask=cat(1,zeros(length(A),1),ones(length(B),1),2*ones(length(C),1),3*ones(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{rewardtask,size2},'varnames',{'reward and task','ens size'},'model','full');


size=[];
for i=1:length(NumNeurons)
    size=cat(1,size,i*ones(repetitions,1));
end
size1=cat(1,size,size,size,size,size,size,size,size);
shuffd=cat(1,zeros(length(A)+length(B)+length(C)+length(D),1),ones(length(A)+length(B)+length(C)+length(D),1));
reward=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1),ones(length(A)+length(B),1),zeros(length(C)+length(D),1));
task=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1),ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D,A1,B1,C1,D1),{reward,task,size1,shuffd},'varnames',{'reward','task','ens size','shuffled vs real'},'model','full');


rewardtask=cat(1,zeros(length(A),1),ones(length(B),1),2*ones(length(C),1),3*ones(length(D),1));
[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{rewardtask,size2},'varnames',{'reward and task','ens size'},'model','full');