/* pain2.sas */
title 'Covariance structure repeated measures analysis of the pain data';

/* Read data from an old .xls spreadsheet */
proc import datafile="/home/u1407221/441s24/data/Pain.xls" 
            out=ouch0 dbms=xls replace;
            getnames=yes;
proc print;

/*           Dose 1         Dose 2         Dose 3
             ------         -----          ------
Drug 1     Drug1Dose1     Drug1Dose2     Drug1Dose3
Drug 2     Drug2Dose1     Drug2Dose2     Drug2Dose3   */

data hurt;
     set ouch0;
     /* Need each response on a separate line for proc mixed. */
     array response{6} Drug1Dose1 -- Drug2Dose3;
     index = 0;
     do drug = 1 to 2;
         do dose = 1 to 3;
             index = index + 1;
             pain = response{index};
             output;
         end;
     end;
     keep patient drug dose pain;

proc print data = hurt (obs=12); /* Just the first 12 lines */

proc mixed data=hurt;
     title2 'Unknown covariance structure';
     class drug dose;
     model pain = drug|dose;
     repeated / type=un subject=patient;
     /* Pairwise comparisons of marginal means for Dose */
     lsmeans dose / pdiff tdiff adjust = bon; 
run;

proc mixed data=hurt;
     title2 'Compound symmetry';
     class drug dose;
     model pain = drug|dose;
     repeated / type=cs subject=patient;
     lsmeans dose / pdiff tdiff adjust = bon;
run;

/* F-tests from cs do not match up with the proc glm univariate tests
   in pain1.sas that are produced as a by-product of the multivariate
   approach. If you know how to read it, output from the code below 
   shows that for this (simple and quite common) example, exact classical
   F-tests for the main effects do not exist unless the interaction is 
   assumed to be absent. The proc mixed tests match what you get by making 
   this assumption. It seems that the "by-product" univariate repeated
   measures analysis produced by proc glm is doing something else. I do
   not know exactly what it is doing, except that in every example I can 
   find where all the classical tests DO exist, proc glm matches them.
   If we want the compound symmetry assumption, we will use proc mixed. */

proc glm data=hurt noprint; /* Not showing the output */
     title2 'Trying classical the old fashioned way';
     class patient drug dose;
     model pain = patient drug|dose;
     random patient / test;
run;
quit;