options ls=79 ps=50 nocenter;

data tukey;
   infile 'tukey.dat' firstobs=2;
   input temp humid y;

/* fit additive model */
/* saving humid margins as lsmean in lshumid file */
proc glm;
   class temp humid;
   model y = temp humid;
   lsmeans humid / stderr out=lshumid;
   lsmeans temp / stderr out=lstemp;
   output out=add p=pred;

/* fit tukey interaction model */
data inter; set add;
   inter = pred * pred / 2;
/* only examine the test for "inter" covariate */
proc glm;
   class temp humid;
   model y = temp humid inter / ss1 solution;
   output out=tukpred p=ptukey;

/* Tukey Interaction Plots */
proc plot;
   plot y*humid=temp ptukey*humid='*' / overlay;

/* Tukey Margin Plots -- against humid */
proc sort data=tukpred; by humid;
proc sort data=lshumid; by humid;
data tukplot; merge tukpred lshumid; by humid;
proc plot;
   plot y*lsmean=temp;		/* actual values */
   plot ptukey*lsmean=temp;	/* Tukey estimates */
/* Tukey Margin Plots -- against temp */
proc sort data=tukpred; by temp;
proc sort data=lstemp; by temp;
data tukplot; merge tukpred lstemp; by temp;
proc plot;
   plot y*lsmean=humid;		/* actual values */
   plot ptukey*lsmean=humid;	/* Tukey estimates */

/* Fit Two Mandel Interaction Models */
/* check only the humid*pred interaction */
proc glm data=inter;
   class temp humid;
   model y = temp humid humid*pred;
   output out=mand1 p=pred1;
/* check only the temp*pred interaction */
proc glm data=inter;
   class temp humid;
   model y = temp humid temp*pred;
   output out=mand2 p=pred2;

/* Mandel Margin Plots -- against humid */
proc sort data=mand1; by humid;
proc sort data=mand2; by humid;
proc sort data=lshumid; by humid;
data mandplot; merge mand1 mand2 lshumid; by humid;
proc plot;
   plot pred1*lsmean=temp;	/* Mandel humid*pred Model */
   plot pred2*lsmean=temp;	/* Mandel temp*pred Model */
/* Mandel Margin Plots -- against temp */
proc sort data=mand1; by temp;
proc sort data=mand2; by temp;
proc sort data=lstemp; by temp;
data mandplot; merge mand1 mand2 lstemp; by temp;
proc plot;
   plot pred1*lsmean=humid;	/* Mandel temp*pred Model */
   plot pred2*lsmean=humid;	/* Mandel humid*pred Model */