/* An example for calling logtwo SAS macro
%include 'logtwo.sas';
%logtwo(n1=10, n2=10,xb1=1.5,xb2=4.5, s1_sq=2.0, s2_sq=3.0, conlev=0.95);

run;


/****************************************************************************
 Name:  logtwo.sas

 This program compute the p-value for testing the difference between means
 of two lognormal distributions. The hypotheses are

         H_0: mean1 <= mean2  vs  H_a: mean1 > mean2

 Note: Consider two lognormal populations with parameters (u1, sigma1^2) and
 (u2,sigma2^2). Then mean1 = exp(u1+sigma1^2/2) and mean2 = exp(u2+sigma2^2/2)
 Let xbi = mean of the logged data from the ith population, i=1,2.
     si_sq = variance of the logged data from the ith population, i=1,2.
 For a given sample sizes n1, n2, xb1, xb2, s1_sq and s2_sq, this program
 computes the p-value for testing the hypotheses given above and confidence
 interval for the difference mean1 - mean2.


*****************************************************************************/
%macro logtwo(n1=, n2=, xb1=, s1_sq=,
        xb2=, s2_sq=, conlev=, m=100000);

options ls = 64 ps = 45 nodate nonumber;
title1 'Output of p-value and confidence interval ';
title2 ' for the difference of two lognormal means';

proc iml;

n1=&n1;
n2=&n2;
xb1=&xb1;
xb2=&xb2;
s1_sq=&s1_sq;
s2_sq=&s2_sq;
cl=&conlev;
m = &m;

alpha = 1.0-cl;
L1 = int(alpha*m/2.0);
L2 = int((1.0-alpha/2.0)*m);

df1 = n1 - 1.0;
sqrtdf1 = sqrt(df1);
df2 = n2 - 1.0;
sqrtdf2 = sqrt(df2);

s1 = sqrt(s1_sq);
s2 = sqrt(s2_sq);

const11 = s1/sqrt(n1);
const12 = s1_sq*df1*0.5;
const21 = s2/sqrt(n2);
const22 = s2_sq*df2*0.5;

seed1 = int(time());
seed2 = int(time()+12345);
seed3 = int(time()+13456);
seed4 = int(time()+19867);
pvalue = 0.0;

reset name=noname center=nocenter;
gv = j(m,1,0);
do j = 1 to m;
   z1 = rannor(seed1);
   v1 = 2.0*rangam(seed2,df1*0.5);
   z2 = rannor(seed3);
   v2 = 2.0*rangam(seed4,df2*0.5);

   T31 = xb1 + z1*const11*sqrtdf1/sqrt(v1)+const12/v1;
   T32 = xb2 + z2*const21*sqrtdf2/sqrt(v2)+const22/v2;
   gv[j] =exp(T31)-exp(T32);
   if gv[j] < 0.0 then pvalue = pvalue + 1.0;
end;
pvalue = pvalue/m;

gv0=gv;
gv[rank(gv)]=gv0;
gvl1=gv[L1];
gvl2=gv[L2];

print 'n1=' n1[format=2.0] ',' 'xb1=' xb1[format=6.3] ','
   's1^2=' s1_sq[format=6.3] ;
print 'n2=' n2[format=2.0] ',' 'xb2=' xb2[format=6.3] ','
   's2^2=' s2_sq[format=6.3] ;
print 'cl=' cl[format=4.2] ;

print 'p-value for testing mean_1 <= mean2 vs mean1 > mean2 is '
        pvalue [format=5.3];

print cl [format=4.2] 'Confidence Interval for mean1 - mean2 is',
        '(' gvl1 [format=12.4]', 'gvl2 [format=12.4]')';


quit;


%mend logtwo;