
disp('One-Way Layouts:');
disp('Kruskal-Wallis test for equality of location parameters');
disp('  Model:  X(i,j) dist X + theta_i ---  Test H_0: Theta are identical');

disp('  ')
disp('Example 1:');
disp('Treatment #1: Control   #2: Obstruc.Air.Disease   #3: Asbestosis');
disp('Values relate to efficiency of ciliary dust clearance.');
disp('Table 6.1 in text page 192');

disp(' ');
disp('Data in 3 treatment groups with sizes 5,4,5');

x1 = [ 2.9  3.0  2.5  2.6  3.2 ];
x2 = [ 3.8  2.7  4.0  2.4 ];
x3 = [ 2.8  3.4  3.7  2.2  2.0 ];

n1=length(x1);  n2=length(x2);  n3=length(x3);

disp(' ');
fprintf('  Trtment #1 (n=%d) ',n1);
for j=1:n1;  fprintf ('  %3.1f', x1(j));  end;  fprintf('\n');
fprintf('  Trtment #2 (n=%d) ',n2);
for j=1:n2;  fprintf ('  %3.1f', x2(j));  end;  fprintf('\n');
fprintf('  Trtment #3 (n=%d) ',n3);
for j=1:n3;  fprintf ('  %3.1f', x3(j));  end;  fprintf('\n');

disp(' ');
disp('Number of treatment groups and total number of observations:');
kk=3;  nn=n1+n2+n3;
disp(' ');
disp([kk nn]);

disp('Assemble in a vertical array for analysis:');
disp('(Columns are Ordinal, Value, Treatment group, Midrank, Tiegroups');

dd = [ x1' ones(n1,1);  x2'  2*ones(n2,1);  x3'  3*ones(n3,1)];  %'

nn=size(dd,1);   % Sum n1+n2+n3;

[midr tiegr] = getmidranks(dd);   % Using 1st column of dd;

dd = [ [1:nn]'  dd  midr  tiegr];  %'

fprintf('   %2d  %4.2f   %d  %5.1f  %2d\n', dd');  %'

disp(' ');
disp('Parallel display with midranks:');
r1 = midr(1:n1,1);  r2=midr(n1+1:n1+n2,1);  r3=midr(n1+n2+1:nn,1);
for i=1:5;
  fprintf('  %4.1f (%4.1f)', x1(i), r1(i));
  if (i<5)  fprintf('   %4.1f (%4.1f)', x2(i), r2(i));
  else  fprintf('              ');  end;  % Of if statement
  fprintf('   %4.1f (%4.1f)\n', x3(i), r3(i));
  end;  % Of for-loop

disp(' ');
% Compute Kruskal-Wallis H statistic
% First, accumulate rank sums and group numbers by row
ranksum=zeros(kk,1);  num=zeros(kk,1);
group=dd(:,3);  rank=dd(:,4);
for i=1:nn;  ranksum(group(i)) = ranksum(group(i)) + rank(i);
    num(group(i)) = num(group(i)) + 1;  end;
% ranksum and num are vectors, so we can divide vectors
rankav=ranksum./num;

disp('By Treatment Group: Number, Rank Sum, Rank Average:');
rdd = [ [1:kk]'  num  ranksum  rankav];  %'
fprintf('  Group #%d  n=%d  %5.1f %6.2f\n', rdd');  %'

% Compute Kruskal-Wallis H statistic
hh=0;
ravc = rankav - (nn+1)/2;
for k=1:kk;  hh = hh + num(k)*ravc(k)*ravc(k);  end;
hh = (12*hh)/(nn*(nn+1));

% Alternatively:  hh = (12/(nn*(nn+1))*sum(num.*ravc.*ravc);

fprintf('\nKruskal-Wallis statistic:  H = %g\n', hh);

disp(' ');
disp('Large-sample chi-square approximation for H:');
disp('No ties, so no tie correction needed')

df=kk-1;
Pval = pvchisq(hh,df);
fprintf('\n  P= %6.4f  (df=%d)\n', Pval,df);

disp(' ');
disp('Doing nsims simulations for the true P-value');
nsims=10000

% Initialize random-number generator from the system clock
disp('Initializing the RNG from the system clock at');

startseed=sum(100*clock)

rand('state',startseed);

% Note that H = A_n * HScore + C_n   for  A_N=12/(nn*(nn+1)), C_N=-3(nn+1)
%   for  Hscore = Sum ranksum^2/num
% Thus simulating  Hscore  is the same as simulating  H (hh)

% Observed Hscore:
Hscore=0;
for k=1:kk;  Hscore=Hscore + ranksum(k)*ranksum(k)/num(k);  end;
 
% Faster:  Hscore=sum(ranksum.*ranksum./num);

H_and_Hscore = [hh  Hscore]

% Do nsims permutations for a test of Hscore
Ngr=0;  group=dd(:,3);  midr=dd(:,4);
for ns=1:nsims;
  %% Permuted copy of midranks:
  midrperm=midr(randperm(nn));
  %% Store rank sums for permuted midranks in rrsum
  rrsum=zeros(kk,1);
  for i=1:nn;  rrsum(group(i))=rrsum(group(i))+midrperm(i);  end;
  %% Calculate the hscore for the permuted midranks
  Hrand = sum(rrsum.*rrsum./num);
  if Hrand >= Hscore  Ngr=Ngr+1;  end;
  end;

disp('Number of simulations with values >= Hscore 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 Hscore >= 0, P-values are inherently two-sided.)');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);

disp(' ');
disp('Test for an alternative of monotonic treatment-group medians:');
disp('The Jonckheere-Terpstra test statistic is very similar to');
disp('  the Kendall nonparametric correlation statistic between');
disp('  values and treatment-group numbers');
% Put values and group numbers in column vectors
val=dd(:,2);  group=dd(:,3);
% Compute the Jonckheere-Terpstra score
jj=0;  
for i=1:nn;  for j=1:nn;
  g1=group(i);  g2=group(j);
  if g1<g2  val1=val(i);  val2=val(j);
    if val1<val2 jj=jj+1;
    elseif val1==val2 jj=jj+0.50;  end;  end;  % For nested ifs
  end;  end;  % For double loop

%% Display the statistic J
jj

disp('The Jonckheere-Terpstra test statistic mean, variance,');
disp('  Z-statistic, and one-sided normal P-value (assuming no ties):');

jmean = (nn*nn - sum(num.*num))/4;
jvar = (nn*nn*(2*nn+3) - sum(num.*num.*(2*num+3)))/72;

zz = (jj-jmean)/sqrt(jvar);

Jmean_Jvar_Zval_NormPval = [jmean  jvar  zz  normpval(zz)/2]




