options ps=60 ls=80 pageno=1 formdlim='_';

data tomato;
input yield variety $ density;
cards;
10.5 a 10000
9.2 a 10000
7.9 a 10000
8.1 b 10000
8.6 b 10000
10.1 b 10000
16.1 c 10000
15.3 c 10000
17.5 c 10000
12.8 a 20000
11.2 a 20000
13.3 a 20000
12.7 b 20000
13.7 b 20000
11.5 b 20000
16.6 c 20000
19.2 c 20000
18.5 c 20000
12.1 a 30000
12.6 a 30000
14.0 a 30000
14.4 b 30000
15.4 b 30000
13.7 b 30000
20.8 c 30000
18.0 c 30000
21.0 c 30000
10.8 a 40000
9.1 a 40000
12.5 a 40000
11.3 b 40000
12.5 b 40000
14.5 b 40000
18.4 c 40000
18.9 c 40000
17.2 c 40000
;
run;

proc glm data=tomato;
class variety density;
model yield= variety density variety*density;
means variety density variety*density;
output out=avgs p=ybar;
run;

proc gplot data=avgs;
plot ybar*variety=density;
symbol1 v=triangle l=1 i=join cv=blue;
symbol2 v=circle l=1 i=join cv=orange;
symbol3 v=square l=1 i=join cv=purple;
symbol4 v=triangle l=1 i=join cv=cyan;
run;

proc glm;
class variety density;
model yield= variety density variety*density;
means variety density variety*density/tukey cldiff;
contrast "linear component of density effect" density -3 -1 1 3;
contrast "quadratic component of density effect" density 1 -1 -1 1;
contrast "cubic component of density effect" density -1 3 -3 1;
run;