%
% Applying the Wilcoxon signed-rank test to data with ties
% Using function  getmidranks()  to find midranks
%

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, Abs(Z_i), and 1 if Z_i>0, 0 if Z_i<0');

ry = size(dd,2);   % Last column = Column number of Y
dd = [dd  zz  abs(zz)  sgn];

dd                 % Display matrix

% Use function getmidranks() to add midranks
zval=dd(:,5);      % Column 5 is Abs(Z)
[midr tiegr] = getmidranks(zval);

%% Add midranks to data matrix
dd = [dd  midr];

disp('DATA WITH MIDRANKS (from getmidranks()):');
disp('Columns are now Ordinal, X, Y, Z, Abs(Z), Sgn, and Midrank');

dd                 % Display matrix

disp('DATA SORTED BY ABS(Z) (with Z then Ord as tie-breakers):');
disp('Columns are Ordinal, X, Y, Z, Abs(Z), Sgn, and Midrank');

ddx = sortrows(dd, [5 4 1])    % Define and display matrix

disp('The same table, using fprintf() instead of a simple display:');
disp(' ');

%% Display ddx in a formatted way row by row
%% Enter ddx', NOT ddx  %'
fprintf('  %2d.  %3d  %3d  %4d  %3d    %d  %5.1f\n', ddx');  %'

disp(' ');
disp('WILCOXON SIGNED RANK TEST:');
disp('Have ties, so use midranks for Tplus and Tminus.');
disp('Recall Tminus+Tplus=nn*(nn+1)/2, even with midranks');

Tplus=0;  Tminus=0;  Tmean=nn*(nn+1)/4;
midr=dd(:,7);    % Column 7 is midranks
sgn=dd(:,6);     % Column 6 is Sgn
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]);

disp('Tie correction depends on distribution of tie groups');
disp('Tie group sizes for size > 1:');

%% The tie correction depends on the following Tiesum
%% Also find the number of tie groups
Tiesum=0;  ntg=0;
for i=1:nn;  tt=tiegr(i);
  if tt>1  ntg=i;  Tiesum=Tiesum+(tt-1)*tt*(tt+1);  end;  end;

disp(' ');
disp(tiegr(1:ntg)');    %' Display as row vector

% Display Tiesum and Tie Correction
TieCorrection = Tiesum/(2*nn*(nn+1)*(2*nn+1));
% Incorrect value TieCorrection = Tiesum/((nn-1)*nn*(nn+1));
fprintf('Tiesum=%g   TieCorrection=%6.4f\n', Tiesum,TieCorrection);

% 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:');

Nsims=10000;
fprintf('\nUsing NSIMS=%d random permutations:\n', Nsims);

disp(' ');
disp('INITIALIZING the random-number generator at');
Wseed=sum(100*clock());
rand('state',Wseed);
fprintf('  Wseed=%6d\n', 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

TTplus = 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 TTplus

kk=(1:nn);      % Fixed row vector 1:nn
NGrEq=0;        % Count of simulations with Trand >= TTplus
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>=TTplus NGrEq=NGrEq+1;  end;
  end;   % End of simulation

disp(' ');
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;

% Sort and display the Walsh averages
ww=sort(ww(1:nwa));

nrows=ceil(nwa/10);
fprintf ('\nSORTED Walsh 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;


% Display Zmean and W_Med
Zmean=sum(zz)/nn;   % Mean of zz
MedZ = median(zz);
MedWAv = median(ww(1:nwa));
Zmean_MedZ_MedWav = [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');  %'

