/* SAS code from "Examples of Writing CONTRAST and ESTIMATE Statements */

/* EXAMPLE 1: A Two-Factor Model with Interaction */
data test;
   seed=6342454;
   do a=1 to 5;
      do b=1 to 2;
         do rep=1 to ceil(ranuni(seed)*5)+5;
            y=5 + a + b + a*b + rannor(seed);
            output;
         end;
      end;
   end;
run;

proc mixed data=test;
   class a b;
   model y=a b a*b / solution;
   lsmeans a*b / e;
run;


/* Computing the Cell Means Using the ESTIMATE Statement */
proc mixed data=test;
   class a b;
   model y=a b a*b;
   lsmeans a*b;
   estimate 'AB11' intercept 1 
                           a 1 0 0 0 0 
                           b 1 0
                         a*b 1 0 0 0 0 0 0 0 0 0;
   estimate 'AB12' intercept 1 
                           a 1 0 0 0 0 
                           b 0 1
                         a*b 0 1 0 0 0 0 0 0 0 0;
run;

/* Estimating and Testing a Difference of Means */
proc mixed data=test;
   class a b;
   model y= a b a*b;
   lsmeans a*b / diff;
   lsmestimate a*b 'AB11 - AB12' 1 -1 0 0 0 0 0 0 0 0;
   slice a*b / sliceby(a='1') diff;
   contrast 'AB11 - AB12' b 1 -1
                        a*b 1 -1 0 0 0 0 0 0 0 0;
   estimate 'AB11 - AB12' b 1 -1
                        a*b 1 -1 0 0 0 0 0 0 0 0;
run;

/*A More Complex Contrast*/
proc mixed data=test;
   class a b;
   model y= a b a*b;
   estimate 'avg AB11,AB12' intercept 1 
                                    a 1 
                                    b .5 .5
                                  a*b .5 .5;
   estimate 'avg AB21,AB22' intercept 1 
                                    a 0 1 
                                    b .5 .5
                                  a*b 0 0 .5 .5;
   contrast 'avg AB11,AB12 - avg AB21+AB22' a 1 -1
                                          a*b .5 .5 -.5 -.5;
   lsmestimate a*b 'avg AB11,AB12 - avg AB21+AB22' 1 1 -1 -1 / divisor=2;
run;

/* Comparing one interaction mean to the average of all 
   interaction means */
data test;
   seed=6342454;
   do a=1 to 2;
      do b=1 to 3;
         do rep=1 to ceil(ranuni(seed)*5)+5;
            y=5 + a + b + a*b + rannor(seed);
            output;
         end;
      end;
   end;
run;

proc mixed data=test;
    class a b;
    model y=a b a*b;
    estimate 'avg ABij' intercept 1 
                                a .5 .5 
                                b .333 .333 .333
                              a*b .167 .167 .167 .167 .167 .167;
run;

proc mixed data=test;
    class a b;
    model y=a b a*b;
    estimate 'AB12' intercept 1 
                            a 1 0 
                            b 0 1 0 
                          a*b 0 1 0 0 0 0;
    estimate 'avg ABij' intercept 6 
                                a 3 3 
                                b 2 2 2 
                              a*b 1 1 1 1 1 1 / divisor=6;
    estimate 'AB12 vs avg ABij' a 3 -3 
                                b -2 4 -2 
                              a*b -1 5 -1 -1 -1 -1 / divisor=6;
    lsmestimate a*b 'AB12 vs avg ABij' -1 5 -1 -1 -1 -1 / divisor=6;
run;


/* EXAMPLE 2: A Three Factor Model with Interactions */
data test;
   seed=8422636;
   do a=1 to 5;
      do b=1 to 2;
         do c=1 to 3;
            do rep=1 to ceil(ranuni(seed)*3)+3;
               y=5 + a + b + c + a*b + a*c + b*c + a*b*c + rannor(seed);
               output;
            end;
         end;
      end;
   end;
run;

proc mixed data=test;
   class a b c;
   model y=a|b|c / solution;
   lsmeans a*b*c;
run;

proc mixed data=test;
   class a b c;
   model y=a|b|c;
   estimate 'ABC121' intercept 1 a 1 b 0 1 c 1
                     a*b 0 1 a*c 1 b*c 0 0 0 1
                     a*b*c 0 0 0 1;
   estimate 'ABC212' intercept 1 a 0 1 b 1 c 0 1
                     a*b 0 0 1 a*c 0 0 0 0 1 b*c 0 1
                     a*b*c 0 0 0 0 0 0 0 1;
   contrast 'ABC121 - ABC212' a 1 -1 b -1 1 c 1 -1
                              a*b 0 1 -1 a*c 1 0 0 0 -1
                              b*c 0 -1 0 1
                              a*b*c 0 0 0 1 0 0 0 -1;
   estimate 'ABC121 - ABC212' a 1 -1 b -1 1 c 1 -1
                              a*b 0 1 -1 a*c 1 0 0 0 -1
                              b*c 0 -1 0 1
                              a*b*c 0 0 0 1 0 0 0 -1;
   lsmestimate a*b*c 'ABC121 - ABC212' 0 0 0 1 0 0 0 -1;
run;


/* EXAMPLE 3: A Two-Factor Logistic Model with Interaction Using Dummy
   and Effects Coding
*/
data uti;
   input diagnosis : $13. treatment $ response $ count @@;
   datalines;
complicated    A  cured 78  complicated   A not 28
complicated    B  cured 101 complicated   B not 11
complicated    C  cured 68  complicated   C not 46
uncomplicated  A  cured 40  uncomplicated A not 5
uncomplicated  B  cured 54  uncomplicated B not 5
uncomplicated  C  cured 34  uncomplicated C not 6
;

/* Dummy Coding */
/* Estimating and Testing Odds Ratios with Dummy Coding */
proc freq data=uti;
   where diagnosis = "complicated" and treatment in ("A","C");
   table treatment * response / relrisk norow nocol nopercent;
   weight count;
   run;

proc logistic data=uti;
   freq count;
   class diagnosis treatment / param=glm;
   model response(event='cured') = diagnosis treatment diagnosis*treatment;
   contrast 'trt A vs C in comp' treatment 1 0 -1
            diagnosis*treatment 1 0 -1 0 0 0 / estimate=both;
   output out=out xbeta=xbeta;
   oddsratio treatment / at(diagnosis='complicated');
   lsmeans diagnosis*treatment / ilink exp diff;
   lsmestimate diagnosis*treatment 'A vs C complicated' 1 0 -1 / exp;
   slice diagnosis*treatment / sliceby(diagnosis='complicated') diff exp;
   run;

proc print data=out noobs;
   where diagnosis="complicated" and response="cured" and
         treatment in ("A","C");
   var diagnosis treatment xbeta;
   run;

proc genmod data=uti;
   freq count;
   class diagnosis treatment;
   model response = diagnosis treatment diagnosis*treatment / dist=binomial;
   estimate 'trt A vs C in comp' treatment 1 0 -1
            diagnosis*treatment 1 0 -1 0 0 0 / exp;
   output out=out xbeta=xbeta;
   run;

/* A Nested Model */
proc logistic data=uti;
   freq count;
   class diagnosis treatment / param=glm;
   model response(event='cured') = diagnosis treatment(diagnosis) / expb;
   contrast 'trt A vs C in comp'
            treatment(diagnosis) 1 0 -1 0 0 0 / estimate=both;
   run;

proc genmod data=uti;
   freq count;
   class diagnosis treatment;
   model response = diagnosis treatment(diagnosis) / dist=binomial;
   estimate 'trt A vs C in comp'
            treatment(diagnosis) 1 0 -1 0 0 0 / exp;
   run;

/* Effects Coding */
/* Estimating and Testing Odds ratios with Effects Coding */
proc logistic data=uti;
   freq count;
   class diagnosis treatment;
   model response(event='cured') = diagnosis treatment diagnosis*treatment;
   contrast 'trt A vs C in comp' treatment 2 1
            diagnosis*treatment 2 1 / estimate=both;
   run;

proc genmod data=uti;
   freq count;
   class diagnosis treatment / param=effect;
   model response = diagnosis treatment diagnosis*treatment / dist=binomial;
   estimate 'trt A vs C in comp' treatment 2 1
            diagnosis*treatment 2 1 / exp;
   run;

/* A More Complex Contrast with Effects Coding */
proc logistic data=uti;
   freq count;
   class diagnosis treatment;
   model response(event='cured') = diagnosis treatment diagnosis*treatment;
   contrast 'trt 1 vs avg trt in comp' treatment 1 0
            diagnosis*treatment 1 0 / estimate=both;
   run;

proc genmod data=uti;
   freq count;
   class diagnosis treatment / param=effect;
   model response = diagnosis treatment diagnosis*treatment / dist=binomial;
   estimate 'trt A vs avg trt in comp' treatment 1 0
            diagnosis*treatment 1 0 / exp;
   run;

proc catmod data=uti;
   weight count;
   model response = diagnosis treatment diagnosis*treatment;
   contrast 'trt 1 vs avg trt in comp' treatment 1 0
            diagnosis*treatment 1 0 / estimate=both;
   run;

proc catmod data=uti;
   weight count;
   model response = diagnosis 
                    treatment(diagnosis='complicated')
                    treatment(diagnosis='uncomplicated');
   run;


/* Example 4: Comparing Nested Models */
data detergent;
   input Softness $ Brand $ Previous $ Temperature $ Count @@;
   datalines;
soft X yes high 19   soft X yes low 57
soft X no  high 29   soft X no  low 63
soft M yes high 29   soft M yes low 49
soft M no  high 27   soft M no  low 53
med  X yes high 23   med  X yes low 47
med  X no  high 33   med  X no  low 66
med  M yes high 47   med  M yes low 55
med  M no  high 23   med  M no  low 50
hard X yes high 24   hard X yes low 37
hard X no  high 42   hard X no  low 68
hard M yes high 43   hard M yes low 52
hard M no  high 30   hard M no  low 42
;

ods select modelfit type3;
ods output modelfit=full;
proc genmod data=detergent;
   class Softness Previous Temperature;
   freq Count;
   model Brand = Softness|Previous|Temperature / dist=binomial type3;
run;

ods select modelfit type3;
ods output modelfit=reduced;
proc genmod data=detergent;
   class Softness Previous Temperature;
   freq Count;
   model Brand = Softness Previous Temperature / dist=binomial type3;
run;

/* Writing a contrast to compare models */
/* Using indicator (dummy) coding (not shown in note) */
proc genmod data=detergent;
   class Softness Previous Temperature;
   freq Count;
   model Brand = Softness|Previous|Temperature / dist=binomial;
   contrast 'lrt'

   softness*previous 1 -1 0 0 -1 1 
   softness*previous*temperature 0 1 0 -1 0 0 0 0 0 -1 0 1,
   
   softness*previous 0 0 1 -1 -1 1 
   softness*previous*temperature 0 0 0 0 0 1 0 -1 0 -1 0 1,
   
   softness*temperature 1 -1 0 0 -1 1 
   softness*previous*temperature 0 0 1 -1 0 0 0 0 0 0 -1 1,
   
   softness*temperature 0 0 1 -1 -1 1 
   softness*previous*temperature 0 0 0 0 0 0 1 -1 0 0 -1 1,
   
   previous*temperature 1 -1 -1 1 
   softness*previous*temperature 0 0 0 0 0 0 0 0 1 -1 -1 1,
   
   softness*previous*temperature 1 -1 -1 1 0 0 0 0 -1 1 1 -1,
   softness*previous*temperature 0 0 0 0 1 -1 -1 1 -1 1 1 -1;
run;

/* Effects coding -- LR test */
ods select contrasts;
proc genmod data=detergent;
   class Softness Previous Temperature / param=effect;
   freq Count;
   model Brand = Softness|Previous|Temperature / dist=binomial;
   contrast 'lrt'
   softness*previous 1 0,
   softness*previous 0 1,
   softness*temperature 1 0,
   softness*temperature 0 1,
   previous*temperature 1,
   softness*previous*temperature 1 0,
   softness*previous*temperature 0 1;
run;

/* Effects coding -- Wald test */
ods select contrasttest;
proc logistic data=detergent;
   class Softness Previous Temperature / param=effect;
   freq Count;
   model Brand = Softness|Previous|Temperature;
   contrast 'lrt'
   softness*previous 1 0,
   softness*previous 0 1,
   softness*temperature 1 0,
   softness*temperature 0 1,
   previous*temperature 1,
   softness*previous*temperature 1 0,
   softness*previous*temperature 0 1;
run;