title 'Cortisol Data';
data cort; 
  infile '.\Dropbox\Teaching\1301-StatisticalPackages\datasets\cortisol.dat';
input group cortisol;
label group = 'Psychiatric status';
run;

proc format;
   value group 1 = 'Normal'
               2 = 'Major Depression'
			   3 = 'Bipolar Depression'
			   4 = 'Schizophrenia'
			   5 = 'Atypical';

run;

proc tabulate data = cort; 
  class group;
  var cortisol;
  tables group, cortisol*(mean std n);
  format group group.;
run; 

proc sort data = cort; 
 by group;
proc boxplot data = cort;
  plot cortisol * group; 
  format group group.;
run;

data cort; 
 set cort; 
 lc = log(cortisol);
 label lc = 'Log Cortisol Level';
run;
proc boxplot data = cort;
  plot lc * group; 
  format group group.;
run;
 
proc anova data = cort;
  class group; 
  model lc = group; 
  means group / scheffe;
  format group group.;
run;

proc anova data = cort;
  class group; 
  model lc = group; 
  means group/hovtest tukey scheffe  cldiff;
run;

proc glm data = cort order = data;
  class group; 
  model lc = group; 
  contrast 'Normal vs ave. Disorder' group -1 0.25 0.25 0.25 0.25;
  contrast 'Normal vs ave. Disorder2' group -4 1 1 1 1;
  contrast 'Normal vs Bipolar' group 2 0 -1 0 0 ;
  format group group.;
run;


title 'Treating Hypertension';


data hyper;
  infile 'C:\Users\sungkyu\Dropbox\Teaching\1301-StatisticalPackages\datasets\hypertension.dat';
  input n1-n12;
  if _n_<4 then biofeed='P';
           else biofeed='A';
  if _n_ in(1,4) then drug='X';
  if _n_ in(2,5) then drug='Y';
  if _n_ in(3,6) then drug='Z';
  array nall {12} n1-n12;
  do i=1 to 12;
      if i>6 then diet='Y';
                 else diet='N';
          bp=nall{i};
          cell=drug||biofeed||diet;
          output;
  end;
  drop i n1-n12;
run;

proc print data=hyper;
 format diet $YN.;
run;

proc format;
   value $YN 'Y' = 'Yes'
             'N' = 'No';
   value $PA 'P' = 'Present'
             'A' = 'Absent';
run;
proc print data=hyper;
 format diet $YN. biofeed $PA.;
run;

 
proc tabulate data=hyper;
  class  drug diet biofeed;
  var bp;
  table drug*diet*biofeed,
        bp*(mean std n);
  format diet $YN. biofeed $PA.;
run;

proc anova data=hyper;
  class cell;
  model bp=cell;
  means cell / hovtest;
  run;

title2 'Two-way models';
proc anova data = hyper; 
  class diet drug biofeed;
  model bp = diet drug diet*drug;
  format diet $YN. biofeed $PA.; 
run; 
  means diet*drug;
  ods output means = twowayDIET_DRUG;
run;

proc anova data = hyper; 
  class diet drug biofeed;
  model bp = drug|biofeed;
  format diet $YN. biofeed $PA.;
run; 
  means biofeed*drug;
  ods output means = twowaybiofeed_DRUG;
run;
proc print data=twowaybiofeed_DRUG;
run;
proc sgplot data=twowaybiofeed_DRUG;
  series y=mean_bp x=drug / group=biofeed;
run;

proc sgplot data=twowayDIET_DRUG;
  series y=mean_bp x=drug / group=diet;
run;

proc anova data=hyper;
  class diet drug biofeed;
  model bp=diet|drug|biofeed;
  format diet $YN. biofeed $PA.;
run;
  means diet*drug*biofeed;
  ods output means=outmeans;
run; 

proc sgpanel data=outmeans;
  panelby drug / rows = 1 ;
  series y=mean_bp x=biofeed / group=diet;
run; 

proc anova data=hyper;
  class diet drug biofeed;
  model bp=diet drug biofeed;
  means drug / scheffe;
run;