/*  
*******************************************************************************
This code is useful for placing the output into an tables in an rtf file. Note that 
lines 260 and 261 (bottom) are needed to complete this output process.

The results of the macro runs on line 13 would be automatically stored as
temp_output.rtf in the folder "designs" on the z drive.  Change names to meet your needs.
If you don't, the program still runs and you can save the output by going to the file
command menu (top left).  You will also be also be prompted at the end of the
program if you want to open or save the file.  The save option allows you to save it 
under any name in any folder.  
*/
ods rtf file='z:\designs\temp_output.rtf' style=Journal ;

%macro gampower(n1, n2,alpha, reps, seed1,seed2,a1,b1,a2,b2);

/*
**********************************************************************************
This macro provides for a simulation of the actual power for two-sided p-values 
associated with the Wilcoxon Test, the T-test, and the Satterthwaite corrected Z-test
for the Gamma Distribution.  
************************************************************************************
SAS Programmers who are interested in generating data from other distributions
rad this.

If one is familiar with SAS programming, this Macro can be readily altered to generate
samples from other distributions (Normal, Beta, uniform, binomial,Poisson,and F for 
example).

There are four changes needed. 
1. The Macro call on line 15 may have more or less parameters than two per group.
2. Cosmetic:Means and standard deviations need to be calculated (See lines 95-100, 209,)
3. The data generation must be changed from Gamma to your preference. See
lines 106-108 and 118-120.
4. The Macro call lines 256 and 258 need to be tuned to line 15.
********************************************************************************
The user supplies the two sample sizes, the type I error rate (alpha), and the two 
cumulative distributions where a power calculation is needed.  If the two 
distributions are the same,the macro will estimate the true Type I error for each 
of the three tests, whose intended intended type I error is alpha (say 0.05). 


********************************************************************************
Macro Usage:


%macro gampower(n1, n2,alpha, reps, seed1,seed2,a1,b1,a2,b2);

gamma density for population j (j=1,2)is 

fj(x)=gamma(aj+1)[bj**(aj)[x**(aj-1)]exp(-x/bj)     x>0.

If aj=1, the distribution is exponential.

If 2aj is an integer and bj=2, then we have a chi-square distribution with 
2aj degrees of freedom.


*******************************************************************************
Input parameters: (Used in the example below)
n1=sample size for the first group=30
n2=sample size for the second group=30
alpha=.05 (Intended Type I error rate=5%)
reps=Number of simulations=100000
seed1 and seed2 large odd numbers to seed the randomization (<2**31) (6371 and 7431)
a1, b1 shape and scale parameters for first group (0.5 and 2 respectively)
a2, b2 shape and scale parameters for first group (0.5 and 2 respectively) (Null case)

%gampower(30,30,.05,100000,6371,7431,.5,2,.5,2);
*****************************************************************************
Output: 
It estimates the probability that the Wilcoxon, t-statistic, and Satterthwaite 
corrected Z-statistic in absolute value exceed the tabulated cutoffs for 
the nominal value alpha.  Standard errors for these estimates are also provided.

The population means and standard deviations are given as well.

Under the alternate hypothesis, it also provides information on the sample size
requirements to obtain 80% and 90% power at P=alpha (user supplied), two-sided, 
for the Wilcoxon and T-tests.


Note: If the means under the null and alternate are the same, the asymptotic power
calculations are not done. 

The output also summarizes the input information.




*/
data a;n1=&n1;n2=&n2;alpha=&alpha;reps=&reps;seed1=&seed1;
seed2=&seed2;a1=&a1;b1=&b1;a2=&a2;b2=&b2;

/*
Means and standard deviations for the Gamma are calculated here.
*/
mean=a1*b1;meana=a2*b2;
sd_null=b1*sqrt(a1);
sd_alt=b2*sqrt(a2);

data b;set a;
x=0;
do i1=1 to reps;
do i2=1 to n1;
/*
Samples for Group 1 are generated here
*/
call rangam(seed1,a1,y);
y=b1*y;
keep i1 x y alpha;
output;end;
end;
data b2;set a;
do i1=1 to reps;
do i2=1 to n2;
x=1;
/*
Samples for Group 2are generated here
*/
call rangam(seed2,a2,y);
y=b2*y;
keep i1 x y alpha;
output;end;
end;


data b;set b b2;
proc sort;by i1 x;
data b;set b;
proc rank ties=mean;;by i1;var y;ranks ry;
proc means noprint;by i1 x;var y ry alpha;
output out=new mean=y ry alpha std=sy sry sa n=ny nny na;
data new;set new;
data new0;set new;if x=0;
rybar1=ry;ybar1=y;
srybar1=sry;sybar1=sy;
n1=ny;
keep i1 n1 rybar1 ybar1 srybar1 sybar1 alpha;
data new1;set new;if x=1;
rybar2=ry;ybar2=y;
srybar2=sry;sybar2=sy;
n2=ny;
keep i1 n2 rybar2 ybar2 srybar2 sybar2;
data new;merge new0 new1;by i1;
data new;set new;
n=n1+n2;
wt1=((sybar1**2)/n1);
wt2=((sybar2**2)/n2);
df2=((wt1+wt2)**2)/((wt1*wt1/(n1-1))+(wt2*wt2/(n2-1)));
sse=((sybar1**2)*(n1-1))+((sybar2**2)*(n2-1));
se=sse/(n-2);
se=sqrt(se*((1/n1)+(1/n2)));
if se<.00000001 then se=.00000001;
t=abs(ybar1-ybar2)/se;
sse_r=((srybar1**2)*(n1-1))+((srybar2**2)*(n2-1));
meanr=(n1*rybar1+n2*rybar2)/n;
ss_tr=n1*((rybar1-meanr)**2)+n2*((rybar2-meanr)**2);
ss_tot_r=ss_tr+sse_r;
var=ss_tot_r/(n-1);
ser=sqrt(var*n2/(n1*n));
if ser<.00000001 then ser=.00000001;
z=abs(rybar1-meanr)/ser;
seu=sqrt(wt1+wt2);
if seu<.00000001 then seu=.00000001;
zz=abs(ybar1-ybar2)/seu;
countw=0;countt=0;counts=0;counttz=0;countzu=0;
if z>-probit(alpha/2) then countw=1;
if t>-tinv(alpha/2,n-2) then countt=1;
if zz>-tinv(alpha/2,df2) then counts=1;
if t>-probit(alpha/2) then counttz=1;
if zz>-probit(alpha/2) then countzu=1;
proc means maxdec=5 mean n stderr noprint;
var countw countt counts counttz countzu; 
output out=d mean=p_wil p_t p_z p_t_normal p_zu_normal
stderr=se_wil se_t se_z se_t_normal se_zu_normal;
data d;set d;xxx=1;proc sort;by xxx;;
data a;set a;xxx=1;proc sort;by xxx;
data a;merge a d;by xxx;
if p_t<.99 then
t_factor80=((.84-probit(alpha/2))/(probit(p_t)-probit(alpha/2)))**2;
if p_wil<.99 then
w_factor80=((.84-probit(alpha/2))/(probit(p_wil)-probit(alpha/2)))**2;
if p_t<.99 then
t_factor90=((1.282-probit(alpha/2))/(probit(p_t)-probit(alpha/2)))**2;
if p_wil<.99 then
w_factor90=((1.282-probit(alpha/2))/(probit(p_wil)-probit(alpha/2)))**2;
n1t_80=int(.5+n1*t_factor80);
n1t_90=int(.5+n1*t_factor90);
n1w_80=int(.5+n1*w_factor80);
n1w_90=int(.5+n1*w_factor90);
n2t_80=int(.5+n2*t_factor80);
n2t_90=int(.5+n2*t_factor90);
n2w_80=int(.5+n2*w_factor80);
n2w_90=int(.5+n2*w_factor90);
if abs(mean-meana)>.00001*sd_null then do;
nt80=((.84-probit(alpha/2))/(mean-meana))**2;
nt80=int(nt80*(sd_null**2+sd_alt**2)+.5);
nt90=((1.282-probit(alpha/2))/(mean-meana))**2;
nt90=int(nt90*(sd_null**2+sd_alt**2)+.5);
end;
 

file print; 

put // @10 'Power Simulation of Wilcoxon, T-test, Satterthwaite Corrected T';

put /
@10 'Gamma Distribution with Parameters a1=' a1 '  b1=' b1 '   a2=' a2 '  b2=' b2;
put / @10 'Mean under Null=' mean best6.  @36 'Mean under Alternative=' meana best6.;
put / @10 'SD under Null=' sd_null best6. @36 'SD under Alternative=' sd_alt best6.;
q=n1+n2;
put // @10 'Total sample size=' q '  Group 1 N=' n1 '  Group 2 N=' n2
// @10 'Based on ' reps ' simulations with seeds=' seed1 'and ' seed2
/// @10 'Nominal Alpha=' alpha '   Two-sided'
// @10 'Statistical Test' @30 'Estimated Power' @50 'Standard Error'
/ @10 'Rank Test'
/ @10 'Wilcoxon' @30 p_wil @50 se_wil best7.
// @10 'Parametric Tests'
/@10 'Pooled Variance'
/ @10 'T-test' @30 p_t @50 se_t best7.   @60 '(T Cutpoint)'
/ @10 'T-test' @30 p_t_normal @50 se_t_normal best7. @60 '(Normal Cutpoint)'
//@10 'Unpooled Variance'
/ @10 'Z,Unpooled Var' @30 p_z @50 se_z best7. @60 '(Satterthwaite)'
/ @10 'Z,Unpooled Var' @30 p_zu_normal @50 se_zu_normal best7.
@60 '(Normal Cutpoint)';

 if abs(mean-meana)>.01*sd_alt then do;
 put _page_;
 put // @10 'Sample Size Projections for 80% and 90% power';
 put //@10 'Method=Browniam Motion';
put //@10 'Actual Sample Sizes=' n1 'and ' n2;
put //@10 'Type I error=' alpha ' Two-sided';
put / @10 'Test' @25 '80% Power'  @55 '90% Power';
put / @10 'Wilcoxon' @25 'n1=' n1w_80  ' and  n2=' n2w_80  
@55 'n1=' n1w_90 ' and n2=' n2w_90 ;
put / @10 'T-test' @25 'n1=' n1t_80 ' and  n2=' n2t_80 @55 'n1=' n1t_90 ' and  n2=' n2t_90 ;
if p_t<.20 then put 
// @10 'Warning: Power for t is too low for accurate t-test Brownian Motion projection'
// @10 'Try larger sample sizes, say double, to get new projection'
// @10 'T-requirement by confidence interval inversion=' nt80 'per group for 80% power'
// @10 'T-requirement by confidence interval inversion=' nt90 'per group for 90% power';

if p_wil<.20 then put 
// @10 'Warning: Power is too low for accurate Wilcoxon Brownian Motion projection'
// @10 'Try larger sample sizes, say double, to get new projection';

end;


run;
%mend;


title 'Gamma Macro Null';
%gampower(30,30,.05,100000,6371,7431,.5,2,.5,2);
title 'Gamma Macro Alternative';
%gampower(30,30,.05,100000,6371,7431,.5,2,.5,6.4);

ods rtf close;
run;