*******************
* bloodpress1.sas *
*******************;

options ls=78 pageno=1 nodate ;

data blood; 
  infile 'C:\Users\Dan Hall\Documents\Dans Work Stuff\Courses\STAT8630\Spring13\bloodpress.dat';
  input drug $	exercise $	person	IBP	Time	BP;
run;

proc sort data=blood; 
  by exercise drug time;
run;
TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Plots of blood pressure versus initial blood pressure';
proc plot data=blood;
  where exercise="No" and drug="No" and time=1; *comment this line out to see them all;
  by exercise drug time;
  plot bp*ibp;
run;
 
ods listing close; *This suppresses all of the output from the following calls to
                    PROC MIXED that fit various different var-cov models;
TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 1a: Split-plot type model with linear covariate.';
title3 'Different intercept and slope for each trt * time combo';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr; *cell means type parameterization;
  /* would have been equivalent to use the following model statement: */
  *Model BP= IBP|Drug|Exercise|Time/ddfm=kr;
  Random person(Drug*Exercise);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit1a;
RUN;
title3;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 1b: Same as model 1a but compound symmetry specified on REPEATED statement';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  Repeated time/type=cs subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit1b;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 2: Same mean structure as model 1, heteroscedastic CS var-cov structure';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  Repeated time/type=csh subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit2;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 3: Same mean structure as model 1, ar(1) structure for R plus random person effects';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=ar(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit3;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4: Same mean structure as model 1, heteroscedastic ar(1) structure for R plus random person effects';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 5: Same mean structure as model 1, antedependence of order 1 structure for R, no random person effects';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  *Random person(Drug*Exercise);
  Repeated time/type=ante(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit5;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 6: Same mean structure as model 1, Toeplitz structure for R';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  Repeated time/type=toep subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit6;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 7: Same mean structure as model 1, heteroscedastic Toeplitz structure for R';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  Repeated time/type=toeph subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit7;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 8: Same mean structure as model 1, unstructured R matrix';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug*Exercise*Time  IBP*Drug*Exercise*Time/noint ddfm=kr;
  Repeated time/type=un subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit8;
RUN;
 
TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4a: ARH(1) + random subject effects wins (by AIC). Now refit model 4 to simplify mean structure';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= IBP|Drug|Exercise|Time/ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4a;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4b: Same var-cov structure as in model 4a. Remove ibp*drug*exercise*time from mean.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP IBP*drug ibp*exercise ibp*time ibp*drug*exercise ibp*drug*time 
            ibp*exercise*time/ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4b;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4c: Same var-cov structure as in model 4a. Remove ibp*exercise*time from mean.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP IBP*drug ibp*exercise ibp*time ibp*drug*exercise ibp*drug*time 
            /ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4c;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4d: Same var-cov structure as in model 4a. Remove ibp*exercise*drug from mean.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP IBP*drug ibp*exercise ibp*time ibp*drug*time 
            /ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4d;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4e: Same var-cov structure as in model 4a. Remove ibp*time*drug from mean.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP IBP*drug ibp*exercise ibp*time /ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4e;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4f: Same var-cov structure as in model 4a. Remove ibp*exercise from mean.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP IBP*drug ibp*time /ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4f;
RUN;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4g: Same var-cov structure as in model 4a. Remove ibp*drug from mean.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP ibp*time /ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  ods exclude modelinfo dimensions nobs iterhistory lrt classlevels;
  ods output fitstatistics=fit4g;
RUN;
ods listing; * Re-activate the output;

data aicbic;
  set fit1a;
  model="1a";
  output;
  set fit1b;
  model="1b";
  output;
  set fit2;
  model="2";
  output;
  set fit3;
  model="3";
  output;
  set fit4;
  model="4";
  output;
  set fit4a;
  model="4a";
  output;
  set fit4b;
  model="4b";
  output;
  set fit4c;
  model="4c";
  output;
  set fit4d;
  model="4d";
  output;
  set fit4e;
  model="4e";
  output;
  set fit4f;
  model="4f";
  output;
  set fit4g;
  model="4g";
  output;
run;
proc sort data=aicbic;
  by model;
run;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Fit statistics for model 1a-4g';
proc print data=aicbic;
run;
TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Summary statistics for initial blood pressure';
proc means data=blood q1 median q3 mean std;
  var ibp;
run;


TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4g: Model 4g is final model. Now refit for additional analyses.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP ibp*time /ddfm=kr s;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  lsmeans exercise*drug*time / slice=drug*time at ibp=136.8 ; * at mean of ibp - 1 sd;
  lsmeans exercise*drug*time / slice=drug*time at ibp=143.8 ; * at mean of ibp;
  lsmeans exercise*drug*time / slice=drug*time at ibp=150.7 ; * at mean of ibp + 1 sd;
  *lsmeans exercise*drug*time / slice=exercise*time at ibp=136.8 ; * at mean of ibp - 1 sd;
  *lsmeans exercise*drug*time / slice=exercise*time at ibp=143.8 ; * at mean of ibp;
  *lsmeans exercise*drug*time / slice=exercise*time at ibp=150.7 ; * at mean of ibp + 1 sd;
  ods output lsmeans=lsm;
  ods  exclude Lsmeans modelinfo dimensions nobs iterhistory lrt classlevels;

RUN;

data mns;
  set lsm;
  if (drug='Yes' and exercise='Yes') then drugex='YY';
  else if (drug='Yes' and exercise='No') then drugex='YN';
  else if (drug='No' and exercise='Yes') then drugex='NY';
  else if (drug='No' and exercise='No') then drugex='NN';
run;

goptions reset=all;
symbol1 value=circle color=red line=1 interpol=join;
symbol2 value=square color=blue line=2 interpol=join;
symbol3 value=triangle color=green line=3 interpol=join;
symbol4 value=diamond color=black line=4 interpol=join;
TITLE 'Blood Pressure Example - RM-ANCOVA, Results from model 4g';
title2 'Profile plot of mean BP over time for each group at 3 levels of IBP';
proc gplot data=mns;
  by ibp;
  plot estimate*time=drugex;
run; 
quit;


data blood2;
  set blood;
  time1=(time=1);
  time2=(time=2);
  time3=(time=3);
  time4=(time=4);
  time5=(time=5);
run;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4g alt parmztion. Comparisons with control at each timept and at avg ibp';
Proc Mixed  DATA=blood2;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|time1 Drug|Exercise|time2 Drug|Exercise|time3 Drug|Exercise|time4 Drug|Exercise|time5 IBP 
  ibp*time1 ibp*time2 ibp*time3 ibp*time4 ibp*time5  /ddfm=kr ;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  lsmeans exercise*drug / diff=control('No' 'No') at (ibp time1 time2 time3 time4 time5)=(143.8 1 0 0 0 0) adjust=dunnett; * at mean of ibp and time 1;
  lsmeans exercise*drug / diff=control('No' 'No') at (ibp time1 time2 time3 time4 time5)=(143.8 0 1 0 0 0) adjust=dunnett; * at mean of ibp and time 2;
  lsmeans exercise*drug / diff=control('No' 'No') at (ibp time1 time2 time3 time4 time5)=(143.8 0 0 1 0 0) adjust=dunnett; * at mean of ibp and time 3;
  lsmeans exercise*drug / diff=control('No' 'No') at (ibp time1 time2 time3 time4 time5)=(143.8 0 0 0 1 0) adjust=dunnett; * at mean of ibp and time 4;
  lsmeans exercise*drug / diff=control('No' 'No') at (ibp time1 time2 time3 time4 time5)=(143.8 0 0 0 0 1) adjust=dunnett; * at mean of ibp and time 5;
  lsmeans exercise*drug / diff=control('No' 'No') at (ibp time1 time2 time3 time4 time5)=(143.8 0 0 0 0 0) adjust=dunnett; * at mean of ibp and time 6;
  ods  exclude Lsmeans modelinfo dimensions nobs iterhistory lrt classlevels tests3;
RUN;



/*
Alternatively, one could follow the approach used by Milliken
and Johnson and perform all pairwise comparisons
between the joint means (the means for each exercise by drug by time combination).
Given the significance effects of the covariate, this should be done
at slected values of the covariate. This generates a very large number
of pairwise comparisons, so M&J used the "simulate" option to 
control the maximum experimentwise type i error rate. Then they only
printed out the comparisons between exercise levels for each combination 
of drug and time, and comparisons between drug levels for each combination of
exercise and time. I find this approach overkill - they are making too many 
comparisons in my view, and conducting a bit of a "fishing expedition"
for significant differences. I prefer the pairwise comparisons with
the control approach that I implemented above, but what comparisons/contrasts
one does should really be based upon the questions of scientific interest 
(preferably the a priori questions of interest), which is not a statistical
question as much as it is one to be answered by the subject-matter researchers
who collected the data.

Below is the (somewhat modified) code used by Milliken and Johnson in
their presentation of this example.



TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Model 4g: Model 4g is final model. Now refit for additional analyses.';
Proc Mixed  DATA=blood;
  Class Drug Exercise person Time;
  Model BP= Drug|Exercise|Time IBP ibp*time /ddfm=kr;
  Random person(Drug*Exercise);
  Repeated time/type=arh(1) subject=person(exercise*drug);
  lsmeans exercise*drug*time / diff at ibp=136.8 adjust=simulate; * at mean of ibp - 1 sd;
  lsmeans exercise*drug*time / diff at ibp=143.8 adjust=simulate; * at mean of ibp;
  lsmeans exercise*drug*time / diff at ibp=150.7 adjust=simulate; * at mean of ibp + 1 sd;
  ods output lsmeans=lsm;
  ods listing exclude Lsmeans;
  ods output Diffs=dif;
  ods listing exclude Diffs;
RUN;

proc sort data=lsm; 
  by ibp exercise drug time;
run;

TITLE 'Blood Pressure Example - RM-ANCOVA';
title2 'Results from model 4g';
title3 'Adjusted Means';
proc print data=lsm;
  by ibp ; 
  var exercise drug time estimate stderr df tvalue probt;
run;
title3;

Proc sort data=dif; 
  by  drug Time ibp;
run;

data exer; 
  set dif;
  if drug=_drug and Time=_Time;
  keep drug Time ibp exercise _exercise estimate stderr df tvalue adjp;
run;
TITLE 'Blood Pressure Example - RM-ANCOVA, Results from model 4g';
title2 'Comparisons of levels of Exercise at each level of drug x time and ibp';
Proc print;
  by drug time ibp; 
  where (ibp=143.8);
  var    exercise _exercise estimate stderr df tvalue adjp;
run;

Proc sort data=dif; 
  by  exercise Time ibp;
run;

data drug; 
  set dif;
  if exercise=_exercise and Time=_Time;
  keep exercise  Time ibp drug _drug estimate stderr df tvalue adjp;
run;

TITLE 'Blood Pressure Example - RM-ANCOVA, Results from model 4g';
title2 'Comparisons of levels of Drug at each level of exercise x time and ibp';
Proc print; 
  by exercise time ibp;
  where (ibp=143.8);
  var drug _drug estimate stderr df tvalue adjp ;
run;
*/