data iridium;
   infile 'iridium.txt';
   input iridium strata depthcat;
   
proc glm;
   class strata depthcat;
   model iridium = strata depthcat strata*depthcat;
   output out = resids r = resid p = yhat;
   title 'Analysis on iridium';
   estimate 'diff btwn strata' strata 1 -1;
   
   lsmeans depthcat / stderr pdiff adjust = tukey;
     
proc plot;
   plot resid*yhat;