function [EpsHF EpsList EpsGG]=GenCalcHFEps(Y,BTFacs,WInFacs,S)
%function [EpsHF EpsList EpsGG]=GenCalcHFEps(Y,BTFacs,WInFacs,S)
%
% This will calculate the Geisser-Greenhouse and Huynh-Feldt epsilon value
% for the general case when given a univariate dataset with any amount of 
% between or within subject factors. 
%
% After calling this function, the intention is to multiply the degrees of
% freedom for tests in a repeated measures ANOVA by the corresponding
% epsilon value given by this code to obtain a corrected p-value. The
% F-statistic does not change as a result of this.
%
% Inputs:   
%           Y-  A column vector of the dependent variable (the value
%               measured) for each data point.
%      BTFacs-  A matrix of between subject factors. Must have the same
%               number of rows as Y. Each factor is down a different column
%               of BTFacs, and the value in each row denotes the level of 
%               that factor for each corresponding datapoint in Y. If there
%               are no between subject factors, input an empty matrix, [].
%     WInFacs-  A matrix of within subject factors. Must have the same
%               number of rows as Y. Each factor is down a different column
%               of WInFacs, and the value in each row denotes the level of 
%               that factor for each corresponding datapoint in Y.
%           S-  A column vector of subject numbers corresponding to each
%               datapoint in Y. Must have the same number of rows as Y. If
%               left empty, the program assumes all subjects are entered in
%               the same order for all combinations of factors.
%
% Outputs:
%       EpsHF-  A row vector of all the possible values for the Huynh-Feldt
%               epsilon, correspondding to all main effects for ecah within
%               subject factor, and all interactions of within subject
%               factors.
%     EpsList-  A cell array of a text list of the effects corresponding to
%               the positions in EpsHF. For example 'A' means a main effect
%               for the first (within subject) factor, 'AB' means a 2way 
%               interaction for the first two (within subject) factors, 
%               etc.
%       EpsGG-  The Geisser-Greenhouse epsilon values, in teh same form as
%               EpsHF.
%
% The Huynh-Feldt epsilon value is less conservative while maintaining the
% proper Type I error rate, which is why the program title focuses on that.
% The Geisser-Greenhouse epsilon is calculated first as part of the 
% process of calculating the Huynh-Feldt value, so I've allowed it to be 
% optionally returned as well. 
%
% This code follows a procedure described in Huynh(1978), but includes a 
% modification of the last step, as described in Chen & Dunlap(1994).
%
% I've tested this with results from SPSS for a couple of sample datasets
% I've found in various websites, and it gives results that match
% According to what I've seen and read though, SPSS and SAS give results
% for EpsHF that are wrong when both between and with-in factors are
% present, and thus they slightly differ from the results this program,
% although the results for EpsGG are identical. To reproduce the SPSS 
% results, uncomment line 153. See other comments near there for references
% on why I think SPSS/SAS is wrong.
%
%
% References:
%
% Huynh H. "Some approximate tests for repeated measurement designs", 
% Psychometrika (1978) 
%
% Chen, RS and Dunlap, WP "A MonteCarlo STudy on the Performance of a 
% Corrected Formula for eps(tilda) suggested by Lecoutre", Journal of
% Educational Statistics (1994)
%
% written on 090202 by Matthew Nelson


if nargin<2 || isempty(BTFacs);    
    BT.nFacs=1;
    BT.nTreats=1;
    g=1;
else
    BT.Facs=BTFacs;
    BT.nFacs=size(BTFacs,2);
    BT.GNums=cell(1,BT.nFacs);
    BT.nTreats=repmat(0,1,BT.nFacs);
    %BT.dfs=BT.nTreats;
    for iBF=1:BT.nFacs
        BT.GNums{iBF}=unique(BT.Facs(:,iBF));
        BT.nTreats(iBF)=length(BT.GNums{iBF});
    end
    %BTdfs=BTnTreats-1; %I don't think this is ever used...
    
    g=prod(BT.nTreats);     
end

if nargin<3 || isempty(WInFacs);
    disp('In GenCalcHFEps- no Within Factors entered; leaving without calculating any epsilons')
    return
else
    WIn.Facs=WInFacs;
    WIn.nFacs=size(WIn.Facs,2);
    WIn.GNums=cell(1,WIn.nFacs);
    WIn.nTreats=repmat(0,1,WIn.nFacs);
    %WIn.dfs=WIn.nTreats;
    WIn.FacMult=WIn.nTreats;
    for iBF=1:WIn.nFacs
        WIn.GNums{iBF}=unique(WIn.Facs(:,iBF));
        WIn.nTreats(iBF)=length(WIn.GNums{iBF});
    end
    %WIn.dfs=WIn.nTreats-1; %I don't think this is ever used...
    
    WIn.nCombs=prod(WIn.nTreats);
    %Calc FactMult for determining curComb later
    for iFac=1:WIn.nFacs-1
        WIn.FacMult(iFac)=prod(WIn.nTreats(iFac+1:end));
    end
    WIn.FacMult(end)=1;
end

if nargin<4 || isempty(S)
    S=[];
    N=length(Y)/ (WIn.nCombs*prod(BT.nTreats));     %if S is not input, we assume that there is data for each subject in each cell... 
else  N=length(unique(S))/prod(BT.nTreats);   %divide the num of unique subjects by the num of BT subj treatments to get the N for each ind. cell
end

%Now Calc Z (SSP) matrix using recursive loops
Z=repmat(0,WIn.nCombs,WIn.nCombs);
Z=BTRLoop(1, repmat(1,length(Y),1),Z, N,Y,BT,WIn,S );

%Calc M matrices using other program
[M EpsList]=GenOrthogComps(WIn.nTreats);

%combine M and Z to calc epsilons
EpsHF=repmat(0,1,length(M));
EpsGG=EpsHF;

%calc totN, used for HF calc below
totN=N*prod(BT.nTreats);
for im=1:length(M)      
    S=M{im}*Z*M{im}';
    r=length(S);
    
    %calc EpsGG
    NumSS=0;
    for is1=1:r
        NumSS=NumSS+S(is1,is1);
    end
    DenSS=0;
    for is1=1:r
        for is2=1:r
            DenSS=DenSS+S(is1,is2)^2;
        end
    end           
    EpsGG(im)=(NumSS^2)/(r*DenSS);
    
    %Calc EpsHF
    EpsHF(im)=( (totN-g+1)*r*EpsGG(im)-2)/( r*(totN-g-r*EpsGG(im)) ); 
    %EpsHF(im)=( N*r*EpsGG(im)-2)/( r*(N-g-r*EpsGG(im)) ); %I have observed that to match SAS and SPSS, you would want to use this line of code instead of the line above it  
                                                           %BUT... According to what I've read, SPSS and SAS are wrong in this respect...
                                                           %see: http://archives.devshed.com/forums/development-94/huynh-feldt-r-vs-sas-bug-293023.html  
                                                           % and: Chen, RS and Dunlap, WP "A MonteCarlo STudy on teh Performance of a Corrected Formula for eps suggested by Lecoutre",
                                                           % Journal of Educational Statistics (1994) 
end

%these vals can't be more than 1
EpsGG=min(EpsGG,1);
EpsHF=min(EpsHF,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Z=BTRLoop(curFac,curSelInds,Z, N,Y,BT,WIn,S )

for iT=1:BT.nTreats(curFac)
    if BT.nTreats(curFac)==1 && BT.nFacs==1
        nextSelInds=curSelInds;
    else
        nextSelInds=curSelInds & BT.Facs(:,curFac)==BT.GNums{curFac}(iT);
    end
    if curFac==BT.nFacs                     %for DevMat         for FacLevList
        DevMat= WInRLoop(1, nextSelInds, repmat(0,WIn.nCombs,N),repmat(0,1,WIn.nFacs), N,Y,WIn,S );
        Z=Z+DevMat*DevMat';
    else
        Z=BTRLoop(curFac+1, nextSelInds,Z ,N,Y,BT,WIn,S  );
    end
end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DevMat=WInRLoop(curFac,curSelInds,DevMat,FacLevList, N,Y,WIn,S )

%if curFac==1;       DevMat=repmat(0,WIn.nCombs,N);      end

for iT=1:WIn.nTreats(curFac)
    FacLevList(curFac)=iT;
    if curFac==WIn.nFacs
        finSelInds = curSelInds & WIn.Facs(:,curFac)==WIn.GNums{curFac}(iT);
        curComb=sum( (FacLevList-1).*WIn.FacMult )+1;
                
        if isempty(S)
            DevMat(curComb,:)=Y(finSelInds)-mean(Y(finSelInds));    %works if subjects are all input in the same order...
        else
            [junk finSelInds2]=sort(S(finSelInds));     %still assuming one val for each subj, but not assuming that they've been input in teh proper order... NOTE- this could be changed later if it'es ever needed to account for more than one val per subject...
            tmpInds=find(finSelInds);
            DevMat(curComb,:)=Y(tmpInds(finSelInds2))-mean(Y(finSelInds));
        end
    else        
        DevMat=WInRLoop(curFac+1, curSelInds & WIn.Facs(:,curFac)==WIn.GNums{curFac}(iT) , DevMat,FacLevList, N,Y,WIn,S );
    end
end


