data one;
  input density yield;
  linear=density;
  quadratic=density**2;
  cubic=density**3;
  quartic=density**4;
  lack_of_fit=density;
  cards;
10 12.2
10 11.4
10 12.4
20 16.0
20 15.5
20 16.5
30 18.6
30 20.2
30 18.2
40 17.6
40 19.3
40 17.1
50 18.0
50 16.4
50 16.6
;

proc glm;
  class density;
  model yield=density;
run;

proc glm;
  class lack_of_fit;
  model yield=linear lack_of_fit;
run;

proc glm;
  class lack_of_fit;
  model yield=linear quadratic lack_of_fit;
run;

proc glm;
  class lack_of_fit;
  model yield=linear quadratic cubic lack_of_fit;
run;

proc glm;
  class lack_of_fit;
  model yield=linear quadratic cubic quartic lack_of_fit;
run;

proc glm;
  model yield=density density*density / solution clparm;
  estimate 'Mean Yield at Density 25' intercept 1 density 25 density*density 625;
run;

proc glm;
  model yield=density density*density;
  output out=two
    predicted=yhat
    stdp=stderrm
    lclm=l95mean
    uclm=u95mean
    / alpha=0.05;
run;

proc sort data=two out=two;
  by density;
run;

proc means data=two noprint;
  var yield yhat stderrm l95mean u95mean;
  by density;
  output out=tablep91 mean=yield yhat stderrm l95mean u95mean;
run;

proc print data=tablep91;
  var density yield yhat stderrm l95mean u95mean;
run;