function [p,t,stats,terms] = anovan(y,group,model,sstype,gnames,displayopt)
%ANOVAN N-way analysis of variance (ANOVA).
%   ANOVAN performs multi-way (N-way) analysis of variance for comparing
%   means that are grouped in several ways.
%
%   P=ANOVAN(Y,GROUP,MODEL,SSTYPE,GNAME,DISPLAYOPT) performs anova on the
%   vector Y grouped by entries in the cell array GROUP.
%   GROUP must be a cell array, with each cell containing a grouping
%     variable.  The grouping variable can be a numeric vector, a
%     character matrix, or a single-column cell array of strings.
%     Each grouping variable must have the same number of items as Y. 
%   MODEL is an indication of the model to be used:
%     'linear' to use only main effects of all factors (default)
%     'interaction' for main effects plus two-factor interactions
%     'full' to include interactions of all levels
%     an integer representing the maximum interaction order, for example
%        3 means main effects plus two- and three-factor interactions
%     a vector V of integers, with each element describing a term, so
%        the Ith term includes the Jth grouping variable if
%        BITGET(V(I),J)==1
%   SSTYPE is 1, 2, or 3 for the type of sum of squares to use (default 3).
%   GNAMES is a character matrix of names for the grouping variables
%     (default has rows 'X1','X2',...).
%   DISPLAYOPT is 'on' (default) to display table, 'off' to omit display
%   P is a returned vector of p-values, one for each term.
%
%   [P,T,STATS,TERMS]=ANOVAN(...) also returns
%     T, the anova table as a cell array,
%     STATS, a structure of statistics used by the MULTCOMPARE function,
%     TERMS, a vector of integers describing terms used (suitable
%        for use as the MODEL input argument if you run ANOVAN again).
%
%   Example:
%	  [p, atab] = anovan(MPG, {Cylinders Origin Year}, 2, 2, ...
%                       strvcat('Cylinders', 'Origin', 'Year'))
%    This example performs three-way anova on MPG using the factors
%    Cylinders, Origin, and Year.  The model will have all two-factor
%    interactions but not the three-factor interaction.  Sum of
%    squares are Type 2.
%
%   See also ANOVA1, ANOVA2, MANOVA1.

%   References:
%      Dunn, O.J., and V.A. Clark (1974), "Applied Statistics:
%         Analysis of Variance and Regression," Wiley, New York.
%      Milliken, G.A., and D.E. Johnson, "Analysis of Messy Data,"
%         Chapman & Hall, New York.
%      Seber, G.A.F. (1977), "Linear Regression Analysis," Wiley,
%         New York.

%   Author:  Tom Lane, 5-19-99
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 1.8 $  $Date: 2000/05/26 17:28:09 $

% First:  argument checking and selection of defaults
error(nargchk(2,6,nargin));
if (prod(size(y))~=length(y)),error('The first argument must be a vector.'),end
if (size(y,1)==1), y=y(:); end
n = length(y);

if (~iscell(group)), error('The second argument must be a cell array.'), end
group = group(:);
ng = length(group);
if (ng>53), error('No more than 53 groups allowed.'), end

if (nargin < 5 | isempty(gnames))
   gnames = cellstr([repmat('X',ng,1) strjust(num2str((1:ng)'),'left')]);
else
   if (iscell(gnames)), gnames = gnames(:); end
   if (size(gnames,1) ~= ng)
      error(sprintf('Fifth argument must have %d rows', ng));
   elseif (~(ischar(gnames) | iscellstr(gnames)))
      error('Fifth argument must be a character matrix or cell array of strings.');
   elseif (ischar(gnames))
      gnames = cellstr(gnames);
   end
end

if (nargin < 4 | isempty(sstype))
   sstype = 3;
elseif (sstype ~= 1 & sstype ~= 2 & sstype ~=  3)
   error('Fourth argument must be 1, 2, or 3.');
end

if (nargin < 3 | isempty(model))
   model = 1;
elseif (~isnumeric(model))
   if strcmp(model, 'linear') | strcmp(model, 'l')
      model = 1;
   elseif strcmp(model, 'interaction') | strcmp(model, 'i')
      model = 2;
   elseif strcmp(model, 'full') | strcmp(model, 'f')
      model = ng;
   else
      error('Third argument is invalid.');
   end
end
if (any(model~=floor(model)|model<=0)), error('Third argument is invalid.'),end
if (length(model)==1)
   termlist = makemodel(model, ng);
else
   termlist = model(:);
end
if (nargin<6), displayopt = 'on'; end

% STEP 1:  Remove NaNs and prepare grouping variables.
% First pass, make sure all groups are rows and find NaNs.
nanrow = isnan(y);
for j=1:ng
   gj = group{j};
   if (size(gj,1) == 1), gj = gj(:); end
   if (size(gj,1) ~= n)
      error(sprintf('Group variable %d must have %d elements.',j,n));
   end
   if (ischar(gj)), gj = cellstr(gj); end
   group{j} = gj;
   if (isnumeric(gj))
      nanrow = (nanrow | isnan(gj));
   else
      nanrow = (nanrow | strcmp(gj,''));
   end
end

% Second pass, remove NaNs and create group index and name arrays.
y(nanrow) = [];
n = length(y);
dcode = 1;
gdum = cell(ng,1);
dfvar = zeros(ng,1);
if (nargout>2)
   allgrps = zeros(n, ng);
   allnames = cell(ng,1);
end
for j=1:ng
   gj = group{j};
   gj(nanrow,:) = [];
   [gij,gnj] = grp2idx(gj);
   dfvar(j) = size(gnj,1) - 1;
   if (nargout>2)
      allgrps(:,j) = gij;
      allnames{j} = gnj;
   end

% STEP 2:  Create dummy variable arrays for each grouping variable.
   gdum{j} = idummy(gij, dcode);
end

% STEP 3:  Create dummy variable arrays for each term in the model.
nterms = length(termlist);
[sterms,sindex] = sort(termlist);
ncols = 1;

termdum = cell(nterms, 1);
tnames = cell(nterms, 1);
dfterm0 = zeros(nterms, 1);   % nominal d.f. for each term

% For each term,
for j=1:nterms
   % Loop over elements of the term
   df0 = 1;
   tm = sterms(j);
   tdum = [];         % empty term so far
   tn = '';           % blank name so far
   while(tm)
      % Get last-numbered effect remaining
      lne = 1 + floor(log2(tm));
      tm = bitset(tm, lne, 0);
      df0 = df0 * dfvar(lne);

      % Combine its dummy variable with the part computed so far
      tdum = termcross(tdum, gdum{lne});

      % Get the name of this variable
      if (isempty(tn))
         tn = gnames{lne};
      else
         tn = [tn '*' gnames{lne}];
      end
      
      % If the rest of this term is computed, take advantage of that
      oldtm = find(tm == sterms(1:j-1));
      if (length(oldtm) > 0)
         k = oldtm(1);
         tdum = termcross(tdum, termdum{sindex(k)});
         tn = [tnames{sindex(k)} '*' tn];
         df0 = df0 * dfterm0(sindex(k));
         break;
      end
   end

   % Store this term's dummy variables and name
   k = size(tdum, 2);
   termlength(sindex(j),1) = k;
   ncols = ncols + k;
   sj = sindex(j);
   termdum{sj} = tdum;
   if (isempty(tn)), tn = 'Constant'; end
   tnames{sj,1} = tn;
   dfterm0(sj) = df0;
end
tnames{length(tnames)+1,1} = 'Error';

% STEP 4:  Create the full design matrix
dmat = ones(n, ncols);
termname = zeros(ncols,1);
termcode = zeros(ncols,1);
termstart = cumsum([2; termlength(1:end-1)]);
termend = termstart + termlength - 1;
for j=1:nterms
   clist = termstart(j):termend(j);
   dmat(:, clist) = termdum{j};
   termname(clist) = j;
   termcode(clist) = termlist(j);
end

% STEP 5:  Fit the full model and compute the residual sum of squares
mu =  mean(y);
y = y - mu;                   % two passes to improve accuracy
mu2 = mean(y);
mu = mu + mu2;
if (nargout>2) y0 = y; end
sst = sum((y - mu2).^2);
[ssx,dfx,dmat,y,termname,stats] = initfit(dmat, y, termname, sst, mu);
sse = sst - ssx;
dfe = n - dfx;

% STEP 6:  Determine which models to compare for testing each term
ssw  = -ones(nterms, 1);      % sum of squares with this term
sswo = ssw;                   % sum of squares without this term
dfw  = ssw;                   % residual d.f. with this term
dfwo = ssw;                   % residual d.f. without this term

switch sstype
 case 1
   modw = tril(ones(nterms)); % get model with this term
   k = nterms;                % locations of model with all terms
 case 2
   modw = termsnotcontained(termlist);
   k = (sum(modw,2) == nterms);
 case 3
   modw = ones(nterms);
   k = 1:nterms;
end
modw = logical(modw);
modwo = logical(modw - eye(nterms));   % get model without this term

% STEP 7:  Fit each model, get its residual SS and d.f.
ssw(k) = ssx;        % for full model we already know the results
dfw(k) = dfx;

% For now fit each model separately; may add updating later
dfboth = [dfw; dfwo];
allterms = 1:nterms;
while 1
   % Find the next model index to fit
   k = find(dfboth == -1);
   if (isempty(k)), break, end
   k = k(end);
   
   % Find the model with this index
   if (k > nterms)
      thismod = modwo(k-nterms, :);
   else
      thismod = modw(k, :);
   end
   
   % Get the design matrix for this model
   keepterms = allterms(thismod);
   clist = ismember(termname, [0 keepterms]);
   X = dmat(:,clist);

   % Fit this model
   [ssx0,dfx0] = anovanfit(X, y);
   
   % Use these results for each term that requires them
   mod0 = repmat(thismod, nterms, 1);
   k = find(all(modw == mod0, 2));
   ssw(k) = ssx0;
   dfw(k) = dfx0;
   dfboth(k) = 0;
   k = find(all(modwo == mod0, 2));
   sswo(k) = ssx0;
   dfwo(k) = dfx0;
   dfboth(nterms+k) = 0;
end
clear mod0

% STEP 8:  Compute the sum of squares attributed to each term
ssterm = ssw - sswo;
dfterm = dfw - dfwo;
ssterm(dfterm==0) = 0;    % make this exact

% STEP 9:  Compute the mean square for each term
msterm = ssterm ./ max(1, dfterm);
mse = sse * (dfe>0) / max(1, dfe);

% STEP 10:  Compute an F statistic for each term
if (mse > 0)
   fstat = msterm ./ mse;
else
   fstat = repmat(Inf, size(msterm));
end

% STEP 11:  Compute the p-value for each term
pval = repmat(NaN, size(fstat));
if (dfe > 0)
   ok = (dfterm > 0);
   if (any(ok)), pval(ok) = 1 - fcdf(fstat(ok), dfterm(ok), dfe); end
end

p = pval;

% Create anova table as a cell array
sing = (dfterm < dfterm0);
tbl = [ssterm dfterm sing msterm fstat pval; sse dfe 0 mse NaN NaN];
tbl = num2cell(tbl);
tbl(end,end-1:end) = {[]};
tbl = [{'Source' 'Sum Sq.' 'd.f.' 'Singular?' 'Mean Sq.' 'F' 'Prob>F'}; ...
       tnames tbl];

% Add total information
tbl{end+1,1} = 'Total';
tbl{end,2} = sst;
tbl{end,3} = n-1;
tbl{end,4} = 0;

if (nargout > 1), t = tbl; end
if (nargout > 3), terms = termlist; end

% For display, mark terms that are rank deficient
star = ' *';
tbl{1,1} = ['  ' tbl{1,1}];
for r=2:size(tbl,1)
   x = tbl{r,1};
   if (tbl{r,4})
      x = ['# ' x];
   else
      x = ['  ' x];
   end
   tbl{r,1} = x;
end
tbl = tbl(:,[1 2 3 5 6 7]);

if (isequal(displayopt,'on'))
   switch(sstype)
    case 1,    cap = 'Sequential (Type I) sums of squares.';
    case 2,    cap = 'Hierarchical (Type II) sums of squares.';
    otherwise, cap = 'Constrained (Type III) sums of squares.';
   end
   if (any(sing)), cap = [cap '  Terms marked with # are not full rank.']; end
   digits = [-1 -1 -1 -1 2 4];
   statdisptable(tbl, 'N-Way ANOVA', 'Analysis of Variance', cap, digits);
end

% Store additional information for multiple comparisons
stats.terms = termlist;
stats.nlevels = dfvar+1;
stats.termcols = [1; termlength];

% --------------------------
function m = termsnotcontained(t)
%TERMSNOTCONTAINED Creates a logical matrix indicating term containment
%   m(i,j)==1  iff  t(i) is not contained by t(j)
%   set diagonals to 1 because we want proper containment

n = length(t);
a = t(:,ones(n,1));
m = (bitand(a',a) ~= a);
m(1:n+1:n^2) = 1;

% --------------------------
function v = makemodel(p,ng)
%MAKEMODEL Helper function to make model vector from scalar
%   P = max model term size, NG = number of grouping variables
vgen = 2.^(0:ng-1)';
v = vgen;
for j=2:p
   v = [v; sum(nchoosek(vgen,j),2)];
end

% --------------------------
function [ssx,dfx] = anovanfit(dmat, y)
%ANOVANFIT Fit one model during an n-way anova calculation

% Tolerance for computing rank from diag(R) after QR decomposition
TOL = 100 * eps * size(dmat,2);      

[q,r,e] = qr(dmat,0);
if (size(r,2) > 1)
   d = abs(diag(r));
else
   d = abs(r(1,1));
end
p = sum(d > TOL);

y = y' * q(:,1:p);
ssx = sum(y(1:p).^2);
dfx = p;

% --------------------------
function [ssx,dfx,dmat2,y2,termname2,stats] = initfit(dmat,y,termname,sst,mu)
%INITFIT Do initial fit and reduce data for subsequent fits

% Tolerance for computing rank from diag(R) after QR decomposition
ncols = size(dmat,2);
TOL = 100 * eps * ncols;

% Do qr decomposition on current design matrix
[q,r,e] = qr(dmat,0);
if (size(r,2) > 1)
   d = abs(diag(r));
else
   d = abs(r(1,1));
end
p = sum(d > TOL);

% Get information for fit of y on this matrix
yy = (y' * q(:,1:p))';
ssx = sum(yy.^2);
dfx = p;

% Return reduced design matrix.
n = length(y);
m = min(n, ncols);
dmat2 = r(1:p,:);

% Rotate y into the column space of dmat and keep only that portion.
y2 = yy(1:min(n,p));

% Re-arrange term names to correspond to new design matrix
termname2 = termname(e);

% Prepare for multiple comparisons if requested
if (nargout > 5)
   stats.source = 'anovan';
   coeffs = dmat2\y2;
   t = zeros(ncols,1);
   t(:) = NaN;
   t(e(1:length(coeffs))) = coeffs;
   t(1) = t(1) + mu;
   stats.coeffs = t;
   RR(:,e) = dmat2;
   stats.Rtr = RR';
   stats.dfe = length(y) - dfx;
   stats.mse = (sst - ssx) / max(1, stats.dfe);
   stats.permvec = e;
end
