options ls=75 nodate;

data one;
   drop i;
   input a b @;
   b=10*b;
   bs=b*b;
   do i=1 to 3;
      input y @;
	  output;
   end;
   cards;
1 1 7.9 9.2 10.5
2 1 8.1 8.6 10.1
3 1 15.3 16.1 17.5
1 2 11.2 12.8 13.3
2 2 11.5 12.7 13.7
3 2 16.6 18.5 19.2
1 3 12.1 12.6 14.0
2 3 13.7 14.4 15.4
3 3 18.0 20.8 21.0
1 4 9.1 10.8 12.5
2 4 11.3 12.5 14.5
3 4 17.2 18.4 18.9
;
run;
proc print data=one;run;

proc glm;
   class b a ;
   *model y=a b b*b/solution;
   model y=a|b;
   contrast 'linear effect' b -3 -1 1 3;
   contrast 'quadratic effect' b 1 -1 -1 1;
   contrast 'cubic effect' b -1 3 -3 1;
   means a b / tukey;
   means a*b;
   lsmeans a*b/slicediff=(a b);
   output out=two p=p r=r;
run;
/*symbol1 i=rq value=dot; symbol2 i=rq value=plus;
symbol3 i=rq value=diamond; symbol2 i=rq value=star;*/
symbol value=dot;

proc gplot;
   plot p*b=a;
run;

symbol i=join value=dot r=4;
proc gplot;
   plot y*b=a;
run;

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

symbol value=dot i=join;
proc gplot;
   title "interaction plot";
   plot ymean*b=a;
run;


symbol value=dot i=rq;
proc gplot;
   title "interaction plot";
   plot ymean*b=a;
run;


proc glm data=one;
   class a;
   model y=a b bs;
   output out=quadfit p=p r=r;
   estimate "v1 - v2         " a 1 -1;
   estimate "v1 - v3" a 1 0 -1;
   estimate "v3 - v2" a 0 -1 1;
   estimate "v3-mean of v1,v2" a -1 -1 2/divisor=2;
   contrast "v3-mean of v1,v2" a -1 -1 2;
run;
symbol value=dot i=rq r=4;
proc gplot data=quadfit;
   plot p*b=a;
run;