function [s, cfg] = pd_statfun_indepsamplesANCOVA_2var(cfg, dat, design)

% PD_STATFUN_INDEPSAMPLESANCOVA calculates the independent samples F-statistic on the
% biological data in dat (the dependent variable), using the information on the
% independent variable (ivar) and the continuous covariate (cvar) in design.

% $Id$

% set the defaults
if ~isfield(cfg, 'computestat'),       cfg.computestat='yes';     end
if ~isfield(cfg, 'computecritval'),    cfg.computecritval='no';   end
if ~isfield(cfg, 'computeprob'),       cfg.computeprob='no';      end
if ~isfield(cfg, 'alpha'),             cfg.alpha=0.05;            end
if ~isfield(cfg, 'tail'),              cfg.tail=1;                end


% perform some checks on the configuration
if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no')
  ft_error('P-values can only be calculated if the test statistics are calculated.');
end
if ~isfield(cfg,'cvar')
  ft_error('Please specify two continuous covariates for ANCOVA using cfg.cvar');
end


ncond=length(unique(design(cfg.ivar,:)));
nrepl=0;
for condindx=1:ncond
    nrepl=nrepl+length(find(design(cfg.ivar,:)==condindx));
end
if nrepl<size(design,2)
  ft_error('Invalid specification of the independent variable in the design array.');
end
if nrepl<=ncond
    ft_error('The must be more trials/subjects than levels of the independent variable.');
end
dfnum = ncond - 1;
dfdenom = nrepl - ncond - 2;
dimdat = size(dat);
dimdat = dimdat(1);

s.stat = NaN(dimdat,1);
if strcmp(cfg.computestat, 'yes')
    for ipt = 1:dimdat % looping over the data input
        cvar1 = design(cfg.cvar(1),:); % first covariate
        cvar2 = design(cfg.cvar(2),:); % second covariate
        ivar = design(cfg.ivar,:);
        X = [cvar1(:),cvar2(:),ivar(:)]; % set up variables
        mdl = fitlm(X,dat(ipt,:)'); % use MATLAB functions to run ANCOVA
        anmdl = anova(mdl);

        s.stat(ipt) = anmdl.F(3); % extract F-statistic for effect of ivar, covarying for cvar1 and cvar2
    end
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  s.dfnum   = dfnum;
  s.dfdenom = dfdenom;
  if cfg.tail==-1
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a left tail critical value.');
  end
  if cfg.tail==0
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a two-sided critical value.');
  end
  if cfg.tail==1
    s.critval = finv(1-cfg.alpha,s.dfnum,s.dfdenom);
  end
end

if strcmp(cfg.computeprob,'yes')
  % also compute the p-values
  s.dfnum   = dfnum;
  s.dfdenom = dfdenom;
  if cfg.tail==-1
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a left tail p-value.');
  end
  if cfg.tail==0
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a two-sided p-value.');
  end
  if cfg.tail==1
    s.prob = 1-fcdf(s.stat,s.dfnum,s.dfdenom);
  end
end

end