%
% Applying the Wilcoxon signed-rank test to data with ties
%

disp('DATA FOR N=11 paired differences:');

dd = [  1.  26   37 
        2.  25   29 
        3.  59   62 
        4.  27   25 
        5.  34   47 
        6.  38   40 
        7.  51   58 
        8.  48   52 
        9.  39   42 
       10.  37   34 
       11.  30   42 ];

xx=dd(:,2);     % Isolate the 2nd and 3rd columns
yy=dd(:,3);
zz=yy-xx;       % After-Before vector difference
sgn=(xx<yy);    % = 1 if After>Before, = 0 if After<Before
                % (Condition)=1 if TRUE, =0 if FALSE,
                %   and is vectorized over the 9 patients
nn=size(dd,1);  % Number of rows = number of observations

disp('Expanded data matrix:');
disp('Columns are Ordinal, Before (X_i), After (Y_i),');
disp('  Z_i=Y_i-X_i, Sgn(Z_i>0), and Abs(Z_i)');

dd = [dd  zz  sgn  abs(zz)];

dd                 % Display matrix

disp('DATA SORTED BY ABS(Z) (with Z then Ord as tie-breakers):');
disp('After dropping columns X and Y');
disp('Add Tentative Ranks after sort');
disp('Columns are now Ordinal, Z, Sgn, Abs(Z), and Tent.Rank');

ddx = [dd(:,1)  zz  sgn  abs(zz)];

ddx = sortrows(ddx, [4 2 1]);

ddx = [ddx  (1:nn)']      %' Define and display matrix

% Find ranks and midranks for data in sorted order
% Abs(ZZ) and Midranks
zval=ddx(:,4);      % Sorted Abs(Z)
midr=zeros(nn,1);   % Buffer for midranks as nnx1 column vector
tiegr=zeros(nn,1);  % Buffer for sizes of tie groups with tsize>1
% Initial values for four variables
% Abs(Z) is sorted, so unique initial value
oldval=zval(1)-1;  g=0;  tsize=-1;  firsti=-1;
% Loop over possibly-tied Abs(Z) values
% Note zval(1)!=oldval
for i=1:nn;
  if (zval(i)==oldval) tsize=tsize+1;  % Same tie group, so add one
  else
    % Beginning of tie group, possibly end of another
    % If tsize>1, record values
    if tsize>1 g=g+1;  tiegr(g)=tsize;  mrval=(firsti+i-1)/2;
      for j=firsti:i-1;  midr(j)=mrval;  end;
    % If tsize==1, only one previous value
    elseif tsize==1  midr(firsti)=firsti;  end;
    % Initialize for the next tie group
    tsize=1;  oldval=zval(i);  firsti=i;  end;
  end;
%% The loop leaves the last tie group open
%% Thus must handle it as well
if tsize>1 g=g+1;  tiegr(g)=tsize;  mr=(firsti+nn-1)/2;
  for j=firsti:nn;  midr(j)=mr;  end;
elseif tsize==1 midr(nn)=nn;  end;
ntg=g;   % More descriptive number for number of tie groups

% List tie-group sizes as a row vector
Tiegroup_Sizes = tiegr(1:ntg)'     %'

% Add midranks to ddx and display

disp('COMPUTED MIDRANKS IN SORTED ORDER');
disp('Columns are Ordinal, Z, Sgn, Abs(Z), Tent.Rank, and Midrank');
ddx = [ddx  midr]

disp(' ');
disp('MIDRANKS IN THE ORIGINAL DATA ORDER');
% The original ordinal values are still in Column 1

ddx = sortrows(ddx,1)    % Define and display matrix

disp('  ');
disp('WILCOXON SIGNED RANK TEST:');
disp('Have ties, so have to use midranks.');

midr=ddx(:,6);    % Column 6 is midranks
sgn=ddx(:,3);     % Column 3 is Sgn
Tplus=0;  Tminus=0;  Tmean=nn*(nn+1)/4;
for i=1:nn;
  if sgn(i)==1 Tplus=Tplus+midr(i);
  else Tminus=Tminus+midr(i);  end;  end;

disp('  ');
disp('Tminus, Tplus, Tmean,  N');
disp([Tminus Tplus Tmean nn]);

disp(' ');
disp('LARGE-SAMPLE NORMAL APPROXIMATION for two-sided Wilcoxon');
disp('  signed-rank P-value for H_0:theta=0:');

% Variance if no ties (all midranks are ranks)
TvarNoTC = nn*(nn+1)*(2*nn+1)/24;
TstdNoTC = sqrt(TvarNoTC);
ZvalNoTC = (Tplus-Tmean)/TstdNoTC;
NormPvalNoTC = normpval(ZvalNoTC);

disp(' ');
disp('IGNORING tie correction (TstdDev then Z then two-sided normal P):');
tt=sqrt(TvarNoTC);
disp([tt  ZvalNoTC  NormPvalNoTC]);

% The tie correction depends on the following Tiesum
Tiesum=0;
for g=1:ntg;  tsize=tiegr(g);
  Tiesum=Tiesum+(tsize-1)*tsize*(tsize+1);  end;
% Display Tiesum and Tie Correction
Tiesum
TieCorrection = Tiesum/(2*nn*(nn+1)*(2*nn+1))
% Incorrect value TieCorrection = Tiesum/((nn-1)*nn*(nn+1))

% Large-sample normal approximation with tie correction
TvarTC = TvarNoTC*(1 - TieCorrection);
TstdTC = sqrt(TvarTC);
ZvalTC = (Tplus-Tmean)/TstdTC;
NormPvalTC = normpval(ZvalTC);
disp(' ');
disp('CORRECT APPROXIMATE P-VALUE WITH tie correction:');
disp('  (TstdDev then Z then two-sided normal P):');
tt=sqrt(TvarTC);
disp([tt  ZvalTC  NormPvalTC]);


disp('SIMULATING THE TRUE Wilcoxon Signed-Rank P-VALUE:');
disp('  ');
disp('Using NSIMS random permutations:');

Nsims=10000

disp('INITIALIZING the random-number generator at');
Wseed=sum(100*clock())
rand('state',Wseed);

% The distribution of Tplus under permutations is symmetric
%   about Tmean.
% If Tplus is on the lower tail, replace Tplus with the corresponding
%   value on the upper tail

Wplus = max(Tplus, 2*Tmean-Tplus);

% Find the distribution of Tplus under random permutations, each of which
%   here is a set of N independent switchings of Before,After values.
% Find the frequency ratio of random permutations for which
%   the random Tplus >= the observed Wplus

kk=(1:nn);      % Fixed row vector 1:nn
NGrEq=0;        % Count of simulations with Trand >= Wplus
for ns=1:Nsims;    % Loop over Nsims simulations
  rr = floor(2*rand(1,nn));   % rand(1,nn) is a row vector of random U(0,1)s
                              % rr is a 1xnn row vector of random 0s and 1s
  Trand = sum(rr.*kk);        % Componentwise multipication
                              %   yielding Tplus for that permutation
  if Trand>=Wplus NGrEq=NGrEq+1;  end;
  end;   % End of simulation

disp('FINISHED NSIMS simulations');

phat=NGrEq/Nsims;     % one-sided P-value
SimPval = 2*phat;     % two-sided P-value

disp(' ');
disp('SIMULATED value for true 2-sided P-value:')

fprintf('NGrEq=%d  Nsims=%d   Pvalue=%7.5f (2-sided)\n',NGrEq,Nsims,SimPval);

disp(' ');
disp('95%% Confidence Interval for true P-value (with Estimate in middle):');
% phat, not SimPval, is distributed as Binom(N,p)/N 
pwid=1.960*sqrt(phat*(1-phat)/Nsims);
Plow=2*(phat-pwid);  Phigh=2*(phat+pwid);

Plow_PEst_Phigh = [Plow  SimPval  Phigh]

disp(' ');
disp('The Hodges-Lehmann estimator for theta for After=Before+theta')
disp('  is the median of the Walsh averages W_k=(Z_i+Z_j)/2:')
% Collect the Walsh sums in a row vector of length nn*(nn+1)/2
ww=zeros(1,nn*(nn+1)/2);
r=0;
for i=1:nn;  for j=i:nn;  r=r+1;
   ww(r)=(zz(i)+zz(j))/2;  end;  end;
nwa=r;   % Number of Walsh averages

nrows=ceil(nwa/10);
fprintf ('\nWalsh averages (n=%d, nwa=%d, nrows=%d):\n', nn,nwa,nrows);
for i=1:nrows;  r=(i-1)*10;  r2=r+10;
  if r2>nwa r2=nwa;  end;
  if r < nwa fprintf(' (%2d)', r+1);
  else break;  end;  % of if statement
  for j=r+1:r2;  fprintf(' %5.1f', ww(j));  end;
  fprintf('\n');  end;

fprintf ('\nSORTED Walsh averages (n=%d, nwa=%d, nrows=%d):\n', nn,nwa,nrows);
ww=sort(ww(1:nwa));
nrows=ceil(nwa/10);
for i=1:nrows;  r=(i-1)*10;  r2=r+10;
  if r2>nwa r2=nwa;  end;
  if r < nwa fprintf(' (%2d)', r+1);
  else break;  end;  % of if statement
  for j=r+1:r2;  fprintf(' %5.1f', ww(j));  end;
  fprintf('\n');  end;

% Display Zmean and MedWAv
Zmean=sum(zz)/nn;   % Mean of zz
MedZ = median(zz);
MedWAv = median(ww(1:nwa));
Zmean_MedZ_MedWalshAverages = [Zmean  MedZ  MedWAv]

disp('EXACT CONFIDENCE INTERVALS based on the Wilcoxon signed-rank statistic:');

fprintf('\nEstimator: %g\n', MedWAv);

% For exact confidence intervals for n=11:
TabA4 = [ 55 0.027;  56 0.021;  57 0.016;  58 0.012;
          59 0.009;  60 0.007;  61 0.005;  62 0.003 ];

fprintf('\nTable A4 for n=%d  (n*(n+1)/2=%d)  has:\n',nn,nwa);
fprintf ('   P(T+ ge %2d)=%5.3f\n', TabA4');  %'
nTab=size(TabA4,1);

disp(' ');
disp('Exact symmetric confidence intervals for theta:');
xhi=TabA4(:,1);  xlo=nwa + 1 - xhi;

CITab = zeros(nTab,2);
for i=1:nTab;
  CITab(i,1)=ww(xlo(i));  CITab(i,2)=ww(xhi(i));  end;
clevel=100*(1-2*TabA4(:,2));
CITab = [CITab  xhi  clevel];

fprintf('  (%5.1f, %5.1f)  x=%2d  %4.1f%% CI\n', CITab');  %'

