powerMAOV1

function [Power] = powerMAOV1
%  Procedures to Estimate the Statistical Power of a Performed (a posteriori) 
%  Single-Factor Multivariate Analysis of Variance.
%
%   Syntax: function [Power] = powerMAOV1 
%      
%     Inputs:
%      In a prompt dialogue you must to choose among 2 options
%      Option 1
%          F - multivariate observed F-statistic.
%         v1 - numerator degrees of freedom.
%         v2 - denominator degrees of freedom.
%      alpha - significance (default = 0.05).
%
%      Option 2
%          L - observed Wilks' lambda.
%          p - number of variables.
%          k - number of samples (groups).
%          N - total number of data.
%      alpha - significance (default = 0.05).
%
%       [You must to enter the parameters among brackets (vectors)]
%
%     Output:
%          Power - the output power of the performed MANOVA test.
%  
%    Example: For an unequal sample size multivariate design with a 
%             Wilks'L (lambda) = 0.0385, p = 2, k = 3, N = 8, and 
%             alpha = 0.05. We are interested to estimate the statistical
%             power in the performed MANOVA test. 
%
%     Calling on Matlab the function: 
%             powerMAOV1
%
%  Immediately it gives you some statements as:
%    -For the performed Multivariate Analysis of Variance, this file has
%    2 options to evaluate the statistical power.
%
%    -For choose option 1 you must to input vector:
%    [F-multivariate observed statistic, v1-numerator df, v2-denominator df,
%    and alpha-significance]
% 
%    -For choose option 2 you must to input vector:
%    [Wilks'L(lambda), p-variables, k-samples(groups), N-total of data, and
%    alpha-significance]
%
%       Then it asks:
% 
%    -Do you need option 1 or option 2?: 2
%    -Give me the vector of parameters:[.0385,2,3,8]
%
%              Answer is:
%
%                ans= 0.9960
%           

%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             [e-mail]
%
%  August 14, 2003.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). powerMAOV1: Estimation of statistical power of a performed 
%    single-factor AMOVA. A MATLAB file. [WWW document]. URL http://www.mathworks.com/matlabcentral/fileexchange/
%    loadFile.do?objectId=3864&objectType=FILE
%
%  References:
%  Cohen, J. (1977), Statistical Power Analysis for the Behavioral Sciences.
%              New York:Academic Press. Chapter 10.
%  Trujillo-Ortiz, A. On the Statistical Power of One-Way 
%              Analysis of Variance Test Model I. The American
%              Statistician. (Submitted).
%  Winer, B.J. (1971), Statistical Concepts in Experimental Design.
%              New York:McGraw-Hill. pp. 220-222; 225-228.
%

disp('For the performed Multivariate Analysis of Variance, this file gives you');
disp('2 options to evaluate the statistical power.');
disp(' ')
disp('For choose option 1 you must to input vector:');
disp('[F-multivariate observed statistic, v1-numerator df, v2-denominator df,');
disp('and alpha-significance]');
disp(' ')
disp('For choose option 2 you must to input vector:');
disp('[Wilks''L(lambda), p-variables, k-samples(groups), N-total of data, and'); 
disp('alpha-significance]');
disp(' ')
o = input('Do you need option 1 or option 2?: ');
if o == 1;
   r = input('Give me the vector of parameters:');
   nr = length(r);
   if nr == 4;
      F = r(1);
      v1 = r(2);
      v2 = r(3);
      alpha = r(4);
      if (alpha <= 0 | alpha >= 1)
         fprintf('Warning: significance must be between 0 and 1\n');
         return;
      end;
   else nr ~= 4;
      if nr == 3;
         F = r(1);
         v1 = r(2);
         v2 = r(3);
         alpha = 0.05;
      else nr < 3;
         error('Requires at least tree input arguments.');
      end;
   end;
   f2 = F*v1/v2;  %effect size index
   l = f2*(v1 + v2 + 1);  %estimated noncentrality parameter
else (o == 2);
   r = input('Give me the vector of parameters:');
   nr = length(r);
   if nr == 5;
      L = r(1);
      p = r(2);
      k = r(3);
      N = r(4);
      alpha = r(5);
      if (alpha <= 0 | alpha >= 1)
         fprintf('Warning: significance must be between 0 and 1\n');
         return;
      end;
   else nr ~= 5;
      if nr == 4;
         L = r(1);
         p = r(2);
         k = r(3);
         N = r(4);
         alpha = 0.05;
      else nr < 4;
         error('Requires at least four input arguments.');
      end;
   end;
   v1 = p*(k - 1);  %numerator degrees of freedom
   if p == 2 | (k-1) == 2
       s = 2;
   else
       s = sqrt((p^2*(k-1)^2-4)/(p^2+((k-1)^2)-5));
   end;
   v2 = (s*((N-1)-((p+k)/2)))-(((p*(k-1))-2)/2);  %denominator degrees of freedom
   v=L^(1/s);
   F=((1-v)/v)*(v2/v1);  %approximation to Rao's F-statistic
   f2 = F*v1/v2;  %effect size index
   l = f2*(v1 + v2 + 1);  %estimated noncentrality parameter
end;
%
% Aproximation of the noncentral F distribution to central F distribution.
%
P = 1 - alpha;
Fc=finv(P,v1,v2)/(1+(l/v1)); %Expected F-statistic value adjusted to the estimated
%noncentrality parameter.
v1m=((v1+l)^2)/(v1+(2*l)); %Numerator degrees of freedom adjusted to the estimated
%noncentrality parameter.
%
% Because the numerator degrees of freedom corrected by the noncentrality parameter
% could results a fraction, the probability function associated to the F distribution
% function is resolved by the Simpson's 1/3 numerical integration method.
%
x=linspace(.00001,Fc,10001);
DF=x(2)-x(1);
y=((v1m/v2)^(.5*v1m)/(beta((.5*v1m),(.5*v2))));
y=y*(x.^((.5*v1m)-1)).*(((x.*(v1m/v2))+1).^(-.5*(v1m+v2)));
N1=length(x);
Power=1-(DF.*(y(1)+y(N1) + 4*sum(y(2:2:N1-1))+2*sum(y(3:2:N1-2)))/3.0);