%reference  Keppel, G. & Wickens, T. D.,(2004). Design and analysis: A researcher's handbook. Pearson Prentice-Hall.
%Cohen, J. (1977). Statistical power for the behavioral sciences (revised edition).
classdef collectivePowerAnalysis
    
    methods
        function  [n, df,criticalT, lambda, p2]=powerCohensd(self,d, n,withinOrBetween,tail)%observedD=observedt*sqrt(2/N), which is sqrt(2)*dz for within design; withinOrBetween=1 for within, 2 for between
            %d=observedD;

            lambda=d*sqrt(n/2); %lambda is observedt
            df=withinOrBetween*(n-1); % for within design

            if tail==1
                criticalT=tinv(0.95,df); %two tailed
                p2=1-nctcdf(criticalT,df,lambda) %for one tailed
            elseif tail==2
                criticalT=tinv(0.975,df); %two tailed
                p2=nctcdf(-criticalT,df,lambda)+1-nctcdf(criticalT,df,lambda) %the percentage x>criticalT or x<- criticalT, should be same as the value in the table
            else
                fprintf("tail should 1 or 2");
            end

            [n, df,criticalT, lambda, p2]
            
            showFigure=0;
            if showFigure==1
                figure 
                plot(n,p2,'-')
            end
        end
        
        function [n, dferror,critical_F, F, lambda, p2]= eta_singleVariable(self,etaSqr, k, n, withinOrBetween)

% 
%             %between subject design
%             % df=k*(Sub_NO-1);
%             % dfeffect=k-1
% 
%             %within subject design
%             withinORbetween=0; %0 for within, 1 for between
% 
            dferror=(k-2+withinOrBetween)*(n-1);
            dfeffect=k-1;
            critical_F=finv(0.95,dfeffect,dferror);

            F=(etaSqr./(1-etaSqr)).*(dferror/dfeffect); % which should not affect power analyses     
            lambda=F.*dfeffect;%lambda=F*dfeffect
            p2=1-ncfcdf(critical_F,dfeffect,dferror,lambda)%pf=observed power= .085
            %[Sub_NO, df, pft]
            [n, dferror,critical_F, F, lambda, p2]
            % from the population
            
            showFigure=0;
            if showFigure==1
                plot(n,p2,'-')

                figure
    
                %visualize critical F under null (unconcentral F) and alternative
                %(concentral F) only for one N
    
                I=find(n== 25);
                x=(0:0.1:20)';
                Lf0=ncfpdf(x,dfeffect,dferror(I),0);
    
                Lf1=ncfpdf(x,dfeffect,dferror(I),lambda(I));
    
                plot (x, Lf0)
                hold on
                plot (x, Lf1)
                hold on
                a=[critical_F(I),critical_F(I)];  b=[0,0.4]; %draw a line on the crtical F
                plot(a,b,'-k','LineWidth',1)
            end 
        end
        
        function [n, dferror,critical_F, F, lambda, p2]= eta_twoBetweenVariables(self,etaSqr, groupA, groupB, n, effectOrInteraction)
            %configure 
            %first, the range of subject numbers, 
            %second, the number of group in each variable, eg 4 and 3 for 4by3 design. the one to be calcualted always is listed as
            %variable A, if we are interested in the power for the varialbe with 3 groups then list groupA=3
            %last whether the interested one is interaction or main effect (i.e. effect A as specified already)
            %groupA=4
            %groupB=3
            %effectOrInteraction=2 %1 for 'effect', 2 for interaction
            %etaSqr=0.0531206757765412 %from a spefic effect
            %finGpower=etaTrue/(1-etaTrue)

            %


            if effectOrInteraction==1
                dfeffect=groupA-1
            else
                dfeffect=(groupA -1)*(groupB -1)
            end

            dferror=groupA*groupB*(n-1);

            critical_F=finv(0.95,dfeffect,dferror);

            F=(etaSqr./(1-etaSqr)).*(dferror/dfeffect);
            %lambda=F*dfeffect
            lambda=(etaSqr./(1-etaSqr)).*(dferror); 
            p2=1-ncfcdf(critical_F,dfeffect,dferror,lambda);%pf=observed power= .085
            %[Sub_NO, df, pft]
            [n, dferror,critical_F, F, lambda, p2]
            % from the population
            
            showFigure=0;
            if showFigure==1
                figure 
                plot(n,p2,'-')
                
                %visualize critical F under null (unconcentral F) and alternative
                %(concentral F) only for one N

                I=find(n== 28);
                x=(0:0.1:20)';
                Lf0=ncfpdf(x,dfeffect,dferror(I),0);

                Lf1=ncfpdf(x,dfeffect,dferror(I),lambda(I));
                figure
                plot (x, Lf0)
                hold on
                plot (x, Lf1)
                hold on
                a=[critical_F(I),critical_F(I)];  b=[0,0.4]; %draw a line on the crtical F
                plot(a,b,'-k','LineWidth',1)
            end
        end   
         
        function [n, dferror,critical_F, F, lambda, p2]= eta_twoWithinVariables(self,etaSqr, measurmentA, measurmentB, n, effectOrInteraction)
            %within design which means all variables are within if there are multiple variables
                        
            
            %configure 
            %first, the range of subject numbers, 
            %second, the number of measurment in each variable, eg 4 and 3 for 4by3 design. the one to be calcualted always is listed as
            %variable A, if we are interested in the power for the varialbe with 3 measurement then list measurmentA=3
            %last whether the interested one is interaction or main effect (i.e. effect A as specified already)
            %Note for all effect, main, interaction, df(error)= df(effect)*(n-1) , same as in one variable within design, 
            %so we can either treat it as one variable and make measurement = df(effect)+1 in the one variable within case 
            %for example for 4by3 design if we are interested in interaction, then we
            %can use measument=(4-1)*(3-1)+1 in the single variable power analysis

            %measurmentA=2 
            %measurmentB=2
            %effectOrInteraction=1 %1 for 'effect', 2 for interaction
            %
            %etaSqr=0.3760 %given an eta squared=0.11 we can change it to any priori value
            %finGpower=sqrt(etaTrue/(1-etaTrue))


            if effectOrInteraction==1
                dfeffect=measurmentA-1;
            else
                dfeffect=(measurmentA -1)*(measurmentB -1);
            end

            dferror=dfeffect*(n-1);

            critical_F=finv(0.95,dfeffect,dferror);

            F=(etaSqr./(1-etaSqr)).*(dferror/dfeffect) % which should not affect power analyses
            %lambda=F*dfeffect
            lambda=(etaSqr./(1-etaSqr)).*(dferror) 
            p2=1-ncfcdf(critical_F,dfeffect,dferror,lambda)%pf=observed power= .085
            %[Sub_NO, df, pft]
            [n, dferror,critical_F, F, lambda, p2]
            
            showFigure=1;
            if showFigure==1
                figure 
                plot(n,p2,'-')
                
                %visualize critical F under null (unconcentral F) and alternative
                %(concentral F) only for one N

                I=find(n== 28);
                x=(0:0.1:20)';
                Lf0=ncfpdf(x,dfeffect,dferror(I),0);

                Lf1=ncfpdf(x,dfeffect,dferror(I),lambda(I));
                figure
                plot (x, Lf0)
                hold on
                plot (x, Lf1)
                hold on
                a=[critical_F(I),critical_F(I)];  b=[0,0.4]; %draw a line on the crtical F
                plot(a,b,'-k','LineWidth',1)
            end
        end
        
        function [n, dferror,critical_F, F, lambda, p2]= eta_twoMixedVariables(self,etaSqr, groupA, measurementB, n, effectOrInteraction)
            %mixed design one between and one within design  tested
           

            %configure 
            %first, the range of subject numbers, 
            %second, the number of levels in each variable, eg 4 and 3 for 4by3 design. the between one is listed as
            %variable A
            %last whether the interested one is between main, within main or interaction or main effect (i.e. effect A as specified already)

            %n=(4:1:100)'
            %groupA=4 
            %measurementB=3
            %effectOrInteraction=2 %2 for between main, %1 for 'within effect', 3 for interaction
            %etaSqr=0.058823529 %from observed partial eta squared for any specifc effect, change effectOrInteraction as well according to the effect

            %powerAnalysis=collectivePowerAnalysis
            %[n, dferror,critical_F, F, lambda, p2]= powerAnalysis.eta_twoMixedVariables(etaSqr, groupA, measurementB, n, effectOrInteraction);


            if effectOrInteraction==2 % between main
                dfeffect=groupA-1;
                dferror=groupA*(n-1);
            elseif effectOrInteraction==1 % within main
                dfeffect=(measurementB -1);
                dferror=groupA*(measurementB -1)*(n-1);
            elseif effectOrInteraction==3 % interaction
                dfeffect=(groupA -1)*(measurementB -1);
                dferror=groupA*(measurementB -1)*(n-1);
            else
                fprintf("effectOrInteraction=1 for between main, 2 for 'within effect', 3 for interaction");
            end

            critical_F=finv(0.95,dfeffect,dferror);

            F=(etaSqr./(1-etaSqr)).*(dferror/dfeffect) % which should not affect power analyses
            %true_lambda=TrueF*dfeffect
            lambda=(etaSqr./(1-etaSqr)).*(dferror) 
            p2=1-ncfcdf(critical_F,dfeffect,dferror,lambda)%pf=observed power= .085
            %[Sub_NO, df, pft]
            [n, dferror,critical_F, F, lambda, p2]
            
            showFigure=0;
            if showFigure==1
                figure 
                plot(n,p2,'-')
                
                %visualize critical F under null (unconcentral F) and alternative
                %(concentral F) only for one N

                I=find(n== 28);
                x=(0:0.1:20)';
                Lf0=ncfpdf(x,dfeffect,dferror(I),0);

                Lf1=ncfpdf(x,dfeffect,dferror(I),lambda(I));
                figure
                plot (x, Lf0)
                hold on
                plot (x, Lf1)
                hold on
                a=[critical_F(I),critical_F(I)];  b=[0,0.4]; %draw a line on the crtical F
                plot(a,b,'-k','LineWidth',1)
            end
        end

        
    end
end