
disp('Testing Parallelism for Several Lines');
disp('(See Section 9.5 in text (attributed to Sen, Adichie))');
disp('The general model is for `parallel regressions`');
disp(' ');
disp('    Y_ij = beta_i X_ij + mu_i + error_ij         (1)');
disp(' ');
disp('where 1 le i le kk are different TREATMENT GROUPS or BLOCKS');
disp('and 1 le j le nn_i are OBSERVATIONS WITHIN THE i-th GROUP.');
disp('(A variant of (1) assumes mu_i=mu across groups.');
disp(' ');
disp('The hypothesis to test is');
disp(' ');
disp('    H_0: beta_i = beta_1 = same for all blocks   (2)');
disp(' ');
disp('This tests equality of incremental responses for a convariate');
disp('  across different groups or treatment groups.');
disp(' ');
disp('In traditional linear model theory (with normal errors), the model');
disp(' ');
disp('    Y_ij = beta X_ij + mu_i + error_ij           (3)');
disp(' ');
disp('(with constant beta) is called an ANALYSIS OF COVARIANCE');
disp('(or ANOCOVA or ANCOVA) model. The name comes from that fact that');
disp('the original methods for estimating beta and the error variance');
disp('in (3) depended on sample covariances of Y_ij vs X_ij.');
disp(' ');
disp('  Nowadays, linear models with normal errors are analyzed by');
disp('setting up and inverting large ``design`` matrices, with the same');
disp('methods being used for all linear models. Thus the term');
disp('`analysis of covariance` is somewhat obsolete.');
disp(' ');
disp('  In linear models, testing (2) for the more general ANOCOVA model (1)');
disp('with treatment-group-dependent beta_i is called testing for an');
disp('``interaction`` between the slope beta and the treatment groups.');

disp(' ');
disp('EXAMPLE:');
disp('Given kk=4 core samples from Apalachicola Bay, Florida,');
disp('  measure the ammonia released (Y_ij) at times X_ij');
disp('See Table 9.4 page 431 in text');

% X data for k=4 groups with n=5 observations per treatment
xx = [ 0   1.5   3   4.5   6 ; 
       0   1.5   3   4.5   6 ; 
       0   1.5   3   4.5   6 ; 
       0   1.5   3   4.5   6  ];

% Ammonia flux data: Y values
yy = [ 0    33.019  111.314  196.205  230.659 ; 
       0   131.831  181.603  230.070  258.119 ; 
       0    33.351   97.463  196.615  217.308 ; 
       0     8.959  105.384  211.392  255.105  ];

[kk nn] = size(yy);
fprintf('\nkk=%d groups and nn=%d (X,Y) pairs per group\n', kk,nn);

disp(' ');
disp('Data:');
for i=1:kk;
  fprintf(' Core#%d:  Y:',i);
  for j=1:nn;  fprintf('   %7.3f', yy(i,j));  end;  fprintf('\n');
  fprintf('          X:');
  for j=1:nn;  fprintf('   %5.1f  ', xx(i,j));  end;  fprintf('\n');
  end;

disp(' ');
disp(' ');
disp('APPROACH I:');
disp('Usual the model described above.');
disp('The H_0 model is  Y_ij = beta X_ij + mu_i + error_ij.');
disp('  Minimizing over the mu_i before the beta shows that the');
disp('  least-squares estimate of beta is the same as for the model');
disp(' ');
disp('     Y_ij - Ybar_i = beta (X_ij - Xbar_i) + error_ij       (4)');
disp(' ');
disp('  Thus, given H_0, the least-squares estimate of beta is the same');
disp('  as for a single linear regression for all (i,j) pairs');
disp('  with covariate and response variables centered at group means.');
disp(' ');
disp('The Sen-Adichie approach is to compare ``aligned`` values');
disp('  Ystar_ij = Y_ij - beta X_ij  for the value of beta estimated');
disp('  from (4).');
disp(' ');
disp('If H_0 is false, then Ystar_ij should tend to decrease in j');
disp('  for fixed i if beta_i<beta and increase if beta_i>beta.');
disp('  The values Ystar_ij are not corrected for the intercepts mu_i,');
disp('  but within-group Friedman-like ranks R_ij = Rank(Ystar_ij)');
disp('  (so that 1 le R_ij le nn) should show the same correlation.');
disp('  Using rank-regression-like statistics');
disp(' ');
disp('     T_i = Sum(j=1,nn) (X_ij-Xbar_i)*R_ij(Ystar)');
disp(' ');
disp('  should pick up this effect more strongly: That is, T_i should');
disp('  be large positive if beta_i<beta and large negative if beta_i>beta');
disp('  The Sen or Sen-Adichie test is based essentially on the sum');
disp('  of the squares of the T_i.');

% The first step is to define  Xc_ij = X_ij - Xbar_i 
xxc=zeros(kk,nn);   % Same shape as xx
for i=1:kk;   % For each row or group
  % Find the mean of the i-th group (the i-th row mean)
  xsum=0;  for j=1:nn;  xsum=xsum+xx(i,j);  end;
  xmean=xsum/nn;
  % Let xxc(i,j) be xx(i,j) centered at the row mean
  for j=1:nn;  xxc(i,j)=xx(i,j)-xmean;  end;
  end;
% Can also be done by
%      xm = xx*ones(nn,nn)/nn;   % kkxnn matrix of row means
%      xxc = xx - xm;            % rows centered at row means
% Estimate of overall beta = beta_bar
bsum=0;  xsum=0;  % Numerator and denominator
for i=1:kk;  for j=1:nn;
  bsum=bsum+xxc(i,j)*yy(i,j);
  xsum=xsum+xxc(i,j)*xxc(i,j);
  end;  end;
betabar = bsum/xsum;
% Can also say  betabar = sum(sum(xxc.*yy))/sum(sum(xxc.*xxc));
%   sum(xx)  gives a 1xnn row vector of column sums
%   sum(sum(xx))  sums all the entries of a kkxnn matrix xx
fprintf('\nOverall beta (beta_bar) = %.3f\n',betabar);
% For comparison, compute the within-group betas as well
% These are not used in the test procedure
bbw=zeros(kk,1);
for i=1:kk;
  % The traditional least-squares regression slope
  bsum=0;  xsum=0;
  for j=1:nn;
    bsum=bsum+xxc(i,j)*yy(i,j);
    xsum=xsum+xxc(i,j)*xxc(i,j);
    end;
  bbw(i)=bsum/xsum;  end;
fprintf('Beta values within groups:  %.3f  %.3f  %.3f  %.3f\n', bbw);

% Define ``aligned observations``
yystar = yy - betabar*xx;
fprintf('\nAligned observations  Y_ij - betabar X_ij:\n');
for i=1:kk;
  fprintf(' Core#%d:',i);
  for j=1:nn;  fprintf('   %7.3f', yystar(i,j));  end;  fprintf('\n');
  end;

% Find within-group ranks for yystar
yvals=zeros(nn,1);  % Column-vector buffer for rows
rr=zeros(kk,nn);    % kkxnn matrix of within-group (row) ranks
for i=1:kk;
  % Copy a row (group) of yystar to column-vector yvals
  for j=1:nn;  yvals(j)=yystar(i,j);  end;
  % Find within-group midranks of column-vector yvals
  [midr tiegr] = getmidranks(yvals);
  % Copy midranks to the i-th row of rr
  for j=1:nn;  rr(i,j)=midr(j);  end;
  end;
% Within-group ranks are listed in the next table

% Under permutation, ranks are approximately uniformly distributed
%   with range (1,nn), mean (nn+1)/2, and variance approx nn^2/12
% Define two column vectors for row (within-group) values
tt=zeros(kk,1);  ci2=zeros(kk,1);
for i=1:kk;
    tsum=0;  csum=0;
    for j=1:nn;
      tsum=tsum+xxc(i,j)*rr(i,j);
      csum=csum+xxc(i,j)*xxc(i,j);
      end;
    ti(i) = tsum/(nn+1);
    ci2(i) = csum;
    end;
disp(' ');
disp('Within-group ranks, Ti=Sum_j Xc(i,j)*Rank(i,j), and');
disp('   Ci2=Sum(j=1,nn_i) (X_ij - Xbar_i)^2:\n');
for i=1:kk;
  fprintf(' Core#%d:',i);
  for j=1:nn;  fprintf('   %3.1f', rr(i,j));  end;
  fprintf('    Ti=%5.2f   Ci2=%5.2f\n', ti(i), ci2(i));
  end;

% Sen-Adichie`s test:
% Will have a chi-square distribution given H_0
vsum=0;
for i=1:kk;  vsum=vsum+ti(i)*ti(i)/ci2(i);  end;
vv = 12*vsum;
df=kk-1;
fprintf('\nSen-Adichie test:  X=%6.4f   P=%6.4f  (df=%d)\n',...
  vv, pvchisq(vv,df), df);

disp(' ')
disp('Compare with a permutation test of within-group ranks:');
disp('(Friedman-like permutations.)');
disp(' ');
disp('Initializing the random-number generator at');
Startseed=12345678;
fprintf('    %d\n',Startseed);
rand('state',Startseed);

nsims=10000;
fprintf('\nCarrying out ns=%d sets of within-group permutations:\n',nsims);

Ngr=0;
% A computationally-friendly version of the statistic vv
% These tt() are like ti() but are not divided by nn+1
tt=zeros(nn,1);
for i=1:kk;  for j=1:nn;
  tt(i)=tt(i)+xxc(i,j)*rr(i,j);  end;  end;
tscore=0;
for i=1:kk;  tscore=tscore+tt(i)*tt(i)/ci2(i);  end;
% Can also say
%   tt=sum(xxc.*rr,2);
%   tscore=sum(tt.*tt);
fprintf('Observed Tscore = %g\n',tscore);
% Define two buffers
rstar=zeros(kk,nn);  rvals=zeros(1,nn);
for ns=1:nsims;
    for i=1:kk;
        %% Randomly permute the i-th row
        rstar(i,:) = rr(i,randperm(nn));  end;
    % Since this is within a loop, use the faster MATLAB code
    % sum(xx,2) for a kkxnn matrix gives a kkx1 vector of row sums
    tt=sum(xxc.*rstar,2);
    ttstar=sum(tt.*tt./ci2);
    if ttstar>=tscore Ngr=Ngr+1;  end;
    end;

disp('Number of simulations with values >= Tscore and total number:');
disp(' ');
disp([Ngr nsims]);

phat=Ngr/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
disp('95% CI for true P-value bracketing estimate of true pvalue:');
disp('(Since H and Fscore >= 0, P-values are inherently two-sided.)');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);


disp(' ');
disp(' ');
disp('APPROACH II:');
disp('Since the within-group regressions appear to have no intercepts,');
disp('  it makes more sense to apply the same arguments with');
disp('  no-intercept regressions and ignore the observations');
disp('  with j=1 and X_ij=Y_ij=0');

% Delete the first column in  X and Y  and adjust nn
xx1 = xx(:,2:nn);
yy1 = yy(:,2:nn);
nn1 = nn-1;

fprintf(...
  'IGNORING the first columns of X and Y: using %dx%d matrices X1,Y1\n',...
  kk,nn1);

% Estimate of overall beta = beta_bar
% The following is the least-squares estimate of the slope for a 
%   no-intercept regression, which are slightly different than for a
%   regression with an intercept
bsum=0;  xsum=0;  % Numerator and denominator
for i=1:kk;  for j=1:nn1;
  bsum=bsum+xx1(i,j)*yy1(i,j);
  xsum=xsum+xx1(i,j)*xx1(i,j);
  end;  end;
betabar = bsum/xsum;

disp(' ');
disp('Using no-intercept regressions:');
fprintf('Overall beta (beta_bar) = %.3f\n',betabar);
% For comparison, compute the within-group betas as well
% These are not used in the test procedure
bbw=zeros(kk,1);
for i=1:kk;
  % The traditional least-squares regression slope
  bsum=0;  xsum=0;
  for j=1:nn1;
    bsum=bsum+xx1(i,j)*yy1(i,j);
    xsum=xsum+xx1(i,j)*xx1(i,j);
    end;
  bbw(i)=bsum/xsum;  end;
fprintf('Beta values within groups:  %.3f  %.3f  %.3f  %.3f\n', bbw);

% Define ``aligned observations``
yystar1 = yy1 - betabar*xx1;
fprintf('\nAligned observations  Y1_ij - betabar X1_ij:\n');
for i=1:kk;
  fprintf(' Core#%d:',i);
  for j=1:nn1;  fprintf('   %7.3f', yystar1(i,j));  end;  fprintf('\n');
  end;

% Find within-group ranks for yystar
yvals=zeros(nn1,1);   % Column-vector buffer for rows of Ystar
rr1=zeros(kk,nn1);    % kkxnn1 matrix of within-group (row) ranks
for i=1:kk;
  % Copy a row (group) of yystar to column-vector yvals
  for j=1:nn1;  yvals(j)=yystar1(i,j);  end;
  % Find within-group midranks of column-vector yvals
  [midr tiegr] = getmidranks(yvals);
  % Copy midranks to the i-th row of rr
  for j=1:nn1;  rr1(i,j)=midr(j);  end;
  end;
% Within-group ranks are listed in the next table

% We still need centered X values Xc1 for the Ti statistics
xxc1=zeros(kk,nn1);
for i=1:kk;
    xs=0;  for j=1:nn1;  xs=xs+xx1(i,j);  end;
    xmean=xs/nn1;
    for j=1:nn1;  xxc1(i,j) = xx1(i,j) - xmean;  end;
    end;

% Carry out the same procedure as before with Xc,R,nn
%   replaced by Xc1,R1,nn1
tt=zeros(kk,1);  ci2=zeros(kk,1);
for i=1:kk;
    tsum=0;  csum=0;
    for j=1:nn1;
      tsum=tsum+xxc1(i,j)*rr1(i,j);
      csum=csum+xxc1(i,j)*xxc1(i,j);
      end;
    ti(i) = tsum/(nn1+1);
    ci2(i) = csum;
    end;
disp(' ');
disp('Within-group ranks, Ti=Sum_j Xc1(i,j)*Rank(i,j)/(nn1+1), and');
disp('   Ci2=Sum(j=1,nn_i) (X1_ij - X1bar_i)^2:\n');
for i=1:kk;
  fprintf(' Core#%d:',i);
  for j=1:nn;  fprintf('   %3.1f', rr(i,j));  end;
  fprintf('    Ti=%5.2f   Ci2=%5.2f\n', ti(i), ci2(i));
  end;

% Sen-Adichie`s test:
% Will have a chi-square distribution given H_0
vsum=0;
for i=1:kk;  vsum=vsum+ti(i)*ti(i)/ci2(i);  end;
vv = 12*vsum;
df=kk-1;
fprintf('\nSen-Adichie test:  X=%6.4f   P=%6.4f  (df=%d)\n',...
  vv, pvchisq(vv,df), df);

disp(' ')
disp('Compare with a permutation test of within-group ranks:');
disp('(Friedman-like permutations.)');
disp(' ');
disp('Initializing the random-number generator at');
Startseed=12345678;
fprintf('    %d\n',Startseed);
rand('state',Startseed);

nsims=10000;
fprintf('\nCarrying out ns=%d sets of within-group permutations:\n',nsims);

Ngr=0;
% A computationally-friendlier version of the statistic vsum
% These tt() are like ti() but are not divided by nn1+1
tt=zeros(nn1,1);
for i=1:kk;  for j=1:nn1;
  tt(i)=tt(i)+xxc1(i,j)*rr1(i,j);  end;  end;
tscore=0;
for i=1:kk;  tscore=tscore+tt(i)*tt(i)/ci2(i);  end;
% Can also say
%   tt=sum(xxc1.*rr1,2);
%   tscore=sum(tt.*tt./ci2);

fprintf('Observed Tscore = %g\n',tscore);
% Define two buffers
rstar1=zeros(kk,nn1);  rvals=zeros(1,nn1);
for ns=1:nsims;
    for i=1:kk;
        %% Randomly permute the i-th row
        rstar1(i,:) = rr1(i,randperm(nn1));  end;
    % Since this is within a loop, use the faster MATLAB code
    % sum(xx,2) for a kkxnn1 matrix gives a kkx1 vector of row sums
    tt=sum(xxc1.*rstar1,2);
    ttstar=sum(tt.*tt./ci2);
    if ttstar>=tscore Ngr=Ngr+1;  end;
    end;

disp('Number of simulations with values >= Tscore and total number:');
disp(' ');
disp([Ngr nsims]);

phat=Ngr/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
disp('95% CI for true P-value bracketing estimate of true pvalue:');
disp('(Since H and Fscore >= 0, P-values are inherently two-sided.)');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);


disp(' ');
disp(' ');
disp('APPROACH III:');
disp('Since within-group regressions appear to have no intercepts,');
disp('  and X_ij>0 except for j=1 with X_ij=Y_ij=0, one could also apply');
disp('  the Kruskal-Wallis test for  Z_ij=Y_ij/X_ij  for 2 le j le nn');

%% Set Z_ij=Y_ij/X_ij
nn1=nn-1;   % As in Approach II
zz=zeros(kk,nn1);
for i=1:kk;  for j=2:nn;
  zz(i,j-1) = yy(i,j)/xx(i,j);  end;  end;
disp('Data for Y, X, and Z for X_ij>0:');
for i=1:kk;
  fprintf(' Core#%d:  Y:',i);
  for j=2:nn;  fprintf('   %7.3f', yy(i,j));  end;  fprintf('\n');
  fprintf('          X:');
  for j=2:nn;  fprintf('   %5.1f  ', xx(i,j));  end;  fprintf('\n');
  fprintf('          Z:');
  for j=1:nn1;  fprintf('   %7.3f', zz(i,j));  end;
  zsum=0;  for j=1:nn1;  zsum=zsum+zz(i,j);  end;
  fprintf('   Average: %7.3f\n', zsum/nn1);
  end;

% Find Kruskall-Wallis ranks for Z_ij
zvals=zeros(kk*nn1);  r=0;
for j=1:nn1;  for i=1:kk;
  r=r+1;  zvals(r)=zz(i,j);  end;  end;
% zvals=reshape(zz,kk*nn1,1);  gives exactly the same vector
[midr tiegr] = getmidranks(zvals);
% Let rr=kkxnn1 matrix of Kruskal-Wallis ranks
rr=zeros(kk,nn1);  r=0;
for j=1:nn1;  for i=1:kk;
  r=r+1;  rr(i,j)=midr(r);  end;  end;
% rr=reshape(midr,kk,nn1);  gives exactly the same matrix
disp(' ');
disp('Kruskal-Wallis ranks, rank sums, and rank means:');
rksum=zeros(kk,1);  rkmean=zeros(kk,1);
for i=1:kk;
  fprintf(' Core#%d:  Y:',i);
  for j=1:nn1;  fprintf('   %4.1f', rr(i,j));  end;
  rsum=0;  for j=1:nn1;  rsum=rsum+rr(i,j);  end;
  rksum(i)=rsum;  rkmean(i)=rsum/nn1;
  fprintf('    Rsum=%5.2f   Rmean=%6.3f\n', rksum(i),rkmean(i));
  end;
% Count number of ties
zg=0;
for i=1:kk*nn1;
  if tiegr(i)>0 zg=i;  else break;  end;  end;
if zg==0  disp('No ties.');
else fprintf('%d tie group(s)\n',zg);  end;

disp('  ');
disp('Large-sample Kruskal-Wallis approximation (ASSUMING no ties):');
num=kk*nn1;
rkmeanc=rkmean-(num+1)/2;
hh=0;
for i=1:kk;  hh=hh+rkmeanc(i)*rkmeanc(i);  end;
% Could also say  hh=sum(rkmeanc.*rkmeanc);
hh=(12/(num*(num+1)))*nn1*hh;
df=kk-1;
pval=pvchisq(hh,df);
fprintf('    H=%6.4f   P=%6.4f  (df=%d)\n', hh,pval,df);


disp('  ');

