mpib_memoreye/gaze_shift/MS_group_analysis.m

clear all;

mainfolder='...' % select the main folder

addpath(genpath([mainfolder '...'])) % add the path for the Circular Stats toolbox

load([mainfolder '../good_partis.mat']) %list of participant with a good performance

partis=good_partis';
partis_sub=good_partis';

%Folders of data to use

config{1}.toloadfolder='/lindedomingo/Memoreye_DATA/eyelink_preprocessed/object_1_onset_11_seconds/';
config{2}.toloadfolder='/lindedomingo/Memoreye_DATA/eyelink_preprocessed/object_1_onset_11_seconds/';
config{3}.toloadfolder='/lindedomingo/Memoreye_DATA/eyelink_preprocessed/object_1_onset_11_seconds/';
config{4}.toloadfolder='/lindedomingo/Memoreye_DATA/eyelink_preprocessed/object_1_onset_11_seconds/';
config{5}.toloadfolder='/lindedomingo/Memoreye_DATA/eyelink_preprocessed/2nd_cue_onset_17_5__seconds/';
config{6}.toloadfolder='/lindedomingo/Memoreye_DATA/eyelink_preprocessed/2nd_cue_onset_17_5__seconds/';

config{1}.timofi= [300 1000];%[100 500]; %first object
config{2}.timofi= [1300 2000]%; %second object
config{3}.timofi= [2350 5850]%; %first delay
config{4}.timofi= [2350 5850]%; %first delay
config{5}.timofi= [650 2850]%[350 2850]; %second delay
config{6}.timofi= [650 2850]%[350 2850]; %second delay

config{1}.object=1;
config{2}.object=2;
config{3}.object=1;
config{4}.object=2;
config{5}.object=1;
config{6}.object=2;

config{1}.retrocue=0;
config{2}.retrocue=0;
config{3}.retrocue=1;
config{4}.retrocue=2;
config{5}.retrocue=2;
config{6}.retrocue=1;

%% Creating empty variables
all_direction=[];
all_direction_size=[];
for rot=1:16
    for conf=1:6
directionPerRot{conf,rot}=[];
directionPerRotAve{conf,rot}=[];
directionPerRotAll{conf}=[];
RotationAll{conf}=[];
sizePerRot{conf,rot}=[];

    end
end

for ppp=1:length(partis)   
    for rot=1:16
        for conf=1:6
            
            directionPerRotPerPar{conf,rot,ppp}=[];
            sizePerRotPerPar{conf,rot,ppp}=[];
            RotationAllPerPar{conf,ppp}=[];
            directionPerRotAllPerPar{conf,ppp}=[];
        end
    end
end

for conf= 1:6 % for each possible configuration pattern
    
    for ppp=1:length(partis)
        
        toload =[mainfolder config{conf}.toloadfolder partis{ppp,:} '.mat'];
        load(toload)
        
        %Loading logfiles
        
        subtable=[];
        
        if strcmp( partis_sub{ppp,:},'p11')
            load([mainfolder 'lindedomingo/Memoreye_DATA/logfiles/resultfile_p11_table_repaired.mat']);
            subtable=sub;
        else
            
            sub = tdfread([mainfolder 'lindedomingo/Memoreye_DATA/logfiles/resultfile_' partis_sub{ppp,:} '.txt'],'tab'); %Logfile for this participant
            subtable=struct2table(sub);
        end
        
        %checking when is time window of interest
        
        tf1=find(data_epoch_recenter.time==config{conf}.timofi(1));
        tf2=find(data_epoch_recenter.time==config{conf}.timofi(2));
        
        time=(tf1:tf2);
        
        
        for ttt=1:size(data_epoch_recenter.eyedat,1)
            
            inputA=squeeze(data_epoch_recenter.eyedat(ttt,1:2,time));
             [num_shifts, shift_directions, sizeG] = detect_gaze_shifts_jld_size(inputA);
%             [shift_directions sizes] = detectMS_EK03_jld(inputA, 1);
%            shift_directions=deg2rad(shift_directions);%only when using MS_EK03
            
            all_direction=[all_direction shift_directions(:)']; %independently of the rotation
            all_direction_size = [all_direction_size, sizeG];
        end
        
        rotlist=unique(data_epoch_recenter.trialinfo.object_1_rot); %list of rotation
        
        for k =1:length(rotlist) % loop - rotations
                        
            %trials of interest            
            if config{conf}.object==1
                
                if config{conf}.retrocue==0;
                    triales = find(data_epoch_recenter.trialinfo.object_1_rot==rotlist(k)...
                        );%objlistid(i));
                else
                    triales = find(data_epoch_recenter.trialinfo.object_1_rot==rotlist(k)...
                        &data_epoch_recenter.trialinfo.retrocue==config{conf}.retrocue);%objlistid(i));
                end
                
            elseif config{conf}.object==2
                
                if config{conf}.retrocue==0;
                    triales = find(data_epoch_recenter.trialinfo.object_2_rot==rotlist(k)...
                        );%objlistid(i));
                else
                    triales = find(data_epoch_recenter.trialinfo.object_2_rot==rotlist(k)...
                        &data_epoch_recenter.trialinfo.retrocue==config{conf}.retrocue);%objlistid(i));
                end
                
            end
            
            forAve=[];
            for ttt=1:length(triales)
                
                inputA=squeeze(data_epoch_recenter.eyedat(triales(ttt),1:2,time));
                [num_shifts, shift_directions, sizeG] = detect_gaze_shifts_jld_size(inputA);
%                 [shift_directions sizes] = detectMS_EK03_jld(inputA, 1);
%                 shift_directions=deg2rad(shift_directions');%only when using MS_EK03


                
                directionPerRot{conf,k}=[directionPerRot{conf,k} shift_directions]; %independently of the rotation
                sizePerRot{conf,k}=[sizePerRot{conf,k} sizeG]; %independently of the rotation
                
                %Per participant
                directionPerRotPerPar{conf,k,ppp}=[directionPerRotPerPar{conf,k,ppp} shift_directions]; %independently of the rotation
                sizePerRotPerPar{conf,k,ppp}=[sizePerRotPerPar{conf,k,ppp} sizeG]; %independently of the rotation
                RotationAllPerPar{conf,ppp}=[RotationAllPerPar{conf,ppp} ones(size(shift_directions)).*rotlist(k)];
                directionPerRotAllPerPar{conf,ppp}=[directionPerRotAllPerPar{conf,ppp} shift_directions];


                forAve=[forAve circ_mean(shift_directions,[],2)];
                directionPerRotAll{conf}=[directionPerRotAll{conf} shift_directions];
                RotationAll{conf}=[RotationAll{conf} ones(size(shift_directions)).*rotlist(k)];


            end
            
            directionPerRotAve{conf,k}=[directionPerRotAve{conf,k} circ_mean(forAve,[],2)];
            
        end
        
    end
    
end

%%
figure;
polarhistogram(all_direction,50)

%% Size of MS

figure;
histogram(all_direction_size.*0.0217)
xlim([0,10])
xlabel('° visual angle')
ylabel('count')

%% %% Combining orientations and cue 1 and 2 for both delays (MS size)
for c=1:6; %Creating an empty variable
sizePerRotCom{4}=[];
end
for i=1:16
sizePerRotCom{1}=[sizePerRotCom{1} sizePerRot{1,i}];
sizePerRotCom{2}=[sizePerRotCom{2} sizePerRot{2,i}];
sizePerRotCom{3}=[sizePerRotCom{3} [sizePerRot{3,i} sizePerRot{4,i}]];
sizePerRotCom{4}=[sizePerRotCom{4} [sizePerRot{5,i} sizePerRot{6,i}]];
end

%%
listOfTil={'stim1','stim2','delay1','delay2'}

figure;

for cond =1:4;
    
    subplot(1,4,cond);

histogram(sizePerRotCom{cond}.*0.0217,500)
xlim([0,10])
xlabel('° visual angle')
ylabel('count')

title([listOfTil{cond}])
end

%% Combining cue 1 and 2 for both delays
for i=1:16
directionPerRotCom{1,i}=directionPerRot{1,i};
directionPerRotCom{2,i}=directionPerRot{2,i};
directionPerRotCom{3,i}=[directionPerRot{3,i} directionPerRot{4,i}];
directionPerRotCom{4,i}=[directionPerRot{5,i} directionPerRot{6,i}];
end

RotationAllCom{1}=RotationAll{1};
RotationAllCom{2}=RotationAll{2};
RotationAllCom{3}=[RotationAll{3} RotationAll{4}];
RotationAllCom{4}=[RotationAll{5} RotationAll{6}];

directionPerRotAllCom{1}=directionPerRotAll{1};
directionPerRotAllCom{2}=directionPerRotAll{2};
directionPerRotAllCom{3}=[directionPerRotAll{3} directionPerRotAll{4}];
directionPerRotAllCom{4}=[directionPerRotAll{5} directionPerRotAll{6}];

%% Same combining per participant
for ppp = 1:length(partis)
    for i=1:16
        directionPerRotComPerPar{1,i,ppp}=directionPerRotPerPar{1,i,ppp};
        directionPerRotComPerPar{2,i,ppp}=directionPerRotPerPar{2,i,ppp};
        directionPerRotComPerPar{3,i,ppp}=[directionPerRotPerPar{3,i,ppp} directionPerRotPerPar{4,i,ppp}];
        directionPerRotComPerPar{4,i,ppp}=[directionPerRotPerPar{5,i,ppp} directionPerRotPerPar{6,i,ppp}];
    end
end

for ppp = 1:length(partis)
    
    RotationAllComPerPar{1,ppp}=RotationAllPerPar{1,ppp};
    RotationAllComPerPar{2,ppp}=RotationAllPerPar{2,ppp};
    RotationAllComPerPar{3,ppp}=[RotationAllPerPar{3,ppp} RotationAllPerPar{4,ppp}];
    RotationAllComPerPar{4,ppp}=[RotationAllPerPar{5,ppp} RotationAllPerPar{6,ppp}];
    
    directionPerRotAllComPerPar{1,ppp}=directionPerRotAllPerPar{1,ppp};
    directionPerRotAllComPerPar{2,ppp}=directionPerRotAllPerPar{2,ppp};
    directionPerRotAllComPerPar{3,ppp}=[directionPerRotAllPerPar{3,ppp} directionPerRotAllPerPar{4,ppp}];
    directionPerRotAllComPerPar{4,ppp}=[directionPerRotAllPerPar{5,ppp} directionPerRotAllPerPar{6,ppp}];
end


%% Plotting MS per condition

colores=(colormap(hsv(16)));

listOfTil={'stim1','stim2','delay1','delay2'}
figure;

for conf=1:4;
subplot(1,4,conf);
for i=1:16

polarhistogram(directionPerRotCom{conf,i}+pi,40,'FaceColor',colores(i,:),'FaceAlpha',.5);hold on;

end
a=RotationAllCom{conf};
b=directionPerRotAllCom{conf};

[rho(conf) pval(conf)]=circ_corrcc(deg2rad(a),b);
[rho180(conf) pval180(conf)]=circ_corrcc((deg2rad(a)+pi).*2,(b).*2)

title([listOfTil{conf} ' pval360=' num2str(pval(conf)) ' pval180='  num2str(pval180(conf))])

end

%% Stats
Stats180.rho=rho180'
Stats180.pval=pval180'
Stats180t=struct2table(Stats180)

Stats360.rho=rho'
Stats360.pval=pval'
Stats360t=struct2table(Stats360)

%% Group Rho  vs. 0

for ppp = 1:length(partis);
    for conf=1:4;
        
        a=RotationAllComPerPar{conf,ppp};
        b=directionPerRotAllComPerPar{conf,ppp};
        
        [rho(conf,ppp) pval(conf,ppp)]=circ_corrcc(deg2rad(a),b);
        [rho180(conf,ppp) pval180(conf,ppp)]=circ_corrcc((deg2rad(a)+pi).*2,(b).*2);
        
    end
end

[h,p,ci,stats] = ttest(rho', 0)
[h,p,ci,stats] = ttest(rho180', 0)


%% Re-rotated comb

listOfTil={'stim1','stim2','delay1','delay2'}
figure;

for conf=1:4;
subplot(1,4,conf);

a=deg2rad(RotationAllCom{conf});
b=directionPerRotAllCom{conf};
c=angdiff(b,a)
polarhistogram((c),40,'FaceColor',colores(i,:),'FaceAlpha',.5);hold on;

[Rpval(conf) Rz(conf)] = circ_rtest((c))

title([listOfTil{conf}, ' z=' num2str(Rz(conf)) ' p=' num2str(Rpval(conf))])
end

StatsR.z=Rz'
StatsR.pval=Rpval'
StatsRt=struct2table(StatsR)

%% %% Group Rho  vs. 0

for ppp = 1:length(partis);
    for conf=1:4;
        
        a=RotationAllComPerPar{conf,ppp};
        b=directionPerRotAllComPerPar{conf,ppp};
        c=angdiff(b,a)

        [Rpval(conf,ppp) Rz(conf,ppp)] = circ_rtest((c))
        
    end
end

[h,p,ci,stats] = ttest(Rpval', 0)

%% Re-rotated 

listOfTil={'stim1','stim2','delay1-cue1','delay1-cue2','delay2-cue1','delay2-cue2',}
figure;

for conf=1:6;
subplot(3,2,conf);

a=deg2rad(RotationAll{conf});
b=(directionPerRotAll{conf});
c=angdiff(b,a);
polarhistogram((c),40,'FaceColor',colores(i,:),'FaceAlpha',.5);hold on;

[pval z] = circ_rtest((c))
degmean=rad2deg(circ_mean((c),[],2));

title([listOfTil{conf}, ' z=' num2str(z) ' p=' num2str(pval) ' mean=' num2str(degmean)])

end