%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).
%% given an eta squared, 
%single variable, change df (df for errors) for between or within design 
close all
clear all
k=4
Sub_NO=(4:1:54)'
%Sub_NO=(20:1:20)'

%between subject design
% df=k*(Sub_NO-1);
% dfeffect=k-1

%within subject design
withinORbetween=0; %0 for within, 1 for between

df=(k-1+withinORbetween)*(Sub_NO-1);
dfeffect=k-1

etaTrue=0.058823529 %given an eta squared=0.11 we can change it to any priori value, 0.0544052196291245 is used in relation between F, et a, and cohneds d.excl



critical_F=finv(0.95,dfeffect,df);

TrueF=(etaTrue./(1-etaTrue)).*(df/dfeffect) % which should not affect power analyses
%true_lambda=TrueF*dfeffect
true_lambda=(etaTrue./(1-etaTrue)).*(df) 
pft=1-ncfcdf(critical_F,dfeffect,df,true_lambda)%pf=observed power= .085
%[Sub_NO, df, pft]
[Sub_NO, df,critical_F, TrueF, true_lambda, pft]
% from the population

plot(Sub_NO,pft,'-')
hold on

figure

%visualize critical F under null (unconcentral F) and alternative
%(concentral F) only for one N

I=find(Sub_NO== 25);
x=(0:0.1:20)';
Lf0=ncfpdf(x,dfeffect,df(I),0);

Lf1=ncfpdf(x,dfeffect,df(I),true_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)


%% given an eta sqaured and df for effect, 
%two way between design which means both variables are between 
close all
clear all

%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)

Sub_NO=(2:1:40)'
groupA=4
groupB=3
effectOrInteraction=2 %1 for 'effect', 2 for interaction
etaTrue=0.0531206757765412 %from a spefic effect
%finGpower=etaTrue/(1-etaTrue)

%


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

df=groupA*groupB*(Sub_NO-1)

critical_F=finv(0.95,dfeffect,df);

TrueF=(etaTrue./(1-etaTrue)).*(df/dfeffect) % which should not affect power analyses
%true_lambda=TrueF*dfeffect
true_lambda=(etaTrue./(1-etaTrue)).*(df) 
pft=1-ncfcdf(critical_F,dfeffect,df,true_lambda)%pf=observed power= .085
%[Sub_NO, df, pft]
[Sub_NO, df,critical_F, TrueF, true_lambda, pft]
% from the population

plot(Sub_NO,pft,'-')

%visualize critical F under null (unconcentral F) and alternative
%(concentral F) only for one N

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

Lf1=ncfpdf(x,dfeffect,df(I),true_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)


%% %% given an eta sqaured and df for effect, 
%within design which means all variables are within if there are multiple variables
close all
clear all

%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

Sub_NO=(6:1:95)'
measurmentA=4 
measurmentB=2
effectOrInteraction=2 %1 for 'effect', 2 for interaction
%
etaTrue=0.058823529 %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

df=dfeffect*(Sub_NO-1)

critical_F=finv(0.95,dfeffect,df);

TrueF=(etaTrue./(1-etaTrue)).*(df/dfeffect) % which should not affect power analyses
%true_lambda=TrueF*dfeffect
true_lambda=(etaTrue./(1-etaTrue)).*(df) 
pft=1-ncfcdf(critical_F,dfeffect,df,true_lambda)%pf=observed power= .085
%[Sub_NO, df, pft]
[Sub_NO, df,critical_F, TrueF, true_lambda, pft]
% from the population

plot(Sub_NO,pft,'-')
figure

%visualize critical F under null (unconcentral F) and alternative
%(concentral F) only for one N

I=find(Sub_NO== 34);
x=(0:0.1:20)';
Lf0=ncfpdf(x,dfeffect,df(I),0);

Lf1=ncfpdf(x,dfeffect,df(I),true_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)



%% given an eta sqaured and df for effect, 
%mixed design one between and one within design  tested
close all
clear all

%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)

Sub_NO=(3:1:40)'
groupA=2 
measuremntB=2
effectOrInteraction=2 %0 for between main, %1 for 'within effect', 2 for interaction
%df=groupA*groupB*(Sub_NO-1)
etaTrue=0.268 %from observed partial eta squared for any specifc effect, change effectOrInteraction as well according to the effect
%finGpower=etaTrue/(1-etaTrue)
%


if effectOrInteraction==0 % between main
    dfeffect=groupA-1
    df=groupA*(Sub_NO-1)
elseif effectOrInteraction==1 % within main
    dfeffect=(measuremntB -1)
    df=groupA*(measuremntB -1)*(Sub_NO-1)
else % interaction
    dfeffect=(groupA -1)*(measuremntB -1)
    df=groupA*(measuremntB -1)*(Sub_NO-1)
end




critical_F=finv(0.95,dfeffect,df);

TrueF=(etaTrue./(1-etaTrue)).*(df/dfeffect) % which should not affect power analyses
%true_lambda=TrueF*dfeffect
true_lambda=(etaTrue./(1-etaTrue)).*(df) 
pft=1-ncfcdf(critical_F,dfeffect,df,true_lambda)%pf=observed power= .085
%[Sub_NO, df, pft]
[Sub_NO, df,critical_F, TrueF, true_lambda, pft]
% from the population

plot(Sub_NO,pft,'-')

%visualize critical F under null (unconcentral F) and alternative
%(concentral F) only for one N

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

Lf1=ncfpdf(x,dfeffect,df(I),true_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',0.1)