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