#delim ;
version 6.0;
*
 Certify somersd using the Auto data
*;

* Open log file *;log using cert1.log,replace;
capture noisily{
* Beginning of main capture block *;

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

* Use auto data set and define additional variables *;
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);
tab manuf,missing;

* Show that replay works as expected *;
somersd foreign mpg weight length;
somersd;
somersd,level(99);
esti clear;
rcof "noisily somersd"==301;
rcof "noisily somersd,level(99)"==301;

*
 Show that Somers' D for binary y-variables
 is the difference between proportions
 (with similar variance)
*;
scal tiny=1e-14;scal small=1e-3;
*
 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;
* Create matrices of Somers' D estimates and variances *;
somersd foreign dear heavy bulky guzzler;
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,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,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,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,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;
drop dear heavy bulky guzzler;

* Show that Kendall's tau-a gives the same value as ktau *;
scal tiny=1e-14;
somersd mpg weight length foreign,taua;
matr beta=e(b);matr list beta;
ktau mpg mpg;matr taua=beta[1,"mpg"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;
ktau mpg length;matr taua=beta[1,"length"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;
ktau mpg weight;matr taua=beta[1,"weight"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;
ktau mpg foreign;matr taua=beta[1,"foreign"];
scal taua=trace(taua);assert reldif(taua,r(tau_a))<tiny;

*
 Show that transformed Kendall's tau-a (after inverse transform)
 gives the same value as ktau
*;
* Fisher's z transform *;
scal tiny=1e-14;
somersd mpg weight length foreign,taua tr(z);
matr beta=e(b);matr list beta;
ktau mpg mpg;matr taua=beta[1,"mpg"];
scal taua=trace(taua);scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg length;matr taua=beta[1,"length"];
scal taua=trace(taua);scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg weight;matr taua=beta[1,"weight"];
scal taua=trace(taua);scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg foreign;matr taua=beta[1,"foreign"];
scal taua=trace(taua);scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
assert reldif(taua,r(tau_a))<tiny;
* Daniels' arcsine transform *;
scal tiny=1e-9;
somersd mpg weight length foreign,taua tr(asin);
matr beta=e(b);matr list beta;
ktau mpg mpg;matr taua=beta[1,"mpg"];
scal taua=trace(taua);scal taua=sin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg length;matr taua=beta[1,"length"];
scal taua=trace(taua);scal taua=sin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg weight;matr taua=beta[1,"weight"];
scal taua=trace(taua);scal taua=sin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg foreign;matr taua=beta[1,"foreign"];
scal taua=trace(taua);scal taua=sin(taua);
assert reldif(taua,r(tau_a))<tiny;
* Greiner's rho transform *;
scal tiny=1e-9;
somersd mpg weight length foreign,taua tr(rho);
matr beta=e(b);matr list beta;
ktau mpg mpg;matr taua=beta[1,"mpg"];
scal taua=trace(taua);scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg length;matr taua=beta[1,"length"];
scal taua=trace(taua);scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg weight;matr taua=beta[1,"weight"];
scal taua=trace(taua);scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg foreign;matr taua=beta[1,"foreign"];
scal taua=trace(taua);scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;
* Fisher's z transform of Greiner's rho transform *;
scal tiny=1e-9;
somersd mpg weight length foreign,taua tr(zrho);
matr beta=e(b);matr list beta;
ktau mpg mpg;matr taua=beta[1,"mpg"];
scal taua=trace(taua);
scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg length;matr taua=beta[1,"length"];
scal taua=trace(taua);
scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg weight;matr taua=beta[1,"weight"];
scal taua=trace(taua);
scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;
ktau mpg foreign;matr taua=beta[1,"foreign"];
scal taua=trace(taua);
scal taua=exp(2*taua);scal taua=(taua-1)/(taua+1);
scal taua=(2/_pi)*asin(taua);
assert reldif(taua,r(tau_a))<tiny;

*
 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;

*
 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 if iweights and pweights
*;
capture program drop wtcert;
program define wtcert;
syntax[,TAua TDist];

scal tiny=1e-9;scal small=0.5e-2;
* Unclustered *;
somersd mpg weight length foreign,`taua' `tdist' tr(z);
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);
  disp "Unclustered, iweights scaled by `wsc'";
  somersd mpg weight length foreign [iwei=`wsc'],`taua' `tdist' tr(z);
  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;
  disp "Unclustered, pweights scaled by `wsc'";
  somersd mpg weight length foreign [pwei=`wsc'],`taua' `tdist' tr(z);
  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;
  disp "Unclustered, fweights scaled by `wsc'";
  somersd mpg weight length foreign [fwei=`wsc'],`taua' `tdist' tr(z);
  *
   For unclustered fweights,
   scale nobs and variance appropriately by weight before matching
   and match variance-covariance matrix very approximately
   as the jackknife method produces variance estimates
   which do not scale exactly with fweights even in theory
  *;
  scal nobs=e(N);matr b=e(b);matr vce=`wsc'*e(V);
  assert reldif(nobs,nobs0*`wsc')<tiny;
  assert mreldif(b,b0)<tiny;assert mreldif(vce,vce0)<small;
  matr list vce0;matr list vce;
};
* Clustered *;
somersd mpg weight length foreign,`taua' `tdist' tr(z) 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);
  disp "Unclustered, iweights scaled by `wsc'";
  somersd mpg weight length foreign [iwei=`wsc'],`taua' `tdist' tr(z) 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;
  disp "Unclustered, pweights scaled by `wsc'";
  somersd mpg weight length foreign [pwei=`wsc'],`taua' `tdist' tr(z) 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;
  disp "Unclustered, fweights scaled by `wsc'";
  somersd mpg weight length foreign [fwei=`wsc'],`taua' `tdist' tr(z) 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;
};

end;

wtcert;
wtcert,taua;
wtcert,tdist;
wtcert,taua tdist;


* End of main capture block *;
};
* Close log file and exit *;log close;exit;