options ls=75 nodate;

data one;
   drop i;
   input a b c @;
   do i=1 to 3;
   input y @;
   y0=sqrt(y);
   output;
   end;
   cards;
   25 80 10 86 52 73
   25 80 25 544 371 482
   25 80 40 390 290 397 
   25 160 10 53 73 86
   25 160 25 393 398 208
   25 160 40 249 265 243
   35 80 10 439 436 349
   35 80 25 249 245 330
   35 80 40 247 277 205
   35 160 10 324 305 364
   35 160 25 352 267 316
   35 160 40 188 223 281
;
run;

symbol1 i=join value=dot;
symbol2 i=join value=plus;

proc sort;  by b; run;

proc means noprint nway;  class a b c;  var y; output out=two mean=ymean;run;

proc sort;  by b; run;
proc gplot data=two;
   by b;
   plot ymean*a=c;
run;
proc glm data=one;
   class a b c;
   model y=a|b|c;
   contrast "theta1: ABC1" a*b -1 1 1 -1
                   a*b*c -1 0 0 1 0 0 1 0 0 -1 0 0;
   estimate "theta1: ABC1" a*b -1 1 1 -1
                   a*b*c -1 0 0 1 0 0 1 0 0 -1 0 0;
   contrast "theta2: ABC2" a*b -1 1 1 -1
                   a*b*c 0 -1 0 0 1 0 0 1 0 0 -1 0;
   estimate "theta2: ABC2" a*b -1 1 1 -1
                   a*b*c 0 -1 0 0 1 0 0 1 0 0 -1 0;
   contrast "theta3: ABC3" a*b -1 1 1 -1
                   a*b*c 0 0 -1 0 0 1 0 0 1 0 0 -1;
   estimate "theta3: ABC3" a*b -1 1 1 -1
                   a*b*c 0 0 -1 0 0 1 0 0 1 0 0 -1;
   contrast "t1-av(t2+t3)" a*b*c -2 1 1 2 -1 -1 2 -1 -1 -2 1 1;
   estimate "t1-av(t2+t3)" a*b*c -2 1 1 2 -1 -1 2 -1 -1 -2 1 1/divisor=2;
   estimate "t2-av(t1+t3)" a*b*c 1 -2 1 -1 2 -1 -1 2 -1 1 -2 1/divisor=2;
   *contrast "AC1"  a*c 1 -1 0 -1 1 0;
   *estimate "AC1"  a*c 1 -1 0 -1 1 0;
   *contrast "AC2"  a*c 1 1 -2 -1 -1 2;
   *estimate "AC2"  a*c 1 1 -2 -1 -1 2;
   *contrast "A1BC" b -1 1 
                   a*b 1 -1 0 0
                   b*c 0 0 1 0 0 -1
                   a*b*c 0 0 1 0 0 -1 0 0 0 0 0 0;
run;