options nocenter ps=50 ls=76;
data hardy;
   infile 'hardy.dat' firstobs=2;
   input potato regime temp code photo leak;

/* full model: you may want to consider various reduced models */
proc glm;
   class potato regime temp;
   model photo leak = potato | regime | temp / ss1 ss2 ss3;
   lsmeans potato | regime | temp / stderr pdiff out = lsm;
   output out = diag r = rphot rleak p = pphot pleak;

/* diagnostic plots */
proc plot data = diag;
   plot rphot * pphot = code;
   plot rleak * pleak = code;

/* simplistic interaction plot using marginal means */
proc sort data = hardy; by potato regime;
proc means noprint; by potato regime;
   var photo leak;
   output out = means mean = meanleak meanphot std = sdleak sdphot;
proc print data = means;
proc plot;
   plot meanleak * potato = regime;
   plot meanphot * potato = regime;

/* interaction plots using Least Squares Means */
proc print data = lsm;
proc sort data = lsm; by temp _NAME_;
proc plot; by temp _NAME_;
   plot lsmean * potato = regime;
run;