%
% Applying the Wilcoxon signed-rank test to data in text

disp('Data in Table 3.1 page 39');

dd = [ 1.  1.830  0.878;
       2.  0.500  0.647;
       3.  1.620  0.598;
       4.  2.480  2.050;
       5.  1.680  1.060;
       6.  1.880  1.290;
       7.  1.550  1.060;
       8.  3.060  3.140;
       9.  1.300  1.290  ];

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

disp('EXPANDED DATA MATRIX:');
disp('Columns are Patient number, Before (X_i), After (Y_i),');
disp('  Z_i=Y_i-X_i, and 1 if Z_i>0, 0 if Z_i<0');

dd = [dd  zz  sgn]

disp('SIGN-TEST ANALYSIS:');
disp('  ');

% Set N=Sample size, B=Sign Test Score, 2-sided binomial P-value

nn=size(dd,1);   % Each row is a patient
bb=sum(sgn);     % Number rows with After>Before

% Is bb on the upper or lower tail of the binomial distribution?
% If on the upper tail, replace bb by nn-bb on the lower tail and use
%    P(Binom(N,1/2)>=B) = P(Binom(N,1/2)<=N-B)
% Then set  Pval = 2*P(Binom(N,1/2) >= N-B);

bb2=min(bb,nn-bb);
bpval = 2*binomtail(nn-bb2,nn,0.50);

disp('N=Sample Size, B=Sign Test Score, P=2-sided binomial P-value:');
N_B_Pval = [nn  bb  bpval]      % Display

% Large-sample normal approximation for P-value?
% If X is binom(n,p), then  Z = (X-np)/sqrt(n*p*(1-p))
%   is approximately N(0,1)

Zval = (bb-(nn/2))/sqrt(nn*0.50*0.50);
Pval = normpval(Zval);

disp('Large-sample normal-approximation test for binomial probability:');
Zval_Pval = [Zval  Pval]      % Display
disp('Note that the large-sample approximation for the sign test');
disp('  is not very accurate.');

disp('  ');
disp('The HODGES-LEHMANN estimator for theta for After=Before+theta')
disp('  for the sign test is the median of the Z_i values:')
disp('Zmean is the mean of the Z=Y-X values.');
% Zmean for comparison
Zmean=sum(zz)/nn;
B_Med = median(zz);

Zmean_BMed = [Zmean B_Med]

disp('  ');
disp('WILCOXON SIGNED RANK TEST:');
disp('  ');
disp('Drop X and Y columns, add Abs(Z)');
disp('Columns are Ordinal, Z, Sgn(Z), and Abs(Z)');

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

disp('SORT the matrix by its FOURTH column = Abs(Z):');
disp('Then ADD a FIFTH column that are the signed ranks');

ddx = sortrows(ddx,4);
ddx = [ddx  (1:nn)']     %'

disp('Data in ORIGINAL ORDER (by sorting on Pat.Number):');
disp('Columns are Pat.Number, Z, Sgn, Abs(Z), and Rank');

ddx = sortrows(ddx,1)

% Calculate T+
Tplus=0;
sgn=ddx(:,3);     % Signs are Column 3
rnk=ddx(:,5);     % Ranks are Column 5
for i=1:nn;  if sgn(i)==1 Tplus=Tplus+rnk(i);  end;  end;

% We could also have said  Tplus=sum(sgn.*rnk);

% Find the mean and standard deviation of Tplus 
%   under random permutations assuming no ties
Tmean = nn*(nn+1)/4;
Tstd = sqrt(nn*(nn+1)*(2*nn+1)/24);

Tplus_Tmean_Tstd = [Tplus  Tmean  Tstd]

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

Zval = (Tplus-Tmean)/Tstd;
Pval = normpval(Zval);

Zval_Pval = [Zval  Pval]

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 upper tail, replace Tplus with the corresponding
%   value on the LOWER TAIL

TTplus = min(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
NLessEq=0;      % Count of simulations with Trand le 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 NLessEq=NLessEq+1;  end;
  end;   % End of simulation

disp('FINISHED NSIMS simulations');

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

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

fprintf('NLessEq=%d  Nsims=%d   Pvalue=%7.5f (2-sided)\n',NLessEq,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('The HODGES-LEHMANN estimator for theta for After=Before+theta')
disp('  for the Wilcoxon signed-rank test is the median of');
disp('  the Walsh averages W_k=(Z_i+Z_j)/2:')
% Collect the Walsh averages 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

% Display Walsh averages unsorted then sorted
nrows=ceil(nwa/10);
fprintf ('\nWalsh averages (not in order; 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+1 > nwa break;  end;  % end of if statement
  for j=r+1:r2;  fprintf(' %6.3f', ww(j));  end;
  fprintf('\n');  end;

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(' %6.3f', ww(j));  end;
  fprintf('\n');  end;


% Show sample mean of Z_i and median of Walsh averages
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=9:
TabA4 = [ 39 0.027;  40 0.020;  41 0.014;  42 0.010;  43 0.006 ];

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('  (%6.3f, %6.3f)  x=%2d  %4.1f%% CI\n', CITab');  %'


