myregrcomp(x)

function myregrcomp(x)
%MYREGCOMP: Compare two linear regression.
%This function compares two least-square linear regression.
%Test are implemented as reported by Stanton A. Glantz book
%This routine uses MYREGR function. If it is not present on
%the computer, myregrcomp will try to download it from FEX
%
% SEE also myregr, myregrcomp
%
% Syntax: 	x=myregrcomp(X)
%
%     Inputs:
%           X - data matrix (Size of matrix must be n-by-3;
%           x data=column 1, y data=column 2, regression=column 3). 
%     Outputs:
%           - Summary of MYREGR function
%           - Global test
%           - Test on slopes (eventually)
%           - Test on intercepts (eventually)
%
% Example:
%       a=[775 800 810 800 850 860 925 900 925 945 945 950 975 1050 1075
%       800 810 850 875 850 870 910 925 925 940 945 960 1100 1100 990]./1000;
%       b=log10([21 20 13.5 8.5 10.5 10 12.8 9.0 6.5 11 10.5 9.5 5.5 6 3.8
%       10 5 9.5 2.5 4 5.8 9.8 8 6 4.3 8.5 9 8.5 4.5 2.3]);
%       g=[repmat(1,15,1); repmat(2,15,1)];
%       x=[a' b' g];
%
%   Calling on Matlab the function: 
%             x=myregrcomp(x)
%
%   Answer is:
%
% --------------------------------------------------------------------
%                     REGR1              REGR2              REGR.TOT
% --------------------------------------------------------------------
% Points               15                 15                 30
% --------------------------------------------------------------------
% Slope               -1.7656            -0.1712            -1.0822
% Slope S.E.          0.3684             0.6466             0.4385
% --------------------------------------------------------------------
% Intercept           2.5802             0.9297             1.8661
% Intercept S.E.      0.3352             0.5997             0.4028
% --------------------------------------------------------------------
% R.S.E.              0.1243             0.2175             0.2100
% --------------------------------------------------------------------
%  
% GLOBAL TEST
% F-value : 6.6773    numerator d.f = 2     denominator d.f. = 26
% Critical value at 95% significance level = 3.3690     p-value = 0.0046
% These regressions are different
%  
% TEST ON SLOPES
% t-value : 3.6289    degrees of freedom = 26
% Critical value at 95% significance level = 1.7056     p-value = 0.0006
% The slopes are different
%  
% TEST ON INTERCEPTS
% t-value : 3.7197    degrees of freedom = 26
% Critical value at 95% significance level = 1.7056     p-value = 0.0005
% The intercepts are different
%
%           Created by Giuseppe Cardillo
%           [e-mail]
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2007) MyRegressionCOMP: a simple routine to compare two LS
% regression. 
% http://www.mathworks.com/matlabcentral/fileexchange/15953

idx=x(:,3); %index of the groups
%first regression
try
    [m1,q1,stat1]=myregr(x(idx==1,1),x(idx==1,2),0);
catch ME
    disp(ME)
    disp('I am trying to download the myregr function by Giuseppe Cardillo from FEX')
    [F,Status]=urlwrite('http://www.mathworks.com/matlabcentral/fileexchange/15473-myregression?controller=file_infos&download=true','myregr.zip')
    if Status
        unzip(F)
        [m1,q1,stat1]=myregr(x(idx==1,1),x(idx==1,2),0);
    end
    clear F Status
end

%second regression
[m2,q2,stat2]=myregr(x(idx==2,1),x(idx==2,2),0);
%total regression
x=sortrows(x,1);
[m3,q3,stat3]=myregr(x(:,1),x(:,2),0);

%assemble all
%points
n=[stat1.n,stat2.n];
%slopes
m=[m1.value,m2.value];
%standard errors on slopes
mse=[m1.se,m2.se];
%intercepts
q=[q1.value,q2.value];
%standard errors on intercepts
qse=[q1.se,q2.se];
%regressions standard errors
rse=[stat1.rse,stat2.rse];
%display results
tr=repmat('-',1,60');
disp(tr)
disp('                    REGR1              REGR2              REGR.TOT')
disp(tr)
fprintf('Points               %d                 %d                 %d\n',n,stat3.n)
disp(tr)
fprintf('Slope               %0.4f %18.4f %18.4f\n',m,m3.value)
fprintf('Slope S.E.          %0.4f %18.4f %18.4f\n',mse,m3.se)
disp(tr)
fprintf('Intercept           %0.4f %18.4f %18.4f\n',q,q3.value)
fprintf('Intercept S.E.      %0.4f %18.4f %18.4f\n',qse,q3.se)
disp(tr)
fprintf('R.S.E.              %0.4f %18.4f %18.4f\n',rse,stat3.rse)
disp(tr)
disp(' ')
disp('GLOBAL TEST')
vd=(sum(n)-4); %denominator degrees of freedom
%combined variance of separated regressions
cs=(sum((n-2).*rse.^2))/vd;
%variation in variance of separated regressions
vs=((vd+2)*stat3.rse^2-vd*cs)/2;
%F test
F=vs/cs; cv=finv(0.95,2,vd); p=1-fcdf(F,2,vd);
%display results
fprintf('F-value : %0.4f    numerator d.f = 2     denominator d.f. = %d\n',F,vd);
fprintf('Critical value at 95%% significance level = %0.4f     p-value = %0.4f\n',cv,p)
if F<=cv
    disp('These regressions are equal')
else
    disp('These regressions are different')
    disp(' ')
    disp('TEST ON SLOPES')
    cs=sum((n-2).*mse.^2)/vd;
    cse=realsqrt(cs.*sum(1./((n-1).*mse.^2)));
    t=abs(diff(m))/cse; cv=tinv(0.95,vd); p=1-tcdf(t,vd);
    fprintf('t-value : %0.4f    degrees of freedom = %d\n',t,vd);
    fprintf('Critical value at 95%% significance level = %0.4f     p-value = %0.4f\n',cv,p)
    if t<=cv
        disp('The slopes are equal')
    else
        disp('The slopes are different')
    end
    disp(' ')
    disp('TEST ON INTERCEPTS')
    cs=sum((n-2).*qse.^2)/vd;
    cse=realsqrt(cs.*sum(1./((n-1).*qse.^2)));
    t=abs(diff(q))/cse; cv=tinv(0.95,vd); p=1-tcdf(t,vd);
    fprintf('t-value : %0.4f    degrees of freedom = %d\n',t,vd);
    fprintf('Critical value at 95%% significance level = %0.4f     p-value = %0.4f\n',cv,p)
    if t<=cv
        disp('The intercepts are equal')
    else
        disp('The intercepts are different')
    end
end