options ls=80 nodate pageno=1 NOFMTERR; 


 data key;
   set ph6470.midterm_2011;
   log_aer =  log(aer);
   male = (sex='M');
   if (aer NE . and whr NE . and sex ne '');
* also delete 2 outliers;
   if whr < 1.1;

 * Proc SGplot data=key;
 *  scatter x=whr y=aer / group=sex;

 Proc GLM data=key;
   title3 'Unadjusted comparison';
   class sex;
   model log_aer   = sex;
   lsmeans sex / pdiff CL;
   ods output LSMeanCL=logCI_crude;
* Proc SGplot data=key;
 *  reg x=whr y=log_aer / group=sex;

  Proc GLM data=key;
    title3 'reg adjusted, all data';
   class sex;
 *  model log_aer = whr sex whr*sex;
   model log_aer = whr sex;
   lsmeans sex / pdiff CL ;
   ods output LSMeanCL=logCI_reg;

   title3 'Quintiles of WHR';

 Proc Rank data=key out=key1 groups=5;
   var whr ;
    ranks Qwhr; 

 proc sort data=key1; by  Qwhr;
 proc means min max n data=key1; by  Qwhr;
   var whr;
 
 Proc freq data=key1;
   tables sex*Qwhr / nopercent nocol norow;

 Proc GLM data=key1;
   class sex Qwhr;
   model log_aer = Qwhr sex Qwhr*sex;

 ods graphics on;
 ods select lsmeans meanplot;  

 proc glimmix data=key1;
   class sex Qwhr;
   model log_aer = Qwhr sex Qwhr*sex;
   lsmeans Qwhr*sex / alpha=.32  plots=(meanplot(cl join sliceby=sex));
 run;
 ods graphics off;

 Proc GLM data=key1;
   class sex Qwhr;
   model log_aer = Qwhr sex;
   lsmeans sex / pdiff CL;
   ods output LSMeanCL=logCI_Q5;

 Proc GLM data=key1;
   where Qwhr in (1,2,3);
   class sex Qwhr;
   model log_aer = Qwhr sex;
   lsmeans sex / pdiff CL;
   ods output LSMeanCL=logCI_Q3;

 data results;
   set logCI_crude(in=crude) logCI_reg (in=reg) logCI_Q5 (in=Q5) logCI_Q3 (in=Q3);
   geometric_mean = exp(lsmean);
   left_CI=exp(LowerCL);
   right_CI = exp(UpperCL);
   model = "crude";
   if (reg=1) then model = "reg";
   if (Q5=1) then model="Q5";
   if (Q3=1) then model="Q3";

 proc print data=results;

 proc reg data=key;
   model log_AER=male;
 *                     Parameter       Standard
Variable     DF       Estimate          Error    t Value    Pr > |t|

Intercept     1        2.47319        0.05851      42.27      <.0001
male          1        0.31675        0.07771       4.08      <.0001
   ;
 proc reg data=key;
   model log_AER=male WHR;

 *                      Parameter       Standard
 Variable     DF       Estimate          Error    t Value    Pr > |t|
 Intercept     1        0.49503        0.51504       0.96      0.3367
 male          1        0.03768        0.10573       0.36      0.7216
 WHR           1        2.54856        0.65930       3.87      0.0001
 ;
 proc reg data=key;
   model WHR=male;
*                      Parameter       Standard
 Variable     DF       Estimate          Error    t Value    Pr > |t|
 Intercept     1        0.77619        0.00260     299.02      <.0001
 male          1        0.10950        0.00345      31.76      <.0001
 ;


 data mediation_test;
   a_b = 0.31675 - 0.03768;
   bc = 2.54856 * 0.10950;
   se = sqrt(2.54856*2.54856*0.00345*0.00345 + 0.10950*0.10950*0.65930*0.65930);
   z = bc/se;
   prob = 2.0*(1.0-probnorm(z));

 proc print;

 run; quit;