data one;
  input surface filler p y;
  cards;
1 1 25 194
1 1 25 208
1 1 50 233
1 1 50 241
1 1 75 265
1 1 75 269
1 2 25 239
1 2 25 187
1 2 50 224
1 2 50 243
1 2 75 243
1 2 75 226
2 1 25 155
2 1 25 173
2 1 50 198
2 1 50 177
2 1 75 235
2 1 75 229
2 2 25 137
2 2 25 160
2 2 50 129
2 2 50 98
2 2 75 155
2 2 75 132
;

proc glm;
  class surface filler p;
  model y=surface filler p surface*filler surface*p filler*p surface*filler*p;
  lsmeans surface filler p surface*filler surface*p filler*p surface*filler*p / stderr;
  ods output LSMeans=means;
run;

proc glm data=one;
  class surface filler;
  model y=surface filler p p*p surface*filler surface*p surface*p*p filler*p filler*p*p
          surface*filler*p surface*filler*p*p;
run;

proc glm data=one;
  class surface filler;
  model y=surface filler surface*filler p filler*p / solution;
  estimate 'intercept for S1 F1' intercept 1 surface 1 0 filler 1 0 surface*filler 1 0 0 0;
  estimate 'intercept for S1 F2' intercept 1 surface 1 0 filler 0 1 surface*filler 0 1 0 0;
  estimate 'intercept for S2 F1' intercept 1 surface 0 1 filler 1 0 surface*filler 0 0 1 0;
  estimate 'intercept for S2 F2' intercept 1 surface 0 1 filler 0 1 surface*filler 0 0 0 1;
  estimate 'slope for S1 F1' p 1 p*filler 1 0;
  estimate 'slope for S1 F2' p 1 p*filler 0 1;
  estimate 'slope for S2 F1' p 1 p*filler 1 0;
  estimate 'slope for S2 F2' p 1 p*filler 0 1;
  estimate '(S2 F1 p=25)-(S2 F2 p=25)' filler 1 -1 surface*filler 0 0 1 -1 p*filler 25 -25;
run;