options ls=75 nodate;

data one;
   drop i;
   input variety b @;
   density=10*b;
   if variety=3 then delete;
   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 glm;
   class variety density;
   model y=variety|density;
   contrast 'linear effect' density -3 -1 1 3;
   contrast 'quadratic effect' density 1 -1 -1 1;
   estimate "linear opc variety 1" density -3 -1 1 3 variety*density -3 -1 1 3;
   estimate "linear opc variety 2" density -3 -1 1 3 variety*density 0 0 0 0 -3 -1 1 3;
   estimate "linear X variety" variety*density 3 1 -1 -3 -3 -1 1 3;
   estimate "quadratic opc variety 1" density 1 -1 -1 1 variety*density 1 -1 -1 1;
   estimate "quadratic opc variety 2" density 1 -1 -1 1 variety*density 0 0 0 0 1 -1 -1 1;
   estimate "quadratic X variety" variety*density -1 1 1 -1 1 -1 -1 1;
   contrast "linear X variety" variety*density 3 1 -1 -3 -3 -1 1 3;
   contrast "quadratic X variety" variety*density -1 1 1 -1 1 -1 -1 1;
   output out=two p=p lclm=lcl uclm=ucl;
run;

proc glm;
   class variety;
   model y=variety density density*density variety*density variety*density*density/ss1 solution;
   *model y=variety density density*density; 
run;
goptions colors=(black);

symbol1 i=join value=dot; symbol2 i=join value=plus;
symbol3 i=join value=diamond; symbol2 i=join value=star;

proc gplot data=two;
   plot p*density=variety;
run;