% Description:
%
%     Remove one term at a time from a MANCOVAN model until all p-values are
%     below the specified threshold.
%
% Syntax:
%
%     [ T, p, stats ] = mStepwise(Y, groups, [], threshold)
%     [ T, p, stats ] = mStepwise(Y, [], covariates, threshold)
%     [ T, p, stats ] = mStepwise(Y, groups, covariates, threshold)
%     [ T, p, stats ] = mStepwise(Y, groups, covariates, threshold, options)
%
%     [ T, p, stats ] = mStepwise(Y, X, terms, threshold, options)
%
% Inputs:
%
%     Y          - [ N x M ] (double) - multivariate response
%     groups     - [ N x G ] (int)    - qualitative variables
%     covariates - [ N x C ] (double) - quantitative variables
%     threshold  - [ 1 x 1 ] (double) - p-value at which to stop dropping terms
%     options    - [ 1 x P ] (cell)   - see Options
%
%     X     - [ N x T ] (double) - design matrix
%     terms - [ 1 x T ] (cell)   - model terms
%
% Outputs:
%
%     T - [ (G + C) x 1 ] (double)
%     p - [ (G + C) x 1 ] (double)
%
%     stats.U - [ N x N ] (double) - U resulting from the SVD (option 'SVD').
%     stats.S - [ N x N ] (double) - S resulting from the SVD (option 'SVD').
%     stats.V - [ M x N ] (double) - V resulting from the SVD (option 'SVD').
%
%     stats.BIC - [ 1 x P ] (double) - values of the Bayesian information
%         criterion (BIC) associated with the SVD (option 'SVD').
%
%     stats.PVE - [ N x 1 ] (double) - the percentage of variance in Y explained
%         by each consecutive column of U (option 'SVD').
%
%     stats.Terms - [ 1 x B ] (cell) - full model terms numbering groups from
%         one through size(groups, 2) and covariates from size(groups, 2) + 1
%         through size(groups, 2) + size(covariates, 2) + 1, and interactions
%         according to combinations of these numbers.
%
%     stats.X - [ N x B ] (double) - full model design matrix, including columns
%         for the grand mean, groups, covariates, and interactions.
%
%     stats.Y - [ N x M ] (double) - the multivariate response used for the
%         computation.
%
%     stats.B - [ B x P ] (double) - regression coefficients associated with the
%         full model after stepwise removal of insignificant terms.
%
%     stats.SSE - [ P x P ] (int) - full model sum of squared errors associated
%         with each column of Y.
%
%     stats.DFE - [ 1 x M ] (int) - degrees of freedom used in the computation
%         of stats.MSE.
%
%     stats.MSE - [ P x P ] (int) - full model mean squared errors associated
%         with each column of Y.
%
% Details:
%
% Options:
%
%     'group-group'         - include group-group interactions
%     'covariate-covariate' - include covariate-covariate interactions
%     'group-covariate'     - include group-covariate interactions
%     'over-determined'     - use over-determined coding for the design matrix
%     'sigma-restricted'    - use sigma-restricted coding for the design matrix
%     'SVD'                 - reduce the dimensionality of Y using an SVD
%     'verbose'             - display extra information to the command window
%
% Examples:
%
%     The following example uses a simple additive model with covariates, but no
%     interactions and avoids using the Statistics Toolbox:
%
%         n          = 100; 
%         groups     = round(2 * rand(n, 2) + 0.5);
%         covariates = 10 * randn(n, 2);
%         Y          = groups + covariates + randn(n, 2);
%
%         [ T, p, stats ] = mStepwise(Y, groups, covariates, 0.10, ...
%             { 'group-group' 'covariate-covariate' 'group-covariate' 'verbose' });
%
%     For other examples, refer to mancovan.m.  The setup and syntax is exactly
%     the same in every case except for the addition of a threshold.
%
% Notes:
%
% Author(s):
%
%     William Gruner (williamgruner@gmail.com)
%
% References:
%
%     Refer to the references listed in mancovan.m.
%
% Acknowledgements:
%
%     Many thanks to Dr. Erik Erhardt and Dr. Elena Allen of the Mind Research
%     Network (www.mrn.org) for their continued collaboration.
%
% Version:
%
%     $Author: williamgruner $
%     $Date: 2010-04-15 07:46:14 -0600 (Thu, 15 Apr 2010) $
%     $Revision: 496 $

function [ T, p, stats ] = mStepwise(Y, groups, covariates, threshold, options)
    
    if nargin == 0
        T = BIT(); return
    end
    
    if ~exist('options', 'var')
        options = cell(0);
    end
    
    if iscell(covariates)
        X     = groups;
        terms = covariates;
    else
        [ X, terms ] = mX(groups, covariates, options);
    end
        
    if ~isempty(strmatch('SVD', options, 'exact'))
        
        [ BIC, U, S, V ] = mSVD(Y, options);

        stats.U   = U;
        stats.S   = S;
        stats.V   = V;
        stats.BIC = BIC;
        stats.PVE = cumsum(diag(S) ./ sum(diag(S)));
        
        b = find(diff(BIC) > 0);
        b = b(1);
        
        U = U(  :, 1:b);
        S = S(1:b, 1:b);
        V = V(  :, 1:b);
        
    else
        
        U = Y;
        
    end
    
    if ~isempty(strmatch('verbose', options, 'exact'))
        fprintf('\n')
    end
            
    while length(terms) > 1
    
        if ~isempty(strmatch('over-determined', options, 'exact'))
            M = X * pinv(X' * X) * X';
        else
            M = X * inv(X' * X) * X';
        end
        
        B = mUnique(terms(2:end));
        T = repmat(NaN, length(B), 1);
        p = repmat(NaN, length(B), 1);

        for i = 1 : length(B)

            I0 = mTerms(B{i}, terms);
            X0 = X(:, I0);
            
            if ~isempty(strmatch('over-determined', options, 'exact'))
                M0 = X0 * pinv(X0' * X0) * X0';
            else
                M0 = X0 * inv(X0' * X0) * X0';
            end
            
            if ~isempty(strmatch('verbose', options, 'exact'))
                mDispModels([ {0} B ], mUnique(terms(I0)))
            end

            [ T(i), p(i) ] = mLHT(U, X, X0, M, M0, options);

            if ~isempty(strmatch('verbose', options, 'exact'))
                fprintf('\n')
            end

        end

        if max(p) < threshold
            break
        end
        
        [ q, I ] = sort(p, 'descend');
        
        for i = 1 : length(I)
            if mIsInteractionTerm(B{I(i)})
                break
            elseif mIsMainTerm(B{I(i)})
                if ~ismember(B{I(i)}, cat(2, B{mFindInteractionTerms(B)}))
                    break
                end
            end
        end
        
        if q(i) < threshold
            
            break
            
        else
            
            if ~isempty(strmatch('verbose', options, 'exact'))
                
                if mIsInteractionTerm(B{I(i)})
                    fprintf('Removing %d%d ...\n', B{I(i)}(1), B{I(i)}(2))
                elseif mIsMainTerm(B{I(i)})
                    fprintf('Removing %d ...\n', B{I(i)})
                end
                
                fprintf('\n')
                
            end
        
            X    (:, mFindTerms(B{I(i)}, terms)) = [];
            terms(:, mFindTerms(B{I(i)}, terms)) = [];
            
        end
        
    end
    
    if ~isempty(strmatch('over-determined', options, 'exact'))
        stats.B = pinv(X' * X) * X' * U;
    else
        stats.B = inv(X' * X) * X' * U;
    end
    
    stats.Terms     = terms;
    stats.X         = X;
    stats.Y         = Y;
    stats.SSE       = U' * (eye(size(M)) - M) * U;
    stats.DFE       = size(X, 1) - rank(X);
    stats.MSE       = stats.SSE / stats.DFE;
    stats.Residuals = U - stats.X * stats.B;

function b = BIT()

    b = true;
    
    % Compare the results to those obtained from the previous version.

    s = load('mStepwise-BIT-1.mat');

    [ T, p, stats ] = mStepwise(s.Y, s.groups, s.covariates, 0.10, ...
        { 'sigma-restricted' 'group-group' 'covariate-covariate' 'group-covariate' 'verbose' });

    e = s.T - T;
    
    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end
    
    e = s.p - p;
    
    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end

    % Compare the results to those obtained from the previous version.
    
    s = load('mStepwise-BIT-2.mat');

    [ T, p, stats ] = mStepwise(s.Y, s.groups, s.covariates, 0.01, ...
        { 'sigma-restricted' 'group-group' 'covariate-covariate' 'group-covariate' 'verbose' });
    
    e = s.T - T;
    
    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end
    
    e = s.p - p;

    if any(abs(e(:)) > sqrt(eps))
        b = false;
    end
    