function anovalc(LMatrix, y, group, stats, planned)
% format anovaLC(LMatrix, y, group, stats,planned)
% Input LMatrix
% Each row of the LMatrix should contain a linear contrast
%   LMatrix = [-1 1 0 0;-0.5 0.5 0 0] will return identical contrasts
%   y=data in a column vector
%   group is the column vector indicating group membership
%   stats is output from anova1, anova2 or anovan
%   planned =1 if the contrast was planned a priori
%   planned =0 if the contrast was not planned, in which case Scheffe
%              multipliers will be used.
% Written by Eugene D. Gallagher 12/7/2010
if nargin<5;planned=1;end
[R,C]=size(LMatrix);
% Create placeholder vectors for the output of the data
G=unique(group);   % Contains indices indicating treatment membership
n=zeros(1,C);
meanDATA=zeros(1,C);
sumDATA=zeros(1,C);
SSC=zeros(R,1);
F=zeros(R,1);
Fprob=zeros(R,1);
g=zeros(R,1);
seg=zeros(R,1);
tdf=tinv(0.975,stats.dfe);
for j=1:C
    i=find(group==G(j));
    n(j)=length(i);
    sumDATA(j)=sum(y(i));
    meanDATA(j)=mean(y(i));
end
for i=1:R % do each linear contrast
        sumLM=sum(LMatrix(i,:));
        sumabsLM=sum(abs(LMatrix(i,:)));
        fprintf('\nContrast Result Number %1.0f:\n',i)
        format rat
         disp(LMatrix(i,:));
         format
        if abs(sumLM)>=3*eps
            error('Linear contrasts must sum to 0');
        elseif abs((sumabsLM-2))>eps
            % This corrects an issue that is found in PASW, in which
            % ANOVA doesn't allow fractional linear contrasts and the
            % effects size and standard error are wrong if a contrast
            % such as [-1 -1 2 0] is used, in which case the sum of
            % the absolute value of the contrasts is 4, not 2 and the
            % estimated effect size and standard are 2x too large.
            LMatrix(i,:)=1/(sumabsLM/2)*LMatrix(i,:);
            fprintf(...
                'Linear Contrast %1.0f converted to equivalent form:\n',i)
            format rat
                disp(LMatrix(i,:))
            format
        end
        SSC(i)=sum(LMatrix(i,:).*sumDATA./n)^2/sum(LMatrix(i,:).^2./n); 
        % Calculate the value of the linear contrast g (from Sleuth)
        g(i)=sum(LMatrix(i,:).*meanDATA);
        % The equation for the standard error of the linear contrast
        % can be found in Statistical Sleuth Chapter 6
        seg(i)=sqrt(stats.mse).*sqrt(sum(LMatrix(i,:).^2./n));
        F(i)=SSC(i)/stats.mse;
        Fprob(i)=1-fcdf(F(i),1,stats.dfe);
        if planned==1
            fprintf('The difference in means is %5.2f +/- %5.2f\n',...
                g(i),seg(i)*tdf)
        else
            Scheffe=sqrt((C-1)*finv(1-0.05,C-1,stats.dfe));
            fprintf(...
      'The difference in means is %5.2f +/- %5.2f (Scheffe Interval)\n',...
                g(i),seg(i)*Scheffe)
        end
        fprintf('\n Source      SS     df    MS    F   Prob\n')
        fprintf(...
          'Contrast    %4.1f     1   %4.1f %4.1f %5.3g\n',SSC(i),SSC(i),...
             F(i),Fprob(i))
  fprintf('   Error   %4.1f    %2.0f %5.2g\n',stats.mse*stats.dfe,...
      stats.dfe,stats.mse)
end