options ls=75 nodate;

data one;
   drop i;
   input a b c @;
   do i=1 to 3;
   input 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;
proc glm;
   class a b c;
   model y=a|b|c;
   lsmeans a*b*c/slice=(b*c);

estimate "mu[AB1C1]" a -1 1 a*b -1 0 1 0 a*c -1 0 0  1 0 0  
                                         a*b*c -1 0 0  0 0 0   1 0 0  0 0 0;
estimate "mu[AB1C2]" a -1 1 a*b -1 0 1 0 a*c 0 -1 0  0 1 0 
                                         a*b*c 0 -1 0  0 0 0   0 1 0  0 0 0;
estimate "mu[AB1C3]" a -1 1 a*b -1 0 1 0 a*c 0 0 -1  0 0 1 
                                         a*b*c 0 0 -1  0 0 0   0 0 1  0 0 0;

/* the contrasts below quantify temp-by-salinity interaction
seperately for each lvel of salinity */

   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;
   /*estimate "theta3: ABC3" a*b ? ? ? ?
                   a*b*c ? ? ? ? ? ? ? ? ? ? ? ?; */
   means a|b|c;
   output out=ymeans p=ymean;
run;

proc sort data=ymeans;
   by c;
run;

symbol1 i=join l=1 value=dot;
symbol2 i=join l=2 value=dot;

proc gplot data=ymeans;
   by c;
   plot ymean*b=a;
run;
data ymeans;
   set ymeans;
   abcombo=a||b;
run;
proc gplot data=ymeans;
   plot ymean*c=abcombo;
run;