/********************************************************************
*								    *
*   	Macro and sample program of how to run the various step-down*
* 	procedures from Tamhane and Logan (2000) to identify the    *
* 	minimum effective dose and maximum safe dose (or therapeutic*
*	window).  Specify the following information:		    *
* 		- medmsd: directory where data file is stored	    *
*		- storeres: directory where results should be stored*
*		- datafile: should have 3 columns with dose (as     *
*		  integer from 0 (control) to k), efficacy response,*
*		  and safety response.				    *
*								    *
*	The macros available are:				    *
*		- %runboot, to run the exact bootstrap SD procedures*
*		- %runbonf, to run the approximate Bonferroni SD    *
*		  procedures					    *
*		- %unionsd, to run the SD procedure based on the    *
*		  union hypothesis testing formulation		    *
*								    *
*	The commands to run these macros are 			    *		 
*  %macro runboot(k,nboot,alpha,dosd1,dosd2,dat,delta1,delta2,printres);*
*  %macro runbonf(k,      alpha,dosd1,dosd2,dat,delta1,delta2);     *
*  %macro unionsd(k,      alpha,dat,delta1,delta2);		    *
*								    *
*   	The parameters are:					    *
* 		- k: number of doses in addition to control	    *
*		- nboot: number of bootstrap samples		    *
*		- alpha: type I error rate			    *
*		- dosd1: 1=run SD1 procedure, 0=don't		    *
*		- dosd2: 1=run SD2 procedure, 0=don't		    *
*		- dat: directory and name of datafile, medmsd.filename *
*		- delta1: efficacy mean threshold, so that a dose   *
*		  is effective if it exceeds the mean of the control*
*		  dose by at least delta1.			    *
*		- delta2: safety mean threshold, so that a dose is  *
*		  safe if it does not exceed the mean of the control*
*		  dose by at least delta2			    *
*		- printres: 1=print details of exact SD procedures, *
*		  including means, SD's of each dose, pooled SD's,  *
*		  degrees of freedom, t statistics and p-values for *
*		  each comparison to control, and the adjusted      *
*		  p-values for efficacy and safety at each step of  *
*		  the step-down testing procedures.  0=don't print  *
*								    *
*	The output is a dataset called storeres.winres with the     *
*	methods used, and the minimum effective dose and maximum    *
*	safe dose for each method.  This can be printed out using   *
* 	the command demonstrated below.				    *
*								    *
*	Following is a sample program from Tamhane and Logan (2000):*
*	This should be put in a new file, and will reference the    *
*	medmsdmacro.sas file with the %include command.	  	    *
*								    *
*								    *
* libname medmsd '/.../nwu.edu/fs/home/b/brl866/d99/medmsddata';    *
* libname storeres '/.../nwu.edu/fs/home/b/brl866/d99/medmsddata';  *
* data storeres.winres;						    *
*   input METHOD $ 1-15 MED MSD ;				    *
*   cards;;							    *
* run;								    *
*								    *
* %include '/.../nwu.edu/fs/home/b/brl866/d99/medmsdmacro.sas';	    *
*								    *
* %runboot(4,100,0.05,1,1,medmsd.womz,0.5,3.0,1);		    *
* %runbonf(4,     0.05,1,1,medmsd.womz,0.5,3.0);  		    *
* %unionsd(4,  0.05,medmsd.womz,0.5,3.0);			    *
*								    *
* proc print data=storeres.winres;				    *
* title1 '-----------------------------------------------------------------';*
* title2 'Minimum Effective Dose and Maximum Safe Dose (Therapeutic Window)';*
* run;								    *
*								    *
*								    *
********************************************************************/



/*********************************************************
*							 *
*	Macro to run bootstrap exact versions of SD1 and *
*     	SD2.  Assumes common covariance matrix across    *
* 	dose.						 *
*							 *
*********************************************************/

%macro runboot(k,nboot,alpha,dosd1,dosd2,dat,delta1,delta2,printres);

proc sort data=&dat;
by dose;
run;

/* Create dataset storing sample sizes and residuals*/

proc means data=&dat  n mean stderr;
by dose;
var eff safe;
output out=dosemn n=ne ns mean=mne mns var=vare vars;
run;
proc sort data=dosemn;
by dose;
run;

data resid;merge &dat dosemn(keep=dose mne mns);
by dose;
rese=eff-mne;
ress=safe-mns;
run;
data resid;set resid(keep=dose rese ress);
run;

/*create dataset storing means to compute raw p-values*/

data mntest; set dosemn(keep=dose mne mns vare vars);
run;

/*open IML*/
proc iml;

/* read datasets into IML matrices*/
use mntest;
read all var {dose mne mns vare vars} into means;
use resid;
read all var {dose rese ress} into resid;
use dosemn;
read all var {ne} into n;

/* compute raw p-values*/
rawpe=j(1, &k,.);
rawps=j(1,&k,.);
nmin1=n-j(&k+1,1,1);
df=nmin1[+,1];
sp2e=nmin1`*means[,4]/df;
sp2s=nmin1`*means[,5]/df        ;

if (&printres>0) then do;
  print "--------------------------------------------------------------------";
  print "Pooled Variance for Efficacy" sp2e;
  print "Pooled Variance for Safety" sp2s;
  print "Degrees of Freedom" df;
end;

do i=1 to &k;
  te=(means[i+1,2]-means[1,2]-&delta1)/
       sqrt(sp2e*(1/n[i+1,1]+1/n[1,1]));
  ts=(means[1,3]-means[i+1,3]+&delta2)/
       sqrt(sp2s*(1/n[i+1,1]+1/n[1,1]));
  rawpe[1,i]=1-probt(   te ,df);
  rawps[1,i ]=1-probt(    ts ,df);
  pvaleff=rawpe[1,i];
  pvalsafe=rawps[1,i];
  if (&printres>0) then do;
    print "------------------------------------------------------------------";
    print "Dose" i "comparison to Control:";
    print "  Efficacy t statistic: " te "Efficacy p-value: " pvaleff;
    print "  Safety t statistic: " ts "Safety p-value: " pvalsafe;
  end;
end;

/*Now we generate the bootstrap samples*/
/*first do some steps that we don't want to do in each loop*/
eff=j(&k+1,n[<>,1],.);
safe=j(&k+1,n[<>,1],.);
ntot=n[+,1];
pvale=j(&nboot,2*&k,.);
pvals=j(&nboot,2*&k,.);


do m=1 to &nboot;
  do i=1 to (&k+1);
    do j=1 to n[i,1];
      randnum=ceil(ranuni(0)*ntot);
      eff[i,j]=resid[randnum,2];
      safe[i,j]=resid[randnum,3];
    end;
  end;


  effmean=eff[,:];
  safemean=safe[,:];
  effvar=(eff[,##]-effmean##2#n )/nmin1;
  safevar=(safe[,##]-safemean##2#n )/nmin1;
  sp2e=nmin1`*effvar/df;
  sp2s=nmin1`*safevar/df;

/*compute joint p-values at mth bootstrap*/

  do i=1 to (&k  );
    tie=(effmean[i+1,1]-effmean[1,1]        )/
           sqrt(sp2e*(1/n[i+1,1]+1/n[1,1]));
    tis=(safemean[1,1]-safemean[i+1,1]        )/
           sqrt(sp2s*(1/n[i+1,1]+1/n[1,1]));
    pvale[m,i]=1-probt(tie,df);
    pvals[m,i]=1-probt(tis,df);
  end;
/*compute joint min pvalues at mth bootstrap for Min SD procedure*/
  pvale[m,&k+1]=pvale[m,1];
  pvals[m,2*&k]=pvals[m,&k];
  do i=2 to &k;
    pvale[m,&k+i]=min(pvale[m,&k+i-1],pvale[m,i  ]);
    pvals[m,2*&k-i+1]=min(pvals[m,2*&k-i+2],pvals[m,&k-i+1]);
  end;
end;

window=j(1,2,.);

/* Compute adjusted p-values using SD1 procedure*/

if (&dosd1=1) then do;
reject=1;
teste=1;
tests=1;
/* Base case*/
/* Compute adjusted p-values*/
cnte=0; cnts=0;
/*compute minimums and index them*/
rawpemin=j(2,&k,rawpe[1,1]);
rawpsmin=j(2,&k,rawps[1,&k]);
rawpemin[1,1]=1;
rawpsmin[1,&k]=&k;
do i=2 to &k;
  if (rawpe[1,i]<rawpemin[2,i-1]) then do;
    rawpemin[1,i]=i;
    rawpemin[2,i]=rawpe[1,i];
    end;
  else do;
    rawpemin[1,i]=rawpemin[1,i-1];
    rawpemin[2,i]=rawpemin[2,i-1];
    end;
  if (rawps[1,&k-i+1]<rawpsmin[2,&k-i+2]) then do;
    rawpsmin[1,&k-i+1]=&k-i+1;
    rawpsmin[2,&k-i+1]=rawps[1,&k-i+1];
    end;
  else do;
    rawpsmin[1,&k-i+1]=rawpsmin[1,&k-i+2];
    rawpsmin[2,&k-i+1]=rawpsmin[2,&k-i+2];
    end;
end;


do m=1 to &nboot;
  pmin=min(pvale[m,2*&k],pvals[m,&k+1]);
  if (pmin<rawpemin[2,&k]) then cnte=cnte+1;
  if (pmin<rawpsmin[2,1]) then cnts=cnts+1;
end;
sdadjpe=cnte/&nboot;
sdadjps=cnts/&nboot;

/* Test initial hypothesis*/
rn=&k; sn=1; ro=&k; so=1;
if (&printres>0) then do;
  print "-----------------------------------------------------------------------";
  print "P-values for each stage of the step-down procedure SD1";
  print "Efficacy Dose: " rn "Safety Dose: " sn;
  print  "Adj. p-value for Efficacy: " sdadjpe "Adj. p-value for safety:" sdadjps;
end;
if (sdadjpe<=&alpha) then do;
    ro=rn;
    rn=rawpemin[1,&k]-1;

    end;
else teste=0;
if (sdadjps<=&alpha) then do;
    so=sn;
    sn=rawpsmin[1,1]+1;

    end;
else tests=0;
if ((rn=0)&(sn=(&k+1))) then do;
  teste=0;tests=0;
  end;

/* Loop*/
do while((teste|tests)=1);
  cnte=0; cnts=0;teste=1;tests=1;
if ((rn=0   )&(sn<=&k) ) then do;
    teste=0;
    do m=1 to &nboot;
      pmin=    pvals[m,&k+sn]          ;
      if (pmin<rawpsmin[2,sn]) then cnts=cnts+1;
    end;
    sdadjps=cnts/&nboot;
    if (&printres>0) then do;
      print "Efficacy Dose: " rn "Safety Dose: " sn;
      print  "Adj. p-value for safety:" sdadjps;
    end;
    if (sdadjps<=(&alpha  )) then do;
      so=sn;
      sn=rawpsmin[1,so]+1;
      if sn=&k+1 then tests=0;
      end;
    else tests =0;
  end;
  else if ((sn=(&k+1))&(rn>0 )) then do;
    tests=0;
    do m=1 to &nboot;
      pmin=    pvale[m,&k+rn]          ;
      if (pmin<rawpemin[2,rn]) then cnte=cnte+1;
    end;
    sdadjpe=cnte/&nboot;
    if (&printres>0) then do;
      print "Efficacy Dose: " rn "Safety Dose: " sn;
      print  "Adj. p-value for Efficacy: " sdadjpe ;
    end;
    if (sdadjpe<=(&alpha  )) then do;
      ro=rn;
      rn=rawpemin[1,ro]-1;
      if rn=0 then teste=0;
      end;
    else teste =0;
  end;
  else do;
    do m=1 to &nboot;
      pmin=min(pvale[m,&k+rn],pvals[m,&k+sn]);
      if (pmin<rawpemin[2,rn]) then cnte=cnte+1;
      if (pmin<rawpsmin[2,sn]) then cnts=cnts+1;
    end;
    sdadjpe=cnte/&nboot;
    sdadjps=cnts/&nboot;
    if (&printres>0) then do;
      print "Efficacy Dose: " rn "Safety Dose: " sn;
      print  "Adj. p-value for Efficacy: " sdadjpe "Adj. p-value for Safety: " sdadjps ;
    end;
    if (sdadjpe<=&alpha) then do;
        ro=rn;
        rn=rawpemin[1,ro]-1;

        end;
    else teste=0;
    if (sdadjps<=&alpha) then do;
        so=sn;
        sn=rawpsmin[1,so]+1;

        end;
    else tests=0;
    if ((rn=0)&(sn=(&k+1))) then do;
      teste=0;tests=0;
      end;
  end;
end;
temppr1=rn+1;
temppr2=sn-1;
/*print temppr1 temppr2;*/
method={Bootstrap_SD1};
window[1,1]=rn+1;
window[1,2]=sn-1;

end; /*Do Loop for SD1 procedure*/

edit storeres.winres;
  append from window [rowname=method];

/* Compute adjusted p-values using SD2 procedure*/
if (&dosd2=1) then do;
reject=1;
teste=1;
tests=1;
/* Base case*/
/* Compute adjusted p-values*/
cnte=0; cnts=0;
do m=1 to &nboot;
  pmin=min(pvale[m,&k],pvals[m,1]);
  if (pmin<rawpe[1,&k]) then cnte=cnte+1;
  if (pmin<rawps[1,1]) then cnts=cnts+1;
end;
sdadjpe=cnte/&nboot;
sdadjps=cnts/&nboot;


/* Test initial hypothesis*/
rn=&k; sn=1; ro=&k; so=1;
if (&printres>0) then do;
  print "------------------------------------------------------------------------";
  print "P-values for each stage of the step-down procedure SD2";
  print "Efficacy Dose: " rn "Safety Dose: " sn;
  print  "Adj. p-value for Efficacy: " sdadjpe "Adj. p-value for safety:" sdadjps;
end;
if (sdadjpe<=&alpha) then do;
    ro=rn;
    rn=ro-1;
    end;
else teste=0;
if (sdadjps<=&alpha) then do;
    so=sn;
    sn=so+1;
    end;
else tests=0;


/* Loop*/

do while((teste|tests)=1);
  cnte=0; cnts=0;teste=1;tests=1;
  if ((rn=0)&(sn<=&k)    ) then do;
    teste=0;
    do m=1 to &nboot;
      pmin=    pvals[m,sn]             ;
      if (pmin<rawps[1,sn]) then cnts=cnts+1;
    end;
    sdadjps=cnts/&nboot;
    if (&printres>0) then do;
      print "Efficacy Dose: " rn "Safety Dose: " sn;
      print  "Adj. p-value for safety:" sdadjps;
    end;
    if (sdadjps<=(&alpha  )) then do;
      so=sn;
      sn=so+1;
      if sn=&k+1 then tests=0;
      end;
    else tests =0;
  end;
  else if ((sn=(&k+1))&(rn>0) ) then do;
    tests=0;
    do m=1 to &nboot;
      pmin=    pvale[m,rn]             ;
      if (pmin<rawpe[1,rn]) then cnte=cnte+1;
    end;
    sdadjpe=cnte/&nboot;
    if (&printres>0) then do;
      print "Efficacy Dose: " rn "Safety Dose: " sn;
      print  "Adj. p-value for Efficacy: " sdadjpe ;
    end;
   if (sdadjpe<=(&alpha)) then do;
      ro=rn;
      rn=ro-1;
      if rn=0 then teste=0;
      end;
    else teste =0;
  end;
  else do;
    do m=1 to &nboot;
      pmin=min(pvale[m,rn],pvals[m,sn]);
      if (pmin<rawpe[1,rn]) then cnte=cnte+1;
      if (pmin<rawps[1,sn]) then cnts=cnts+1;
    end;
    sdadjpe=cnte/&nboot;
    sdadjps=cnts/&nboot;
    if (&printres>0) then do;
      print "Efficacy Dose: " rn "Safety Dose: " sn;
      print  "Adj. p-value for Efficacy: " sdadjpe "Adj. p-value for Safety: " sdadjps ;
    end;
    if (sdadjpe<=&alpha) then do;
        ro=rn;
        rn=ro-1;

        end;
    else teste=0;
    if (sdadjps<=&alpha) then do;
        so=sn;
        sn=sn+1;

        end;
    else tests=0;
    if ((rn=0)&(sn=(&k+1))) then do;
      teste=0;tests=0;
      end;
  end;


end;
temppr1=rn+1;
temppr2=sn-1;
/* print temppr1 temppr2;*/
method={Bootstrap_SD2};
window[1,1]=rn+1;
window[1,2]=sn-1;
end; /*Do Loop for SD2 procedure*/
edit storeres.winres;
  append from window [rowname=method];


quit; /*IML*/

%mend runboot;



/************************************************************
*   							    *
*   This macro applies the Bonferroni version of SD1 and    *
*   SD2 procedures.  Uses large sample Dunnett critical     *
*   values, so that modifications are required for smaller  *
*   sample sizes.					    *
*							    *
************************************************************/

%macro runbonf(k,      alpha,dosd1,dosd2,dat,delta1,delta2);

proc sort data=&dat;
by dose;
run;
/* Create dataset storing sample sizes and residuals*/
proc means data=&dat n mean stderr noprint;
by dose;
var eff safe;
output out=dosemn n=ne ns mean=mne mns var=vare vars;
run;
proc sort data=dosemn;
by dose;
run;

/*create dataset storing means to compute raw p-values*/
data mntest; set dosemn(keep=dose mne mns vare vars);
run;

/*open IML*/
proc iml;
/* read datasets into IML matrices*/
use mntest;
read all var {dose mne mns vare vars} into means;
use dosemn;
read all var {ne} into n;
/* compute raw p-values*/
rawpe=j(1, &k,.);
rawps=j(1,&k,.);
nmin1=n-j(&k+1,1,1);
df=nmin1[+,1];
sp2e=nmin1`*means[,4]/df;
sp2s=nmin1`*means[,5]/df        ;
/* print df sp2e sp2s;*/

do i=1 to &k;
  te=(means[i+1,2]-means[1,2]-&delta1)/
       sqrt(sp2e*(1/n[i+1,1]+1/n[1,1]));
  ts=(means[1,3]-means[i+1,3]+&delta2)/
       sqrt(sp2s*(1/n[i+1,1]+1/n[1,1]));
  rawpe[1,i]=1-probt(   te ,df);
  rawps[1,i ]=1-probt(    ts ,df);
  temppr1=rawpe[1,i];
  temppr2=rawps[1,i];
/*  print te ts temppr1 temppr2;*/
end;

/* if sample is large then critical values are */

crit025={1.96 2.212 2.349 2.442 2.512};
/*  print crit025;*/
/* if n=10, then */
/*
crit025={2.0049 2.2715 2.4170 2.5161 2.5907};
*/

/* Otherwise, look up the appropriate Dunnett critical values
and insert them here */

crit05={1.645 1.916 2.062 2.160 2.234};
pcrit025=1-probt(crit025,df);
pcrit05=1-probt(crit05,df);

window=j(1,2,.);


/* Compute adjusted p-values using SD1 procedure*/
if (&dosd1=1) then do;
reject=1;
teste=1;
tests=1;
/* Base case*/
/* Compute adjusted p-values*/
cnte=0; cnts=0;
/*compute minimums and index them*/
rawpemin=j(2,&k,rawpe[1,1]);
rawpsmin=j(2,&k,rawps[1,&k]);
rawpemin[1,1]=1;
rawpsmin[1,&k]=&k;
do i=2 to &k;
  if (rawpe[1,i]<rawpemin[2,i-1]) then do;
    rawpemin[1,i]=i;
    rawpemin[2,i]=rawpe[1,i];
    end;
  else do;
    rawpemin[1,i]=rawpemin[1,i-1];
    rawpemin[2,i]=rawpemin[2,i-1];
    end;
  if (rawps[1,&k-i+1]<rawpsmin[2,&k-i+2]) then do;
    rawpsmin[1,&k-i+1]=&k-i+1;
    rawpsmin[2,&k-i+1]=rawps[1,&k-i+1];
    end;
  else do;
    rawpsmin[1,&k-i+1]=rawpsmin[1,&k-i+2];
    rawpsmin[2,&k-i+1]=rawpsmin[2,&k-i+2];
    end;
end;


/* Test initial hypothesis*/
rn=&k; sn=1; ro=&k; so=1;
if (rawpemin[2,rn]<=pcrit025[1,rn]) then do;
    ro=rn;
    rn=rawpemin[1,&k]-1;
    if (rn=0) then teste=0;
    end;
else teste=0;
if (rawpsmin[2,sn]<=pcrit025[1,&k-sn+1]) then do;
    so=sn;
    sn=rawpsmin[1,1]+1;
    if (sn=&k+1) then tests=0;
    end;
else tests=0;

/* Loop*/
do while((teste|tests)=1);
if ((teste=0)&(tests=1)) then do;
    if (rawpsmin[2,sn]<=pcrit025[1,&k-sn+1]) then do;
      so=sn;
      sn=rawpsmin[1,so]+1;
      if sn=&k+1 then tests=0;
      end;
    else tests =0;
  end;
  else if ((teste=1)&(tests=0)) then do;
    if (rawpemin[2,rn]<=pcrit025[1,rn]) then do;
      ro=rn;
      rn=rawpemin[1,ro]-1;
      if rn=0 then teste=0;
      end;
    else teste =0;
  end;
  else do;
    if (rawpemin[2,rn]<=pcrit025[1,rn]) then do;
        ro=rn;
        rn=rawpemin[1,ro]-1;
        if rn=0 then teste=0;
        end;
    else teste=0;
    if (rawpsmin[2,sn]<=pcrit025[1,&k-sn+1]) then do;
        so=sn;
        sn=rawpsmin[1,so]+1;
        if sn=&k+1 then tests=0;
        end;
    else tests=0;
  end;
end;

method={Bonferroni_SD1};
window[1,1]=rn+1;
window[1,2]=sn-1;

end; /*Do Loop for SD1 procedure*/

edit storeres.winres;
  append from window[rowname=method];

/* Compute adjusted p-values using SD2 procedure*/

if (&dosd2=1) then do;
reject=1;
teste=1;
tests=1;
/* Base case*/
/* Test initial hypothesis*/
rn=&k; sn=1; ro=&k; so=1;
if (rawpe[1,rn]<=pcrit025[1,1 ]) then do;
    ro=rn;
    rn=ro-1;
    end;
else teste=0;
if (rawps[1,sn]<=pcrit025[1,1      ]) then do;
    so=sn;
    sn=so+1;
    end;
else tests=0;
/* Loop*/
do while((teste|tests)=1);
  cnte=0; cnts=0;
  if ((teste=0)&(tests=1)) then do;
    if (rawps[1,sn]<=pcrit025[1,1      ]) then do;
      so=sn;
      sn=so+1;
      if sn=&k+1 then tests=0;
      end;
    else tests =0;
  end;
  else if ((teste=1)&(tests=0)) then do;
    if (rawpe[1,rn]<=pcrit025[1,1 ]) then do;
      ro=rn;
      rn=ro-1;
      if rn=0 then teste=0;
      end;
    else teste =0;
  end;
  else do;
    if (rawpe[1,rn]<=pcrit025[1,1 ]) then do;
        ro=rn;
        rn=ro-1;
        if rn=0 then teste=0;
        end;
    else teste=0;
    if (rawps[1,sn]<=pcrit025[1,1      ]) then do;
        so=sn;
        sn=sn+1;
        if sn=&k+1 then tests=0;
        end;
    else tests=0;
  end;
end;
method={Bonferroni_SD2};
window[1,1]=rn+1;
window[1,2]=sn-1;
end; /*Do Loop for SD2 procedure*/
edit storeres.winres;
  append from window[rowname=method];

quit; /*IML*/

%mend runbonf;


/***********************************************************
*							   *
*	Macro to run SD3 procedure, based on the union     *
* 	hypothesis formulation, which does not assume weak *
*	monotonicity.					   *
*							   *
***********************************************************/

%macro unionsd(k,      alpha,dat,delta1,delta2);

proc sort data=&dat;
by dose;
run;
/* Create dataset storing sample sizes and residuals*/
proc means data=&dat n mean stderr noprint;
by dose;
var eff safe;
output out=dosemn n=ne ns mean=mne mns var=vare vars;
run;
proc sort data=dosemn;
by dose;
run;

/*create dataset storing means to compute raw p-values*/
data mntest; set dosemn(keep=dose mne mns vare vars);
run;

/*open IML*/
proc iml;
/* read datasets into IML matrices*/
use mntest;
read all var {dose mne mns vare vars} into means;
use dosemn;
read all var {ne} into n;
/* compute raw p-values*/
rawpe=j(1, &k,.);
rawps=j(1,&k,.);
maxp=j(1,&k,.);
indmat=j(1,&k,.);
indrev=j(1,&k,.);
minmaxp=j(1,&k,.);
nmin1=n-j(&k+1,1,1);
df=nmin1[+,1];
sp2e=nmin1`*means[,4]/df;
sp2s=nmin1`*means[,5]/df        ;
/*print df sp2e sp2s;*/

med2=6;msd2=0;
window=j(1,2,.);
do i=1 to &k;
  te=(means[i+1,2]-means[1,2]-&delta1)/
       sqrt(sp2e*(1/n[i+1,1]+1/n[1,1]));
  ts=(means[1,3]-means[i+1,3]+&delta2)/
       sqrt(sp2s*(1/n[i+1,1]+1/n[1,1]));
  rawpe[1,i]=1-probt(   te ,df);
  rawps[1,i ]=1-probt(    ts ,df);
  temppr1=rawpe[1,i];
  temppr2=rawps[1,i];

/*  print te ts temppr1 temppr2;*/
end;

maxp=rawpe<>rawps;
indmat=rank(maxp);
/*print maxp indmat;*/
do i=1 to &k;
  minmaxp[1,indmat[1,i]]=maxp[1,i];
  indrev[1,indmat[1,i]]=i;
end;
/*print minmaxp indrev;*/
keeptest=1;
if (minmaxp[1,1]<=&alpha/&k) then do;
  med2=min(med2,indrev[1,1]);
  msd2=max(msd2,indrev[1,1]);
  j=2;
  end;
else keeptest=0;
do while (keeptest>0)   ;
  if (indrev[j]>med2 & indrev[j]<msd2) then j=j+1;
  else do;
    newalph=&alpha/(&k-1-(msd2-med2));
    if (minmaxp[1,j]<=newalph)                 then do;
      med2=min(med2,indrev[1,j]);
      msd2=max(msd2,indrev[1,j]);
      j=j+1;
      end;
    else keeptest=0;
    end;

  if (j=&k+1) then keeptest=0;
  end;

method={SD3};
window[1,1]=med2 ;
window[1,2]=msd2 ;
edit storeres.winres;
  append from window [rowname=method];
*/
quit; /*IML*/

%mend unionsd;