
disp('Theil`s Method for Estimation of Slope and testing');
disp('  H_0:beta=beta_0 for a simple regression');
disp('See Sections 9.1-9.3 in text');

disp(' ');
disp('Data: Amount of pre-packed engine lubricant and engine survival times:');
disp('(Source of data unknown.)');
disp(' ');
disp('Lubricant  Survival          Lubricant Survival');

ww2 = [ 1.   29     23     11.   84    429;
        2.   21     36     12.   37    489;
        3.   30     67     13.   94    504;
        4.   45    104     14.   56    925;
        5.   48    147     15.   67    980;
        6.   92    164     16.   79   1254;
        7.   74    247     17.   70   2961;
        8.   61    304     18.   93   5351;
        9.   99    355     19.   76   9915;
       10.   24    408     20.   72  11301  ];

fprintf (' %3d.  %2d  %5d          %3d.   %2d  %5d\n', ww2');  %'

disp(' ');
disp('Let X=Lubricant and Y=Survival (in days)');

% Rewrite data as one nnx3 array
ww = [ ww2(:,1) ww2(:,2) ww2(:,3);  ww2(:,4) ww2(:,5) ww2(:,6)];
nn=size(ww,1);      % Number of rows

fprintf('\nSorting data into one %dx3 matrix\n', nn);

ww=sortrows(ww,1);  % Sort on ordinal number

% Fetch Xs and Ys as nnx1 column vectors
xx=ww(:,2);  yy=ww(:,3);

% Are there any ties?
[midrx  xties] = getmidranks(xx);
xg=0;
for i=1:nn;  
  if xties(i)>0 xg=i;  else break;  end;  end;
[midry  yties] = getmidranks(yy);
yg=0;
for i=1:nn;  
  if yties(i)>0 yg=i;  else break;  end;  end;
fprintf('\nWithin-sample tie groups:  X: %d   Y: %d\n', xg, yg);
if xg==0 && yg==0  disp('NO TIES in Xs or Ys');  end;
%% NOTE:
%% If there are ties, use the tie correction for Kendall`s test
%%   from Section 8.1 page 366 in text

disp(' ');
disp('Theil`s test of H_0:beta=0  (same as Kendall`s test)');
disp('(NOT assuming that X(i) are increasing, so different from text)');
cc=0;
for i=1:nn;  for j=i+1:nn;
  zz = (yy(j)-yy(i))*(xx(j)-xx(i));
  if (zz>0) cc=cc+1;
  elseif zz<0 cc=cc-1;  end;  end;  end;
varC = nn*(nn-1)*(2*nn+5)/18;
zz = cc/sqrt(varC);
fprintf('  C=%d   Z=%6.4f   P=%6.4f (Assuming no ties)\n', ...
  cc, zz, normpval(zz));

beta_0=5;
fprintf('\nTheil`s test of H_0:beta=%g  (similar to Kendall`s test)\n',beta_0);
%% Let yyx(i) = yy(i) - beta_0*xx(i)
yyx=yy-beta_0*xx;
cc=0;
for i=1:nn;  for j=i+1:nn;
  zz = (yyx(j)-yyx(i))*(xx(j)-xx(i));
  if (zz>0) cc=cc+1;
  elseif zz<0 cc=cc-1;  end;  end;  end;
varC = nn*(nn-1)*(2*nn+5)/18;
zz = cc/sqrt(varC);
fprintf('  C=%d   Z=%6.4f   P=%6.4f (Assuming no ties)\n', ...
  cc, zz, normpval(zz));

%% Compare with least-squares regression estimates
ymean=sum(yy)/nn;  xmean=sum(xx)/nn;
yyc=yy-ymean;    xxc=xx-xmean;
ssxx = sum(xxc.*xxc);
beta_lsq = sum(yyc.*xxc)/ssxx;
mu_lsq = ymean - beta_lsq*xmean;
resid_lsq = yy - beta_lsq*xx - mu_lsq;

disp(' ');
disp('Normal-theory test of H_0:beta=0  (Student-t test)');
mse = sum(resid_lsq.*resid_lsq)/(nn-2);
tt = beta_lsq/sqrt(mse/ssxx);
fprintf('  beta_hat=%6.4f   T=%6.4f   P=%6.4f  (df=%d)\n', ...
  beta_lsq, tt, pvtdist(tt,nn-2), nn-2);
disp('  (For normal-theory test, note huge outliers in Y_i.)');

disp(' ');
disp('SIMPLE REGRESSION OF Y on X  (Y = beta*X + mu)');
disp(' ');
disp('Theil`s estimate of beta is essentially the Hodges-Lehmann');
disp('  estimator of beta based on the Kendall (or Kendall tau) test.');
disp('That is, that value of beta such that the Kendall K statistic');
disp('  for  Y_i - beta X_i  vs  X_i  is 0 (or 1 or -1)');
disp('EXERCISE: If X_i<X_j for i<j, prove that this estimate of beta');
disp('  is the median of the ratios  (Y_j-Y_i)/(X_j-X_i) (for X_i < X_j).');
disp('  This estimator is also called Theil`s estimator of the slope beta.');

% Estimate beta=beta_hat = median of slopes (Y_i-Y_j)/(X_i-X_j)
% Fill array of slopes:
slopes=zeros(1,nn*nn+3);  nd=0;
for i=1:nn;  for j=i+1:nn;
  xdiff=xx(j)-xx(i);
  if xdiff~=0 nd=nd+1;  slopes(nd)=(yy(j)-yy(i))/xdiff;  end;
  end;  end;
%% DON`T COUNT unused values in slopes() before finding median
slopes=slopes(1:nd);   
beta_hat=median(slopes);
slopes=sort(slopes);

fprintf('\nThe nd=%d difference ratios (Y_j-Y_i)/(X_j-X_i)\n',nd);
fprintf('  range between %g (minimum) and %g (maximum)\n',...
  slopes(1), slopes(nd));
fprintf('The median is beta_hat=%6.4f, which is Theil`s estimate\n',beta_hat);
disp('  of the slope');

disp(' ');
disp('Given beta, mu=mu_hat is the median of the Walsh averages');
disp('  of Y_i - beta X_i.');
% Given beta, mu=mu_hat is the median of the Walsh averages
%   of Y(i)-beta_hat*X(i)
% We use the name resid(i) for Y_i - beta_hat*X_i, even though
%   they are not really residuals without the mu.
resid=zeros(1,nn);
for i=1:nn;  resid(i)=yy(i)-beta_hat*xx(i);  end;
wsums=zeros(1,nn*nn+3);  nw=0;
for i=1:nn;  for j=i:nn;
  nw=nw+1;  wsums(nw)=(resid(i)+resid(j))/2;  end;  end;
%% Again, REMOVE unused values before finding median
wsums=wsums(1:nw);
mu_hat=median(wsums);
fprintf('The median of nw=%d Walsh averages is %6.4f\n', nw,mu_hat);
fprintf('Thus Theil`s regression line in this case is\n');
fprintf('    Y = %6.4f + %6.3f\n', beta_hat,mu_hat);

disp(' ');
disp('Theil also derived a nonparametric confidence interval');
fprintf('  for beta based on quantiles of the nd=%d difference ratios\n',...
  nd);
disp('Using a large-sample estimate of the quantiles,');
% Large-sample nonparametric confidence interval for true beta
cca = floor(1.960*sqrt(varC));   %% Least integer <= 1.960*sqrt(..)
npa=nn*(nn-1)/2;
klow=floor((npa-cca)/2);
khigh=floor((npa+cca)/2) + 1;
% The slopes array is already sorted

disp ('  the nonparametric 95% Conf.Int. for beta bracketing beta_hat');
disp ('  and the corresponding offsets in the slopes array are');
fprintf('    (%5.3f(%d)   %5.3f(%g)   %5.3f(%d))\n', ...
  slopes(klow),klow, beta_hat,(npa+1)/2, slopes(khigh),khigh);

% Find MS (Mean-Square) and Sum_Absolute_Value errors
%   for Theil regression
resid_theil = yy - beta_hat*xx - mu_hat;
err_theil_sabs = sum(abs(resid_theil))/nn;
err_theil_ms = sqrt(sum(resid_theil.*resid_theil)/nn);

disp(' ');
disp('Regression line errors:');
disp('  E1 is mean of absolute values of residuals');
disp('  E2 is root-mean-square of residuals');

fprintf('  Theil:          Y = %7.4f X + %8.3f   E1=%7.2f   E2=%7.2f\n', ...
   beta_hat, mu_hat, err_theil_sabs, err_theil_ms);

% Find the corresponding least-squares regression values
% beta_lsq, mu_lsq, resid_lsq computed earlier
err_lsq_sabs=sum(abs(resid_lsq))/nn;
err_lsq_ms=sqrt(sum(resid_lsq.*resid_lsq)/nn);

fprintf('  Least-Squares:  Y = %7.4f X + %8.3f   E1=%7.2f   E2=%7.2f\n', ...
   beta_lsq, mu_lsq, err_lsq_sabs, err_lsq_ms);


disp(' ');
disp('LET`S LOOK AT LOG(Y) on X instead of Y on X:');
disp('SIMPLE REGRESSION OF log(Y) on X  (log(Y) = beta*X + mu)');

logy = log(yy);

% Use Theil`s nonparametric method
% First estimate slope (beta)
slopes=zeros(1,nn*nn+3);  nd=0;
for i=1:nn;  for j=i+1:nn;
  xdiff=xx(j)-xx(i);
  if xdiff~=0 nd=nd+1;  slopes(nd)=(logy(j)-logy(i))/xdiff;  end;
  end;  end;

% Theil`s estimate of beta for log(Y) on X
slopes=slopes(1:nd);   %% DON`T COUNT unused values in slopes()
beta_hat=median(slopes);

% Estimate mu as the median of Walsh sums of log(Y(i))-beta_hat*X(i)
% resid now stands for  log(Y) - beta_hat*X
resid=zeros(1,nn);
for i=1:nn;  resid(i)=logy(i)-beta_hat*xx(i);  end;
wsums=zeros(1,nn*nn+3);  nw=0;
for i=1:nn;  for j=i:nn;
  nw=nw+1;  wsums(nw)=(resid(i)+resid(j))/2;  end;  end;
wsums=wsums(1:nw);    % DON`T COUNT unused values
mu_hat=median(wsums);

% Find MS (Mean-Square) and Sum_Absolute_Value errors
%   for Theil regression
resid = logy - beta_hat*xx - mu_hat;
err_theil_sabs = sum(abs(resid))/nn;
err_theil_ms = sqrt(sum(resid.*resid)/nn);

disp(' ');
disp('Regression line and errors:');

fprintf('  Theil:          log(Y) = %7.4f X + %5.3f   E1=%5.3f   E2=%5.3f\n', ...
   beta_hat, mu_hat, err_theil_sabs, err_theil_ms);

% Find corresponding least-squares estimates
%
ymean=sum(logy)/nn;  xmean=sum(xx)/nn;
logyc=logy-ymean;  xxc=xx-xmean;
beta_lsq = sum(logyc.*xxc)/sum(xxc.*xxc);
mu_lsq = ymean - beta_lsq*xmean;

resid = logy - beta_lsq*xx - mu_lsq;
err_lsq_sabs=sum(abs(resid))/nn;
err_lsq_ms=sqrt(sum(resid.*resid)/nn);

fprintf('  Least-Squares:  log(Y) = %7.4f X + %5.3f   E1=%5.3f   E2=%5.3f\n', ...
   beta_lsq, mu_lsq, err_lsq_sabs, err_lsq_ms);


disp(' ');
disp('For log(Y_i) and X_i:');
disp('Theil`s test of H_0:beta=0  (same as Kendall`s test)');
disp('(NOT assuming that X(i) are increasing, so different from text)');
cc=0;
for i=1:nn;  for j=i+1:nn;
  zz = (logy(j)-logy(i))*(xx(j)-xx(i));
  if (zz>0) cc=cc+1;
  elseif zz<0 cc=cc-1;  end;  end;  end;
varC = nn*(nn-1)*(2*nn+5)/18;
zz = cc/sqrt(varC);
fprintf('  C=%d   Z=%6.4f   P=%6.4f (Assuming no ties)\n', ...
  cc, zz, normpval(zz));


disp(' ');
disp('Normal-theory test of H_0:beta=0  (Student-t test)');
ymean=sum(logy)/nn;  xmean=sum(xx)/nn;
logyc=logy-ymean;    xxc=xx-xmean;
ssxx = sum(xxc.*xxc);
beta_lsq = sum(logyc.*xxc)/ssxx;
mu_lsq = ymean - beta_lsq*xmean;
resid = logy - beta_lsq*xx - mu_lsq;
mse = sum(resid.*resid)/(nn-2);
tt = beta_lsq/sqrt(mse/ssxx);
fprintf('  beta_hat=%6.4f   T=%6.4f   P=%6.4f  (df=%d)\n', ...
  beta_lsq, tt, pvtdist(tt,nn-2), nn-2);


disp(' ');
disp('For log(Y_i) and X_i:');
disp ('Nonparametric 95% Conf.Int. for beta bracketing beta_hat and');
disp ('  the corresponding offsets in the slopes array:');
slopes=sort(slopes);
fprintf('  (%6.4f(%d)   %6.4f(%g)   %6.4f(%d))\n', ...
  slopes(klow),klow, beta_hat,(npa+1)/2, slopes(khigh),khigh);



