#delim ;
*
 Certification script for -multproc- and -smileplot-
*;

* Begin certification script *;
cscript "Certification of multproc and smileplot"
 adofiles multproc smileplot;

*
 Input the sequence of P-values and calculated cut-offs from:
 Benjamini Y, Krieger A, Yekutieli D.
 Two-staged linear step up FDR controlling procedure.
 Downloadable from http://www.math.tau.ac.il/~ybenja/.
*;
clear;
#delim cr;
input prank p cutoff1 cutoff2
 1 .006 .005 .00625
 2 .009 .010 .01250
 3 .018 .015 .01875
 4 .028 .020 .02500
 5 .030 .025 .03125
 6 .045 .030 .03750
 7 .150 .035 .04375
 8 .430 .040 .05000
 9 .625 .045 .05625
10 .810 .050 .06250
end
#delim ;
lab var prank "Rank of p-value";
lab var p "p-value";
lab var cutoff1 "Cutoff First Stage";
lab var cutoff2 "Cutoff Second Stage";
format p cutoff1 %8.3f;format cutoff2 %8.5f;
desc;
list;

* Define fantasy estimate for smile plots *;
gene estimate=invnorm(1-p/2)*(2*mod(_n,2)-1);

* All methods *;

* Error capture *;
cap noi multproc,me(bogus);
assert _rc==498;
cap noi multproc,punc(-1);
assert _rc==498;
cap noi multproc,punc(1.5);
assert _rc==498;
cap noi multproc,pc(-1);
assert _rc==498;
cap noi multproc,pc(1.5);
assert _rc==498;

multproc,pc(0.005);
for M in any bonferroni sidak holm holland liu1 liu2 hochberg rom simes yekutieli krieger:
 smileplot,me(M) t1("Smile plot for method(M)") \ more;

* P-plots and trivial validations *;
for M in any  bonferroni sidak holm holland liu1 liu2 hochberg rom simes yekutieli krieger:
 multproc,me(M) ra(rank) cr(pc_M) nh(nhcred) rej(nhrej) \
 assert nhrej==!nhcred \
 graph7 p pc_M rank,s(O.) c(.L) xlab ylab t1("P-plot for method(M)") \ more \
 drop rank nhcred nhrej;
assert pc_bonferroni<pc_sidak;
assert pc_bonferroni<pc_holm if prank>1;
assert pc_bonferroni==pc_holm if prank==1;
assert pc_sidak<pc_holland if prank>1;
assert pc_sidak==pc_holland if prank==1;
assert pc_yekutieli<pc_simes;
assert pc_krieger>pc_simes;
for M in any bonferroni sidak:summ pc_M \ assert r(min)==r(max);
for M in any holm holland liu1 liu2 hochberg rom simes yekutieli krieger:
 summ pc_M \ assert 0<r(min) \ assert r(min)<r(max) \ assert r(max)<1;
drop pc_*;

* Show that -float- option works *;
multproc,me(simes) cr(cr1) ra(ra1) nh(nh1) rej(rj1) gpu(punc1) gpc(pcor1);
savedresults save resf r();
multproc,me(simes) cr(cr2) ra(ra2) nh(nh2) rej(rj2) gpu(punc2) gpc(pcor2) float;
savedresults compare resf r(),verbose;
foreach X in cr ra nh rj punc pcor {;
  assert reldif(`X'1,`X'2)<1e-6;
};
desc cr? ra? nh? rj? punc? pcor?;
list p cr? ra? nh? rj? punc? pcor?,nod;
savedresults drop resf;
drop cr? ra? nh? rj? punc? pcor?;


* Show that -fast- option works *;
multproc,me(simes) cr(cr1) ra(ra1) nh(nh1) rej(rj1) gpu(punc1) gpc(pcor1) fast;
savedresults save resf r();
multproc,me(simes) cr(cr2) ra(ra2) nh(nh2) rej(rj2) gpu(punc2) gpc(pcor2);
savedresults compare resf r(),verbose;
foreach X in cr ra nh rj punc pcor {;
  assert `X'1==`X'2;
};
desc cr? ra? nh? rj? punc? pcor?;
list p cr? ra? nh? rj? punc? pcor?,nod;
savedresults drop resf;
drop cr? ra? nh? rj? punc? pcor?;

* Simes method (Stage 1) *;
multproc,me(simes) cr(cu1) ra(ra1) nh(nh1) rej(rj1);
format cu1 %8.3f;
assert ra1==prank;
gene rd1=reldif(cutoff1,cu1);
assert rd1<=1e-6;
count if !nh1;assert r(N)==2;
list prank p ra1 cutoff1 cu1 rd1 nh1 rj1;
graph7 cu1 cutoff1;more;

* Krieger method (Stage 2) *;
local q=0.05/(1-0.05);
multproc,me(krieger) punc(`q') cr(cu2) ra(ra2) nh(nh2) rej(rj2);
format cu2 %8.5f;
assert ra2==prank;
gene rd2=reldif(cutoff2,cu2);
assert rd2<=1e-6;
count if !nh2;assert r(N)==5;
list prank p ra2 cutoff2 cu2 rd2 nh2 rj2;
graph7 cu2 cutoff2;more;

/*
 Input the sequence of calculated cut-offs from:
 Rom DM.
 A sequentially rejective test procedure based on a modified Bonferroni inequality.
 Biometrika 1990; 77(3): 663-665.
*/
#delim cr;
input mh0_05 h0_05 mh0_01 h0_01
5.11e-3 5.00e-3 1.01e-3 1.00e-3
5.68e-3 5.56e-3 1.12e-3 1.11e-3
6.39e-3 6.25e-3 1.26e-3 1.25e-3
7.30e-3 7.14e-3 1.43e-3 1.43e-3
8.51e-3 8.33e-3 1.67e-3 1.67e-3
1.02e-2 1.00e-2 2.01e-3 2.00e-3
1.27e-2 1.25e-2 2.51e-3 2.50e-3
1.69e-2 1.67e-2 3.34e-3 3.33e-3
2.50e-2 2.50e-2 5.00e-3 5.00e-3
5.00e-2 5.00e-2 1.00e-2 1.00e-2
#delim ;
lab var mh0_05 "Rom cutoff (alpha=0.05)";
lab var  h0_05 "Hochberg cutoff (alpha=0.05)";
lab var mh0_01 "Rom cutoff (alpha=0.01)";
lab var  h0_01 "Hochberg cutoff (alpha=0.01)";
format mh0_05 h0_05 mh0_01 h0_01 %9.2e;
desc;

*
 Hochberg method
*;
* Alpha=0.05 *;
multproc,me(hochberg) punc(0.05) cr(cu3) ra(ra3) nh(nh3) rej(rj3);
format cu3 %9.2e;
assert ra3==prank;
gene rd3=reldif(h0_05,cu3);
assert rd3<=1e-4;
assert nh3==1-rj3;
list prank ra3 h0_05 cu3 rd3 nh3 rj3;
graph7 cu3 h0_05;more;
* Alpha=0.01 *;
multproc,me(hochberg) punc(0.01) cr(cu4) ra(ra4) nh(nh4) rej(rj4);
format cu4 %9.2e;
assert ra4==prank;
gene rd4=reldif(h0_01,cu4);
assert rd4<=1e-4;
assert nh4==1-rj4;
list prank ra4 h0_01 cu4 rd4 nh4 rj4;
graph7 cu4 h0_01;more;

*
 Rom method
*;
* Alpha=0.05 *;
multproc,me(rom) punc(0.05) cr(cu5) ra(ra5) nh(nh5) rej(rj5);
format cu5 %9.2e;
assert ra5==prank;
gene rd5=reldif(mh0_05,cu5);
assert rd5<=1e-4;
assert nh5==1-rj5;
list prank ra5 mh0_05 cu5 rd5 nh5 rj5;
graph7 cu5 mh0_05;more;
* Alpha=0.01 *;
multproc,me(rom) punc(0.01) cr(cu6) ra(ra6) nh(nh6) rej(rj6);
format cu6 %9.2e;
assert ra6==prank;
gene rd6=reldif(mh0_01,cu6);
assert rd6<=1e-4;
assert nh6==1-rj6;
list prank ra6 mh0_01 cu6 rd6 nh6 rj6;
graph7 cu6 mh0_01;more;

*
 Holm method
 (this should give the same individual critical P-values as the Hochberg method,
 although the overall critical P-value may be different
 because the Holm procedure is step-down
 and the Hochberg procedure is step-up)
*;
* Alpha=0.05 *;
multproc,me(holm) punc(0.05) cr(cu7) ra(ra7) nh(nh7) rej(rj7);
format cu7 %9.2e;
assert ra7==prank;
gene rd7=reldif(h0_05,cu7);
assert rd7<=1e-4;
assert nh7==1-rj7;
list prank ra7 h0_05 cu7 rd7 nh7 rj7;
graph7 cu7 h0_05;more;
* Alpha=0.01 *;
multproc,me(holm) punc(0.01) cr(cu8) ra(ra8) nh(nh8) rej(rj8);
format cu8 %9.2e;
assert ra8==prank;
gene rd8=reldif(h0_01,cu8);
assert rd8<=1e-4;
assert nh8==1-rj8;
list prank ra8 h0_01 cu8 rd8 nh8 rj8;
graph7 cu8 h0_01;more;

/*
 Input the series of P-values and calculated cut-offs from:
 Benjamini Y, Liu W.
 A distribution-free multiple-test procedure that controls the false discovery rate.
 Report, Dept. of Statistics and OR, Tel Aviv University, RP-SOR-99-3.
 Downloadable from Yoav Benjamini's website at http://www.math.tau.ac.il/~ybenja/.
 The P-values belong to comparisons of mean tensile strength
 between 4 mixing techniques for Portland cement.
*/
clear;
#delim cr;
input pval pcliu2 pcyeku
0.000 0.008 0.003
0.003 0.012 0.007
0.006 0.019 0.010
0.017 0.033 0.014
0.039 0.075 0.017
0.651 0.300 0.020
end
#delim ;
gene byte compseq=_n;order compseq;
lab def compseq 1 "1 vs 4" 2 "1 vs 3" 3 "1 vs 2" 4 "2 vs 4" 5 "3 vs 4" 6 "2 vs 3";
lab val compseq compseq;
lab var compseq "Comparison";
lab var pval "P-value";
lab var pcliu2 "Liu (1999b) cutoff";
lab var pcyeku "Yekutieli cutoff";
format pval pcliu pcyeku %8.3f;
desc;
list;

* Liu (1999b) procedure *;
multproc,pv(pval) me(liu2) cr(cu1) ra(ra1) nh(nh1) rej(rj1);
format cu1 %8.3f;
gene rd1=reldif(pcliu2,cu1);
assert ra1==compseq;
assert rd1<=1e-3;
count if nh1;
assert r(N)==1;
assert rj1==1-nh1;
list compseq ra1 pval pcliu2 cu1 rd1 nh1 rj1;
graph7 cu1 pcliu2;more;

* Yekutieli (2001) procedure *;
multproc,pv(pval) me(yekutieli) cr(cu2) ra(ra2) nh(nh2) rej(rj2);
format cu2 %8.3f;
gene rd2=reldif(pcyeku,cu2);
assert ra2==compseq;
assert rd2<=1e-3;
count if nh2;
assert r(N)==3;
assert rj2==1-nh2;
list compseq ra1 pval pcyeku cu2 rd2 nh2 rj2;
graph7 cu2 pcyeku;more;

/*
 Input the series of P-values and cutoffs from:
 Benjamini Y, Liu W.
 A step-down multiple hypothesis testing procedure
 that controls the false discovery rate under independence.
 Journal of Statistical Planning and Inference 1999; 82: 163-170.
 The 4 P-values belong to differences
 in proportions of patients free of breast cancer
 between patients treated with progesterones PF and PFT
 in 4 groups of patients defined by age and progesterone receptor level.
*/
clear;
#delim cr;
input agegp prgp pval pcliu1_10 pcliu1_05
1 1 0.0362 0.0466 0.0227
2 1 0.0972 0.1056 0.0513
1 2 0.4440 0.4000 0.2000
2 2 0.0058 0.0260 0.0127
end
#delim ;
format pval pcliu1_* %8.4f;
lab def agegp 1 "<50" 2 ">50";
lab def prgp 1 "<10" 2 ">10";
lab val agegp agegp;lab val prgp prgp;
lab var pval "P-value";
lab var pcliu1_10 "Liu (1999a) 10% cutoff";
lab var pcliu1_05 "Liu (1999a) 5% cutoff";
desc;
list;

* Liu (1999a) procedure at 10% level *;
multproc,pv(pval) me(liu1) cr(cu1) ra(ra1) nh(nh1) rej(rj1) punc(0.10);
format cu1 %8.4f;
gene rd1=reldif(pcliu1_10,cu1);
assert rd1<=1e-4;
count if nh1;
assert r(N)==1;
assert nh1==1-rj1;
list agegp prgp ra1 pval pcliu1_10 cu1 rd1 nh1 rj1;
graph7 cu1 pcliu1_10;more;

* Liu (1999a) procedure at 5% level *;
multproc,pv(pval) me(liu1) cr(cu2) ra(ra2) nh(nh2) rej(rj2) punc(0.05);
format cu2 %8.4f;
gene rd2=reldif(pcliu1_05,cu2);
assert rd2<=1e-4;
count if nh2;
assert r(N)==3;
assert nh2==1-rj2;
list agegp prgp ra2 pval pcliu1_05 cu2 rd2 nh2 rj2;
graph7 cu2 pcliu1_05;more;

/*
 Demonstrate the argument of:
 Holm, S.
 A Simple Sequentially Rejective Multiple Test Procedure.
 Scandinavian Journal of Statistics 1979; 6: 65-70.
*/
clear;
#delim cr;
input refdunt seqbont
1.70 1.70
1.99 2.04
2.15 2.23
2.25 2.36
2.33 2.46
2.40 2.54
2.45 2.60
2.50 2.66
2.54 2.71
end
#delim ;
gsort -refdunt;
gene refdunp=ttail(30,refdunt)*2;
gene seqbonp=ttail(30,seqbont)*2;
lab var refdunt "Refined Dunnett t";
lab var refdunp "Refined Dunnett P";
lab var seqbont "Sequential Bonferroni t";
lab var seqbonp "Sequential Bonferroni P";
format refdunt seqbont %8.2f;
format refdunp seqbonp %9.2e;
desc;
list;

gene p=0.1*_n/_N;

* Compare Holm procedure with seqbonp *;
multproc,method(holm) punc(0.1) cr(holmp);
format holmp %9.2e;
gene holmrd=reldif(holmp,seqbonp);
assert holmrd<=1e-3;
desc;
list;
graph7 holmp seqbonp seqbonp,s(o.) c(.L);more;

* Compare Holland-Copenhaver procedure with refdunp *;
multproc,method(holland) punc(0.1) cr(hollp);
format hollp %9.2e;
gene hollrd=reldif(hollp,refdunp);
assert holmrd<=1e-2;
desc;
list;
graph7 hollp refdunp refdunp,s(o.) c(.L);more;

* Show that -by varlist:- works *;
clear;
local nbygp=20;
set obs `nbygp';
set seed 9876543210;
* True mean SND *;
local trumu=1;

expgen npinbygp=(_n*(_n+1))/2,oldseq(bygp) copyseq(pvseq);
lab var npinbygp "N of P-values in by-group";
lab var bygp "By-group";
lab var pvseq "P-value sequence within by-group";
sort bygp pvseq;
gene double ranexp=`trumu'*0.5*invchi2(2,uniform());
gene double rangau=invnorm(uniform());
gene double ranmix=rangau+ranexp;
gene double pvalmix=1-norm(ranmix);
lab var ranexp "Random exponential(mu=`trumu') deviate";
lab var rangau "Random Gaussian(0,1) deviate";
lab var ranmix "Random deviate from mixture distribution";
lab var pvalmix "P-value for testing mu=0";
desc;
for X in var ranexp rangau ranmix pvalmix:graph7 X,xlabel freq bin(20) \ more;

*
 Show that the -saving- option of -smileplot- works as it should
 (ie not in presence of -by varlist:-
 but in absence of -by varlist:-)
*;
tempfile tf0;
conf new file `tf0';
smileplot,pv(pvalmix) est(ranmix) meth(simes) t1(" ") saving(`tf0',replace);
conf file `tf0';
graph use `tf0';
erase `tf0';
conf new file `tf0';
by bygp,rc0:smileplot,pv(pvalmix) est(ranmix) meth(simes) t1(" ") saving(`tf0',replace);
conf new file `tf0';

* Show that -by varlist:- works in the absence of -saving- *;
foreach M in  bonferroni sidak holm holland liu1 liu2 hochberg rom simes yekutieli krieger {;
  disp _n as text "Analysis using by-groups with method: " as result "`M'";
  by bygp:multproc,pv(pvalmix) me(`M') rank(rank_`M') gpunc(gpunc_`M') gpcor(gpcor_`M') crit(crit_`M') rej(rej_`M') nh(nh_`M');
  savedresults save lastby r();
  assert inlist(rej_`M',0,1);assert inlist(nh_`M',0,1);
  assert rej_`M'==1-nh_`M';
  disp _n as text "Analysis using if-groups with method: " as result "`M'";
  forv i1=1(1)`nbygp' {;
    multproc if bygp==`i1',pv(pvalmix) me(`M') rank(rank) gpunc(gpunc) gpcor(gpcor) crit(crit) rej(rej) nh(nh);
    if `i1'==`nbygp' {;
      disp _n as text "Comparing saved results of last if-group with saved results from last by-group:";
      savedresults compare lastby r(),verbose;
    };
    foreach X of var rank gpunc gpcor crit rej nh {;
      assert `X'==`X'_`M' if bygp==`i1';
      assert `X'==. if bygp!=`i1';
    };
    drop rank gpunc gpcor crit rej nh;
  };
  disp _n as text "Repeat the analysis using the -gpcuncor- and -gpcor- variables";
  by bygp,rc0:smileplot,pv(pvalmix) est(ranmix) punc(gpunc_`M') pcor(gpcor_`M') me(`M') rank(rank) gpunc(gpunc) gpcor(gpcor) crit(crit) rej(rej) nh(nh);
  foreach X of var rank gpunc gpcor rej nh {;
        assert `X'==`X'_`M';
  };
  assert crit==gpcor;
  drop rank gpunc gpcor crit rej nh;
  savedresults drop lastby;
};

exit;