%% Documentation
%run statistical analysis on data including model likelihood, ANOVA and
%t-tests
%%
    subNum=unique(data.subData.sub);
    nSub=length(subNum);
    parentpath = cd(cd('..'));
    statsFile=fullfile(parentpath,'data/subdata/stats/summarystatsData.csv');
    statsFile2=fullfile(parentpath,'data/subdata/stats/difstatsData.csv');
    cond=unique(data.subData.sigXcond);

    % save summary data for double check in SPSS/JASP
    for i= 1:length(subNum)
        statsDat(i,1)= data.subData.wei80pThresh(data.subData.sub==subNum(i)&data.subData.sigXcond==cond(1)&data.subData.dim==0);
        statsDat(i,2)=data.subData.wei80pThresh(data.subData.sub==subNum(i)&data.subData.sigXcond==cond(2)&data.subData.dim==0);
        statsDat(i,3)=data.subData.wei80pThresh(data.subData.sub==subNum(i)&data.subData.sigXcond==cond(3)&data.subData.dim==0);
        statsDat(i,4)=data.subData.wei80pThresh(data.subData.sub==subNum(i)&data.subData.sigXcond==cond(1)&data.subData.dim==1);
        statsDat(i,5)=data.subData.wei80pThresh(data.subData.sub==subNum(i)&data.subData.sigXcond==cond(2)&data.subData.dim==1);
        statsDat(i,6)=data.subData.wei80pThresh(data.subData.sub==subNum(i)&data.subData.sigXcond==cond(3)&data.subData.dim==1);
    end
    T=table(statsDat(:,1),statsDat(:,2),statsDat(:,3),statsDat(:,4),statsDat(:,5),statsDat(:,6));
    T.Properties.VariableNames(1:6) = {'Dim2Dis1' 'Dim2Dis2','Dim2Dis3','Dim3Dis1','Dim3Dis2','Dim3Dis3'};
   
   writetable(T,statsFile)
 

    %% 2x3 repeated measures anova on thresholds
    AnovaAna=rm_anova2(data.subData.wei80pThresh,data.subData.sub,data.subData.sigXcond,data.subData.dim,{'sigX','Dim'});
    fprintf('\n..................\n')
    fprintf('Statistical analysis \n \n')
    disp('results of repeated measures ANOVA are')
    disp(AnovaAna)
    
    %% post hoc paired samples t tests between 2D and 3D
    for i = 1:length(cond)
    %get 2D and 3D data, sort by subject.
    dat3D=[data.subData.sub(data.subData.dim==1&data.subData.sigXcond==cond(i));data.subData.wei80pThresh(data.subData.dim==1&data.subData.sigXcond==cond(i))]';
    dat3D=sortrows(dat3D);
    dat2D=[data.subData.sub(data.subData.dim==0&data.subData.sigXcond==cond(i));data.subData.wei80pThresh(data.subData.dim==0&data.subData.sigXcond==cond(i))]';
    dat2D=sortrows(dat2D);
    [~,P(i),~,STATS{i}]=ttest(dat2D(:,2),dat3D(:,2));
    STATS{i}.p = P(i);
    end
     
    fprintf('Results of paired t-test between 2D and 3D\n');
    fprintf('\nsigX 2\n')
    disp(STATS{1})
    fprintf('\nsigX 3\n')
    disp(STATS{2})
    fprintf('\nsigX 4\n')
    disp(STATS{3})
    fprintf('..................\n')
    
    %% Model likelihood
    fprintf('\n..................\n \n')
    fprintf('Model analysis\n \n')
    modelInf={'FE','NEcd','NEid','NV'};
    for i=1:length(modelInf) % run for each model
        [argval(i),argind(i)]=max(data.model.modelLikeli{i});
    end
    [maxlik,maxlikmodel]= max(argval);
    fprintf('Most likely model is model %s \n',modelInf{maxlikmodel})
    for i=1:length(modelInf) % run for each model
        fprintf('relative likelihood for model %s : %10f \n',modelInf{i},data.model.modelLikeli{i}(argind(i))-maxlik)
    end
    for i=1:length(modelInf) % run for each model
        fprintf('best fit parameters for model %s, c: %9f, d:%10f \n',modelInf{i},data.model.parC(argind(i),i),data.model.parD(argind(i),i))
    end
    
      c=data.model.parC(argind(maxlikmodel),maxlikmodel);
      d=data.model.parD(argind(maxlikmodel),maxlikmodel);
    %% Convert parameters to format of other experiments
    fprintf(' \nConverting best fit parameters from most likely model %s \n',modelInf{maxlikmodel})
    modelParConversion(c,d)
        
        
        
    