%% Analysis script for Jaffe and Brainard, eLife (2020)
% This script analyzes the data and makes the components of Figure 1.
% Code for the other analyses is similar. 

%% Setup
all_features = {'pitch','pitch_cv','seq_length','amp'};
feature_labels = {'pitch','pitch c.v.','syllable sequence duration','amplitude'};
all_conditions = {'carb','saline','carb_atrp','carb_mec_mla'};
condition_labels = {'carb','saline','carb+atrp','carb+mec+mla'};

% Load data
data_path = '';
cd(data_path);
load('Figure1Data2.mat');

% Setup figure
figure_1 = figure; hold on
scatter_jitter = 0.04;
label_font_size = 16;
axis_font_size = 12;
edge_alpha = 0.4;
bar_colors = {'r','b','g','y'};

%% Main
for f = 1:length(all_features)
    this_feature = all_features{f};
    this_label = feature_labels{f};
    disp('---------------------------------------------')
    disp(['mean +/- SEM percent change, ',this_label,':'])   
    subplot(1,length(all_features),f); hold on
    for c = 1:length(all_conditions)
        this_condition = all_conditions{c};
        this_vals = data.(this_feature).(this_condition).normed_vals;
        this_n = length(this_vals);
        this_mean = mean(this_vals);
        this_mean_perc = (this_mean-1)*100;
        this_SEM = std(this_vals)/sqrt(this_n);
        this_SEM_perc = this_SEM*100;
        % Display summary statistics
        disp([this_condition,': ',num2str(this_mean_perc),' +/- ',num2str(this_SEM_perc)])
        % Plot
        bar(c,this_mean,'BaseValue',1,'FaceColor',bar_colors{c})
        if this_mean < 1
            er = errorbar(c,this_mean,this_SEM,[],'Color','k');
        else
            er = errorbar(c,this_mean,[],this_SEM,'Color','k');
        end
        x_scatter = normrnd(c,scatter_jitter,1,this_n);
        scatter(x_scatter,this_vals,'ko','MarkerEdgeAlpha',edge_alpha)
    end %c    
    
    % Adjust plot  
    set(gca,'FontSize',axis_font_size)
    ylabel(this_label,'FontSize',label_font_size)
    set(gca,'XTick',1:length(all_conditions))
    set(gca,'XTickLabel',condition_labels)
    xtickangle(45)   
    xlim([0,5])
    
    % carb vs. saline hypothesis test
    carb_vals = data.(this_feature).carb.normed_vals;
    this_n = length(carb_vals);
    sal_vals = data.(this_feature).saline.normed_vals;
    p_carb_vs_sal = signrank(carb_vals,sal_vals); 
    disp(['carb vs. saline, p = ',num2str(p_carb_vs_sal),', n = ',num2str(this_n),' syllables'])  
    % carb+atrp vs. carb+mec+mla hypothesis test
    carb_atrp_vals = data.(this_feature).carb_atrp.normed_vals;
    this_n = length(carb_atrp_vals);
    carb_mec_mla_vals = data.(this_feature).carb_mec_mla.normed_vals;
    p_carb_atrp_vs_carb_mec_mla = signrank(carb_atrp_vals,carb_mec_mla_vals); 
    disp(['carb+atrp vs. carb+mec+mla, p = ',num2str(p_carb_atrp_vs_carb_mec_mla),', n = ',num2str(this_n),' syllables'])
end %f
