%Read data associated with Object Bias test, perform stats analysis and generate Figure 3    
% Ying Zheng, 04/03/2024. Version 1.

clear all; close all; home;

OBTdata = xlsread('Data_ObjBiasTest','sheet1','B3:Y10'); %All object bias test data

[m, n]=size(OBTdata); %m=num of rats
Nobj=n/3; % number of object, 3 columns per obj 
for an = 1:m
    id_obj=0;
    for i=1:3:n
        id_obj=id_obj+1;
        T_LR_obj(an,id_obj) = sum(OBTdata(an,i:i+1),2);
        WhichDay(an,id_obj) = OBTdata(an,i+2);
    end
end
%****test for outlier
temp1=filloutliers(T_LR_obj,"previous","mean");  
max(max(abs(temp1-T_LR_obj)))  % zero means no outlier      
%test for normality
for i=1:Nobj
    h1(i)=jbtest(T_LR_obj(:,i)); % h1 zero means normal
end
%******* one way repeated measure anova using RMAOV1
c1=T_LR_obj(:);
c2=ones(8,8)*diag(1:8);c2=c2(:);
v=1:8;c3=repmat(v',8,1);
X=[c1, c2, c3];
RMAOV1(X);  %This function was downloaded from Online Matlab help centre  
            %It was cited in the paper. Web link:
%https://uk.mathworks.com/matlabcentral/fileexchange/5576-rmaov1?s_tid=srchtitle_support_results_1_RMAOV1 

%read data based on which day tested
for an=1:m
    for d=1:Nobj
        idx=find(WhichDay(an,:)==d);
        T_LR_day(an,d) = T_LR_obj(an,idx);
    end
end
%test for outlier
temp2=filloutliers(T_LR_day,"previous","mean");  
max(max(abs(temp2-T_LR_day)))    % zero means no outlier    
%test for normality
for i=1:8
    h2(i)=jbtest(T_LR_day(:,i)); % h2 zero means normal
end
%******* repeated measure anova
c1=T_LR_day(:);
Y=[c1, c2, c3];
RMAOV1(Y);  %downloaded from Online Matlab help centre    

%*******  Confirm the result with anova2, then use multcompare    
[pp,tbl,stats]=anova2(T_LR_day);
C=multcompare(stats,"CriticalValueType","bonferroni");


%plot data for Fig 3 in the paper
figure('pos',[200 200 700 150]);
m_obj=mean(T_LR_obj); sem_obj=std(T_LR_obj)/sqrt(m);
subplot(122);
bar(m_obj,0.7,'FaceColor', [0.8 0.8 0.8]);hold on; errorbar(m_obj,sem_obj,'.r');
axis([0 9 0 90]);xlabel('Object #');ylabel('T_t_o_t [s]');box off
m_day=mean(T_LR_day); sem_day=std(T_LR_day)/sqrt(m);
subplot(121);
bar(m_day,0.7,'FaceColor', [0.8 0.8 0.8]);hold on; errorbar(m_day,sem_day,'.r');
axis([0 9 0 90]);xlabel('Test Day #');ylabel('T_t_o_t [s]');box off
set(gcf,'color','w');