function  [SS,df] = sumsq(DVs,IVs,sets)

% SUMSQ  -- sums of squares and degrees of freedom
%
% [SS,df] = sumsq(DVs,IVs,sets)
%
% DVs and IVs are column-oriented matrices of dependent and independent
% variables, respectively. SETS is a cell array of length K. Each set 
% is a vector of "split indices" in IVs. For example, {[],[1],[2],[1,2]} 
% specifies the two-factor design {total,A,B,AxB}. When there is only
% one split set, SETS=[...] can be used as a shorthand for SETS={[...]}.
% Use PARTITION to generate full-factorial split sets.
%
% SS and DF are matrices with K rows and size(DVS,2) columns. Each split
% set is processed separately in the following way for each DV:
%  1. s = split(DV,IVs(:,sets{k}))    % split into treatment "cells"
%  2. n = aggreg(s,'length')          % sample sizes, possibly unequal
%  3. T = aggreg(s,'sum')             % treatment totals
%  4. SS(k) = sum((T.^2)./n) - CF     % sum of squares around cell means
%  5. df(k) = length(s)-1             % simple-minded degrees of freedom
% The correction factor in line 4 is CF=(sum(DV))^2/N.
% Be sure to correct the DF for interactions: if df(A)=1 and df(B)=2,
% for instance, then df(AB)=1 while this function returns df(AB)=5.
% Both SS and df for an interaction are SScells (see p. 386 and 468 in
% Howell, D. C. (1992). Statistical methods for psychology). The true
% term is calculated by:  SS(AB)=SS(cellsAB)-SS(A)-SS(B); df(AB)=...
%
% Example: One-way ANOVA:
%  [SS,df]=sumsq(DV,IV,{[],[1]});
%  SStotal=SS(1);SStreatment=SS(2);SSerror=SStotal-SStreatment;
%  MStreatment=SStreatment/df(2);dferr=df(1)-df(2);MSerr=SSerr/dferr;
%  F=MStreatment/MSerror; p=1-fcdf(F,df(2),dferr);
%
% See also ANOVA, PARTITION, SPLIT, AGGREG, VAR.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2002/12/20 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

%---  Handle input arguments
[N,N_DVs] = size(DVs) ;
[n,N_IVs] = size(IVs) ;
if (N~=n)
  error('The DVs and IVs matrices must have the same number of rows.') ;
end

if (isnumeric(sets))  sets = {sets} ;  end    % shorthand
K = length(sets) ;

%--- Correction factor
CF = sum(DVs) .^2 / N ;       % [1xN_DVs]

%--- Main loop
SS = zeros(K,N_DVs) ;
df = zeros(K,N_DVs) ;

for k=1:K
  iv = IVs(:,sets{k}) ;
  if  (isempty(iv))   % special case: total SS
    T = DVs ;         % each replication is a "cell"
    SS(k,:) = sum(T.^2) - CF ;
    df(k,:) = N - 1 ;
  else 
    s = split(DVs,iv) ;         % split into treatment "cells", [g*1]
    g = length(s) ;             % number of cells
    n = aggreg(s,inline('size(x,1)')) ;     % sample sizes, [g*N_DVs]
    T = aggreg(s,inline('sum(x,1)')) ;      % treatment totals, [g*N_DVs]
    SS(k,:) = sum(T.^2 ./ n) - CF ;         % sum of squares, [1*N_DVs]
    df(k,:) = g-1 ;             % simple-minded degrees of freedom
  end
end

%--- Return SS and DF
%%%%% End of SUMSQ.M
