options formdlim = '-';

data pygm;
  infile '/home/pdixon/WWW/stat401b/data/pygmalion.txt';
  input company treat score;

proc glm;           
  class company treat;
  model score =  company treat treat*company;
  title 'Fit model with interaction';
  lsmeans treat /stderr;
  means treat;
  
  estimate 'diff btwn treatments' treat 1 -1;
  
proc glm;           
  class company treat;
  model score =  company treat;
  title 'Fit model without interaction';
  lsmeans treat /stderr;
  means treat;
  
  estimate 'diff btwn treatments' treat 1 -1;

/* plot cell means for combination of company and treat */
/*   use company as the X axis.  Assumes that company codes are numbers */
/*   plots WILL NOT WORK if the variable on the X axis is character */

proc sort;
  by company treat;
     
proc means noprint;
  by company treat;
  var score;                  /* this is the name of the response variable */
  output out = means mean = mean;
     
proc plot;
  plot mean*company = treat;  /* plot means for each company and treat */