#delim ;
*
 Certify somersd using the auto data.
 This certification script uses the following STB additions:
 jknife (from sg34 in STB-24)
*;

* Begin certification script *;
cscript "somersd: Somers' D and Kendall's tau-a" adofile somersd;

*
 Input auto data set into memory
 and define additional variables.
*;
clear;set matsize 200;
use auto,clear;desc;
* Create manuf variable *;
gene str18 manuf=ltrim(make);
replace manuf=substr(manuf,1,index(manuf," ")-1);
replace manuf=make if(index(make," ")<=1);
lab var manuf "Manufacturer";
tab manuf,missing;
* Create test weight variable *;
gene wtest=mod(_n,5);
lab var wtest "Test weight variable";
tab wtest,missing;
* Define if and in test expressions *;
global iftest="if(mod(_n,2))";
global intest="in 33/64";
tab wtest $intest $iftest,missing;

*
 Test set 1.
 Show that replay works as expected.
*;
somersd foreign mpg weight length,tr(z);
somersd;
somersd,level(99);
somersd,ci(asci);matr list asci;
somersd,level(99) ci(asci);matr list asci;
esti clear;
rcof "noisily somersd"==301;
rcof "noisily somersd,level(99)"==301;
rcof "noisily somersd,ci(asci)"==301;
rcof "noisily somersd,level(99) ci(asci)"==301;

*
 Test set 2.
 Show that out-of-range levels fail as expected.
*;
somersd foreign mpg weight length,tr(z);
rcof "noisily somersd,level(100)"==198;
rcof "noisily somersd,level(9)"==198;
rcof "noisily somersd foreign mpg weight length,tr(z) level(100)"==198;
rcof "noisily somersd foreign mpg weight length,tr(z) level(9)"==198;

*
 Test set 3.
 Show that predict fails as expected.
*;
somersd foreign mpg weight length,tr(z);
rcof "noisily predict"==498;
rcof "noisily predict yhat"==498;
rcof "noisily predict yhat,xb"==498;

*
 Test set 4.
 Show response to tiny numbers of observations.
*;
scal tiny=1e-9;
rcof "noisily somersd mpg weight length foreign [fweight=0],taua"==2000;
somersd mpg weight length foreign in 1/1,taua;
matr b=e(b);matr vce=e(V);
matr zb=0*b;matr zvce=0*vce;
assert mreldif(b,zb)<tiny;assert mreldif(vce,zvce)<tiny;
somersd mpg weight length foreign in 1/2,taua;
matr b=e(b);matr vce=e(V);
matr zb=0*b;matr zvce=0*vce;
assert mreldif(vce,zvce)<tiny;
somersd mpg weight length foreign in 1/3,taua;

*
 Test set 5.
 Show that transformations and tdist option act as expected.
*;
capture program drop trantes;
program define trantes;
syntax [,TDist];

local tiny=1e-9;
somersd foreign mpg weight length,ci(ci0) `tdist';
matr beta0=ci0[1...,1];
local n1=rowsof(ci0);local n2=colsof(ci0);
* Fisher's z *;
somersd foreign mpg weight length,ci(ci1) tr(z) `tdist';
matr beta1=ci1[1...,1];
assert mreldif(beta1,beta0)<`tiny';
* Daniels' arcsine *;
somersd foreign mpg weight length,ci(ci1) tr(asin) `tdist';
matr beta1=ci1[1...,1];
assert mreldif(beta1,beta0)<`tiny';
* Greiner's rho *;
somersd foreign mpg weight length,ci(ci1) tr(rho) `tdist';
matr beta1=ci1[1...,1];
* Transform rho to a tau scale *;
local i1=0;
while(`i1'<rowsof(ci1)){local i1=`i1'+1;
  scal ai1=beta1[`i1',1];scal ai1=(2/_pi)*asin(ai1);
  matr beta1[`i1',1]=ai1;
};
assert mreldif(beta1,beta0)<`tiny';
* Greiner's rho (z-transformed) *;
somersd foreign mpg weight length,ci(ci1) tr(zrho) `tdist';
matr beta1=ci1[1...,1];
* Transform rho to a tau scale *;
local i1=0;
while(`i1'<rowsof(ci1)){local i1=`i1'+1;
  scal ai1=beta1[`i1',1];scal ai1=(2/_pi)*asin(ai1);
  matr beta1[`i1',1]=ai1;
};
assert mreldif(beta1,beta0)<`tiny';

end;

trantes;trantes,tdist;

*
 Test set 6.
 Show that transformations, tdist and level
 interact as expected.
*;
capture program drop trantl1;
program define trantl1;
* Call somersd and test that confidence interval is correct *;
args level;

local tiny=1e-6;

qui somersd,le(`level') ci(ci);
local tran "`e(transf)'";local tdist "`e(tdist)'";
* Make CI matrix symmetric *;
local i1=0;
while(`i1'<rowsof(ci)){local i1=`i1'+1;local i2=0;
  while(`i2'<colsof(ci)){local i2=`i2'+1;
    scal aij=ci[`i1',`i2'];
    if(("`tran'"=="z")|("`tran'"=="zrho")){
      scal aij=(log(1+aij)-log(1-aij))/2;
    };
    else if("`tran'"=="asin"){scal aij=asin(aij);};
    matr ci[`i1',`i2']=aij;
  };
};
matr hwid1=ci[1...,3]-ci[1...,1];
matr hwid2=ci[1...,1]-ci[1...,2];
* Check that CI is symmetric *;
assert mreldif(hwid1,hwid2)<`tiny';
* Check that CI has appropriate half width *;
matr vari=(vecdiag(e(V)))';matr mult=hwid1;
local i1=0;
while(`i1'<rowsof(vari)){local i1=`i1'+1;
  scal a=vari[`i1',1];scal a=sqrt(a);
  scal b=hwid1[`i1',1];scal a=b/a;
  matr mult[`i1',1]=a;
};
if("`tdist'"!=""){
  * t-distribution *;
  scal mulex=invt(e(df_r),`level'/100);
};
else{
  * Normal distribution *;
  scal mulex=invnorm((`level'+100)/200);
};
matr multexp=mulex*J(rowsof(mult),1,1);
assert mreldif(mult,multexp)<`tiny';

end;

capture program drop trantl2;
program define trantl2;
* Call trantl1 under several different settings *;
syntax [,TRansf(passthru)];

somersd foreign mpg weight length,`transf';
for X in num 90/99,nohe:trantl1 X;
somersd foreign mpg weight length,`transf' tdist;
for X in num 90/99,nohe:trantl1 X;
somersd foreign mpg weight length,taua `transf';
for X in num 90/99,nohe:trantl1 X;
somersd foreign mpg weight length,taua `transf' tdist;
for X in num 90/99,nohe:trantl1 X;

end;

* Call trantl2 under several transformations *;
trantl2;
trantl2,tr(z);
trantl2,tr(asin);
trantl2,tr(rho);
trantl2,tr(zrho);

*
 Test set 7.
 Show that somersd,taua gives same results as ktau.
*;

capture program drop kttest;
program define kttest;
syntax [in] [if];

scal tiny=1e-14;
somersd mpg weight length foreign `if' `in',taua;
matr beta=e(b);matr list beta;
ktau mpg mpg `if' `in';matr taua=beta[1,"mpg"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;
ktau mpg length `if' `in';matr taua=beta[1,"length"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;
ktau mpg weight `if' `in';matr taua=beta[1,"weight"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;
ktau mpg foreign `if' `in';matr taua=beta[1,"foreign"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;

end;

kttest;
kttest $iftest;
kttest $intest;
kttest $iftest $intest;

*
 Test set 8.
 Show that somersd,taua tdist gives same result
 as jknife and ktau.
 (Except that jknife and ktau takes ages,
 because it requires an amount of processing time
 that varies cubically with the number of observations,
 whereas somersd is streamlined, requiring processing time
 varying quadratically with the number of distinct values of varlist.)
*;
capture program drop jktes1;
program define jktes1;
syntax [varname] [in];

scal tiny=1e-6;

jknife clear;
jknife cmd ktau mpg `varlist' `in';
jknife stat 1 taua=r(tau_a);
jknife n=r(N);
jknife do;
drop taua;
matr mid=cimat["`varlist'",1];
matr min=cimat["`varlist'",2];matr max=cimat["`varlist'",3];
scal tau=trace(mid);
scal taumin=trace(min);scal taumax=trace(max);
assert reldif(tau,$S_2)<tiny;
assert reldif(taumin,$S_5)<tiny;assert reldif(taumax,$S_6)<tiny;
 
end;

capture program drop jktes2;
program define jktes2;
* Call jktes1 under various settings *;
syntax [in];

somersd mpg weight length foreign `in',taua tdist ci(cimat);
matr list cimat;
jktes1 mpg `in';
jktes1 weight `in';
jktes1 length `in';
jktes1 foreign `in';

end;

jktes2;
jktes2 $intest;

*
 Test set 9.
 Show that Kendall's tau-b(X,Y) (provided by tabulate)
 is the geometric mean of Somers' D(X,Y) and D(Y,X).
*;
capture program drop tbtes1;
program define tbtes1;
syntax [if] [in] [fweight];

scal tiny=1e-9;

tabulate foreign rep78 `if' `in' [fwei `exp'],taub;
scal taub=r(taub);
somersd foreign rep78 `if' `in' [fwei `exp'];
matr beta=e(b);matr tauxy=beta[1,"rep78"];
scal tauxy=trace(tauxy);
somersd rep78 foreign `if' `in' [fwei `exp'];
matr beta=e(b);matr tauyx=beta[1,"foreign"];
scal tauyx=trace(tauyx);
scal taubxy=sqrt(tauxy*tauyx);
disp "Tau-b (from somersd) = " taubxy;
assert reldif(taubxy,taub)<tiny;

end;

tbtes1 [fwei=1];
qui tbtes1 $iftest [fwei=1];
qui tbtes1 $intest [fwei=1];
qui tbtes1 $iftest $intest [fwei=1];
tbtes1 [fwei=wtest];
qui tbtes1 $iftest [fwei=wtest];
qui tbtes1 $intest [fwei=wtest];
qui tbtes1 $iftest $intest [fwei=wtest];

*
 Test set 10.
 Show that Somers' D for binary y-variables
 is the difference between proportions
 (with similar variance)
*;
scal tiny=1e-14;scal small=1e-2;
*
 Create binary variables to compare
 between domestic and foreign cars
*;
gene byte dear=price>5006.5;
gene byte heavy=weight>3190;
gene byte bulky=displ>196;
gene byte guzzler=mpg<20;
summ dear heavy bulky guzzler;

capture program drop pdtes1;
program define pdtes1;
syntax [if] [in] [fweight];
if("`weight'"==""){local weight="fweight";};
if("`exp'"==""){local exp="= 1";};

* Create matrices of Somers' D estimates and variances *;
somersd foreign dear heavy bulky guzzler `if' `in' [fwei`exp'];
matr beta=e(b);matr vari=vecdiag(e(V));
*
 Do regressions to find differences between proportions
 and check that differences are equal
 to corresponding Somers' D estimates
 and that their variances are similar
 to corresponding Somers' D variances
*;
regr dear foreign `if' `in' [fwei`exp'],robust;
matr breg=e(b);matr vreg=e(V);
matr bs=beta[1,"dear"];matr br=breg[1,"foreign"];
assert mreldif(bs,br)<tiny;
matr vs=vari[1,"dear"];matr vr=vreg["foreign","foreign"];
assert mreldif(vs,vr)<small;
regr heavy foreign `if' `in' [fwei`exp'],robust;
matr breg=e(b);matr vreg=e(V);
matr bs=beta[1,"heavy"];matr br=breg[1,"foreign"];
assert mreldif(bs,br)<tiny;
matr vs=vari[1,"heavy"];matr vr=vreg["foreign","foreign"];
assert mreldif(vs,vr)<small;
regr bulky foreign `if' `in' [fwei`exp'],robust;
matr breg=e(b);matr vreg=e(V);
matr bs=beta[1,"bulky"];matr br=breg[1,"foreign"];
assert mreldif(bs,br)<tiny;
matr vs=vari[1,"bulky"];matr vr=vreg["foreign","foreign"];
assert mreldif(vs,vr)<small;
regr guzzler foreign `if' `in' [fwei`exp'],robust;
matr breg=e(b);matr vreg=e(V);
matr bs=beta[1,"guzzler"];matr br=breg[1,"foreign"];
assert mreldif(bs,br)<tiny;
matr vs=vari[1,"guzzler"];matr vr=vreg["foreign","foreign"];
assert mreldif(vs,vr)<small;

end;

pdtes1;
qui pdtes1 $iftest;
qui pdtes1 $intest;
qui pdtes1 $iftest $intest;
pdtes1 [fwei=wtest];
qui pdtes1 $iftest [fwei=wtest];
qui pdtes1 $intest [fwei=wtest];
qui pdtes1 $iftest $intest [fwei=wtest];

drop dear heavy bulky guzzler;

*
 Test set 11.
 Show that in and if work as they should,
 and produce the same results as keeping if and in
 and then running the command without if or in.
*;
capture program drop iites1;
program define iites1;
syntax [if] [in] [fweight iweight pweight]
 [,TAua CLuster(passthru) TRansf(passthru)];
if("`weight'"==""){local weight="fweight";};
if("`exp'"==""){local exp="= 1";};

scal tiny=1e-14;

somersd foreign mpg weight length `if' `in' [`weight'`exp']
 ,`taua' `cluster' `transf';
matr def b0=e(b);matr def v0=e(V);scal n0=e(N);scal nc0=e(N_clust);
preserve;capture keep `if' `in';
qui somersd foreign mpg weight length [`weight'`exp'],
 `taua' `cluster' `transf';
matr def b1=e(b);matr def v1=e(V);scal n1=e(N);scal nc1=e(N_clust);
assert mreldif(b1,b0)<tiny;assert mreldif(v1,v0)<tiny;
assert reldif(n0,n1)<tiny;assert reldif(nc0,nc1)<tiny;
restore;

end;

capture program drop iites2;
program define iites2;
* Call iiwtes1 with a number of settings *;
syntax[,TAua CLuster(passthru) TRansf(passthru)];

iites1,`taua' `cluster' `transf';
iites1 $intest,`taua' `cluster' `transf';
iites1 $iftest,`taua' `cluster' `transf';
iites1 $iftest $intest,`taua' `cluster' `transf';
qui iites1 [fwei=wtest],`taua' `cluster' `transf';
qui iites1 $intest [fwei=wtest],`taua' `cluster' `transf';
qui iites1 $iftest [fwei=wtest],`taua' `cluster' `transf';
qui iites1 $iftest $intest [fwei=wtest],`taua' `cluster' `transf';
qui iites1 [iwei=wtest],`taua' `cluster' `transf';
qui iites1 $intest [iwei=wtest],`taua' `cluster' `transf';
qui iites1 $iftest [iwei=wtest],`taua' `cluster' `transf';
qui iites1 $iftest $intest [iwei=wtest],`taua' `cluster' `transf';
qui iites1 [pwei=wtest],`taua' `cluster' `transf';
qui iites1 $intest [pwei=wtest],`taua' `cluster' `transf';
qui iites1 $iftest [pwei=wtest],`taua' `cluster' `transf';
qui iites1 $iftest $intest [pwei=wtest],`taua' `cluster' `transf';

end;

capture program drop iites3;
program define iites3;
* Call iites2 with a number of settings *;
syntax[,TRansf(passthru)];

iites2,`transf';
iites2,`transf' cluster(manuf);
iites2,`transf' taua;
iites2,`transf' taua cluster(manuf);

end;

iites3;
qui iites3,tr(iden);
qui iites3,tr(z);
qui iites3,tr(asin);
qui iites3,tr(rho);
qui iites3,tr(zrho);

*
 Test set 12.
 Show that fweights (and iweights and pweights in the clustered case)
 work as they should
 and produce the same results as expanding
 and then running the command without fweights.
*;
capture program drop fwtes1;
program define fwtes1;
syntax [if] [in] [,TAua CLuster(passthru) TRansf(passthru)];

preserve;
capture keep `if' `in';keep if(wtest>0);
expand wtest;
disp "Expanded data, without weights";
somersd foreign mpg weight length,
 `taua' `cluster' `transf';
matr def b0=e(b);matr def v0=e(V);scal n0=e(N);scal nc0=e(N_clust);
restore;
disp "Unexpanded data, with fweights";
somersd foreign mpg weight length `if' `in' [fwei=wtest]
 ,`taua' `cluster' `transf';
matr def b1=e(b);matr def v1=e(V);scal n1=e(N);scal nc1=e(N_clust);
assert mreldif(b1,b0)<tiny;assert mreldif(v1,v0)<tiny;
assert reldif(n0,n1)<tiny;assert reldif(nc0,nc1)<tiny;
if("`cluster'"!=""){
  *
   For clustered data,
   iweights and pweights should act as fweights
   (except that e(N) is lower)
  *;
  disp "Unexpanded data, with iweights";
  somersd foreign mpg weight length `if' `in' [iwei=wtest]
   ,`taua' `cluster' `transf';
  matr def b1=e(b);matr def v1=e(V);scal n1=e(N);scal nc1=e(N_clust);
  assert mreldif(b1,b0)<tiny;assert mreldif(v1,v0)<tiny;
  assert reldif(nc0,nc1)<tiny;
  disp "Unexpanded data, with pweights";
  somersd foreign mpg weight length `if' `in' [pwei=wtest]
   ,`taua' `cluster' `transf';
  matr def b1=e(b);matr def v1=e(V);scal n1=e(N);scal nc1=e(N_clust);
  assert mreldif(b1,b0)<tiny;assert mreldif(v1,v0)<tiny;
  assert reldif(nc0,nc1)<tiny;
};

end;

capture program drop fwtes2;
program define fwtes2;
* Call fwtes1 with a number of settings *;
syntax[,TAua CLuster(passthru) TRansf(passthru)];

fwtes1,`taua' `cluster' `transf';
qui fwtes1 $intest,`taua' `cluster' `transf';
qui fwtes1 $iftest,`taua' `cluster' `transf';
qui fwtes1 $iftest $intest,`taua' `cluster' `transf';

end;

capture program drop fwtes3;
program define fwtes3;
* Call fwtes2 with a number of settings *;
syntax[,TRansf(passthru)];

fwtes2,`transf';
fwtes2,`transf' cluster(manuf);
fwtes2,`transf' taua;
fwtes2,`transf' taua cluster(manuf);

end;

fwtes3;
qui fwtes3,tr(iden);
qui fwtes3,tr(z);
qui fwtes3,tr(asin);
qui fwtes3,tr(rho);
qui fwtes3,tr(zrho);

*
 Test set 13.
 Show that uniformly scaling the weights
 does not affect variances in the presence of clustering
 and, in the absence of clustering,
 scales variances down in the case of fweights
 and has no effect in the case of iweights and pweights.
*;
capture program drop wttes1;
program define wttes1;
syntax [if] [in] [fweight iweight pweight /] [,TAua TRansf(passthru)];
if("`weight'"==""){local weight="fweight";};
if("`exp'"==""){local exp="1";};

scal tiny=1e-9;scal small=1/3;
* Unclustered *;
somersd mpg weight length foreign `if' `in' [`weight'=`exp']
  ,`taua' `transf';
scal nobs0=e(N);matr b0=e(b);matr vce0=e(V);
local i1max=4;
local i1=0;
while(`i1'<`i1max'){local i1=`i1'+1;
  local wsc=4^(`i1');
  if("`weight'"=="fweight"){
    disp "Unclustered, `weight's (`exp') scaled by `wsc'";
    somersd mpg weight length foreign `if' `in' [`weight'=(`exp')*`wsc']
     ,`taua' `transf';
    *
     For unclustered fweights,
     scale nobs and variance appropriately by weight before matching
     and match variance-covariance matrix very approximately
     as the jackknife method produces estimates and variance estimates
     which do not scale exactly with fweights even in theory,
     and may only be in the right binary order of magnitude
    *;
    scal nobs=e(N);matr b=e(b);matr vce=`wsc'*e(V);
    assert reldif(nobs,nobs0*`wsc')<tiny;
    assert mreldif(b,b0)<small;assert mreldif(vce,vce0)<small;
  };
  else{
    disp "Unclustered, `weight's (`exp') scaled by `wsc'";
    somersd mpg weight length foreign `if' `in' [`weight'=(`exp')]
      ,`taua' `transf';
    scal nobs=e(N);matr b=e(b);matr vce=e(V);
    assert reldif(nobs,nobs0)<tiny;
    assert mreldif(b,b0)<tiny;assert mreldif(vce,vce0)<tiny;
  };
};
* Clustered *;
somersd mpg weight length foreign `if' `in' [`weight'=`exp']
  ,`taua' `transf' cluster(manuf);
scal nobs0=e(N);matr b0=e(b);matr vce0=e(V);
local i1max=4;
local i1=0;
while(`i1'<`i1max'){local i1=`i1'+1;
  local wsc=4^(`i1'-1);
  if("`weight'"=="fweight"){
    disp "Clustered, `weight's (`exp') scaled by `wsc'";
    somersd mpg weight length foreign `if' `in' [`weight'=(`exp')*`wsc']
      ,`taua' `transf' cluster(manuf);
    scal nobs=e(N);matr b=e(b);matr vce=e(V);
    *
     For clustered fweights,
     scale nobs appropriately by weight before matching
    *;
    assert reldif(nobs,nobs0*`wsc')<tiny;
    assert mreldif(b,b0)<tiny;assert mreldif(vce,vce0)<tiny;
  };
  else{
    disp "Clustered, `weight's (`exp') scaled by `wsc'";
    somersd mpg weight length foreign `if' `in' [`weight'=(`exp')*`wsc']
      ,`taua' `transf' cluster(manuf);
    scal nobs=e(N);matr b=e(b);matr vce=e(V);
    assert reldif(nobs,nobs0)<tiny;
    assert mreldif(b,b0)<tiny;assert mreldif(vce,vce0)<tiny;
  };
};

end;

capture program drop wttes2;
program define wttes2;
* Call wttes1 with a number of weights *;
syntax [if] [in] [,TAua TRansf(passthru)];

wttes1 `if' `in' [fwei=1],`taua' `transf';
wttes1 `if' `in' [iwei=1],`taua' `transf';
wttes1 `if' `in' [pwei=1];
wttes1 `if' `in' [fwei=wtest],`taua' `transf';
wttes1 `if' `in' [iwei=wtest],`taua' `transf';
wttes1 `if' `in' [pwei=wtest],`taua' `transf';

end;

capture program drop wttes3;
program define wttes3;
* Call wttes2 with a number of settings *;
syntax [,TAua TRansf(passthru)];

wttes2,`taua' `transf';
qui wttes2 $iftest,`taua' `transf';
qui wttes2 $intest,`taua' `transf';
qui wttes2 $iftest $intest,`taua' `transf';

end;

wttes3;
qui wttes3,tr(iden);
qui wttes3,tr(z);
qui wttes3,tr(asin);
qui wttes3,tr(rho);
qui wttes3,tr(zrho);
wttes3,taua;
qui wttes3,taua tr(iden);
qui wttes3,taua tr(z);
qui wttes3,taua tr(asin);
qui wttes3,taua tr(rho);
qui wttes3,taua tr(zrho);

*
 Test set 14.
 Show that somersd works with enormous varlists.
*;
somersd
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 ,tr(z);
somersd
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 ,tr(z) taua;
somersd
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 ,tr(z);
somersd
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 foreign price mpg rep78 hdroom trunk weight length turn displ gratio
 ,tr(z) taua;

exit;