% RegrCompare: Randomized permutation comparison of two or more linear regression 
%              lines for homogeneity.
%
%     Usage: [b,iof_stats,slope_stats,intcpt_stats] = regrcompare(x,y,g,{iter})
%
%         x =    [N x 1] vector for independent variable.
%         y =    matching vector for dependent variable.
%         g =    matching group-identification vector for k groups.
%         iter = optional number of random permutation iterations.
%         ----------------------------------------------------------------------
%         b = [2 x k] matrix of regression coefficients for each of k groups.
%                The first row contains the intercepts and the second the slopes.
%         iof_stats = vector of improvement-of-fit statistics:
%                F       - F statistic;
%                pr_asym - asymptotic significant level;
%                pr_boot - bootstrapped significant level;
%                sse(r)  - sum-of-squared error for reduced (pooled-regr) model;
%                sse(f)  - sse for full (separate-regr) model;
%                df(r)   - degrees of freedom for reduced model;
%                df(f)   - degrees of freedom for full model.
%         slope_stats = vector of anova statistics for heterogeneity of slopes:
%                F       - F statistic;
%                pr_asym - asymptotic significance level;
%                pr_boot - optional randomized significance level;
%                df1     - numerator degrees of freedom;
%                df2     - denominator degrees of freedom.
%         intcpt_stats = vector of anova statistics for heterogeneity of intercepts:
%                F       - F statistic;
%                pr_asym - asymptotic significance level;
%                pr_boot - optional randomized significance level;
%                df1     - numerator degrees of freedom;
%                df2     - denominator degrees of freedom.
%               

% RE Strauss, 10/14/02

function [b,iof_stats,slope_stats,intcpt_stats] = regrcompare(x,y,g,iter)
  if (nargin < 4) iter = []; end;
  
  if (isempty(iter))
    iter = 0;
  end;
  
  [grp_id,n] = uniquef(g);                    % Groups
  ngrps = length(grp_id);
  N = length(x);                              % Sample size
  
  [b_pooled,stats_pooled] = linregr(x,y);     % Pooled regression
  [b_sep,stats_sep] = linregr([x,g],y);       % Separate regressions
  
  b = zeros(2,ngrps);
  s = zeros(2,ngrps);
  ssef = 0;
  for i = 1:ngrps                             % Regressions by group
    ig = find(g==grp_id(i));
    [b(:,i),stats,p,r,s(:,i)] = linregr(x(ig),y(ig));
    s(:,i) = s(:,i)*sqrt(n(i));
    ssef = ssef + stats(2)*stats(5);
  end;
  
  dfr = N-2;                                           % Improvement-of-fit statistics
  dff = N-2*ngrps;
  sser = stats_pooled(2)*stats_pooled(5);     
  F_iof = (((sser-ssef)-(dfr-dff))/(dfr-dff))/(ssef/dff);
  pr_asym = 1-fcdf(F_iof,dfr-dff,dff);
  pr_boot = NaN;
  iof_stats = [F_iof,pr_asym,pr_boot,sser,ssef,dfr,dff]';
  
  [F_slopes,df,pr_asym] = anovaparam(b(2,:),s(2,:),n); % Slope statistics
  pr_boot = NaN;
  slope_stats = [F_slopes,pr_asym,pr_boot,df(1),df(2)]';
  
  [F_intcpt,df,pr_asym] = anovaparam(b(1,:),s(1,:),n); % Intercept statistics
  pr_boot = NaN;
  intcpt_stats = [F_intcpt,pr_asym,pr_boot,df(1),df(2)]';
  
  if (iter)                                   % Random permutations 
    F_iof_null = zeros(iter,1);
    F_slopes_null = zeros(iter,1);
    F_intcpt_null = zeros(iter,1);
    
    for it = 1:iter
      g = g(randperm(length(g)));               % Randomly permute group membership vector
  
      ssef = 0;
      for i = 1:ngrps                           % Regressions by group
        ig = find(g==grp_id(i));
        [b(:,i),stats,p,r,s(:,i)] = linregr(x(ig),y(ig));
        s(:,i) = s(:,i)*sqrt(n(i));
        ssef = ssef + stats(2)*stats(5);
      end;
      
      F_iof_null(i) = (((sser-ssef)-(dfr-dff))/(dfr-dff))/(ssef/dff);
      F_slopes_null(it) = anovaparam(b(2,:),s(2,:),n); % Slope statistics
      F_intcpt_null(it) = anovaparam(b(1,:),s(1,:),n); % Intercept statistics
    end;
    iof_stats(3) = randprob(F_iof,sort(F_iof_null));
    slope_stats(3) = randprob(F_slopes,sort(F_slopes_null));
    intcpt_stats(3) = randprob(F_intcpt,sort(F_intcpt_null));
  end;
  
  return;
  