/* 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;

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;

proc mixed data=test;
   class a b;
   model y= a b a*b;
   lsmeans a*b / pdiff;
   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;

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;
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 glm 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;
quit;

ods select estimates;
proc glm 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;
run;
quit;

/* 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 / solution;
   lsmeans a*b*c / pdiff;
   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;
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 = 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;
   run;

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

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

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

/* A More Complex Contrast with Effects Coding */
proc logistic data=uti;
   freq count;
   class diagnosis treatment;
   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*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 Models -- Likelihood Ratio Test */
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;

/* Constructing a likelihood ratio test */
data lrt;
   lr=2*(-682.2478 - -686.3618);
   df=7;
   p=1-probchi(lr,df);
   run;
proc print noobs;
   format p pvalue.;
   run;

data lrt;
  retain dff dfr LRDF;
  set full end=endf;
  dff=df;
  if endf then llf=value;
  set reduced end=endr;
  dfr=df;
  if endr then llr=value;
  if _n_=1 then LRDF=dfr-dff;
  if endf and endr then do;
    LR=2*(llf-llr);
    p=1-probchi(LR, LRDF);
    keep LR LRDF p;
    output;
  end;
  run;

/* Writing a contrast to compare models */
/* Indicator (dummy) coding */
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;