options nocenter ps=64 ls=80;
data a;
   infile 'bacteria.dat' firstobs=2;
   input id bact $ temp $ run $ cage day bill leg;
/* note: ignore day (consider 21 and 22 to be the same) */

proc freq; /* this table (and others?) may help you sort out design */ 
   tables temp * bact * run / norow nocol nopercent;

proc glm;			/* full model anova */
   class bact temp run cage;
   model bill leg = bact | temp | run;
   lsmeans bact temp / stderr pdiff;
proc glm;			/* simplify model */
   class bact temp run cage;
   model bill leg = bact temp;
   lsmeans bact temp / stderr pdiff;

proc glm;			/* full model ancova */
   class bact temp run cage;
   model bill = leg bact | temp | run;
   lsmeans bact temp / stderr pdiff;
proc glm;			/* simpler ancova */
   class bact temp run cage;
   model bill = leg bact;
   lsmeans bact / stderr pdiff;
   output out = diag r = rbill p = pbill;

proc plot data = diag;
   plot rbill * pbill = temp;

/* now account for cage-to-cage variability */
proc glm;
   class bact temp run cage;
   model bill = leg bact | temp | run cage(bact * temp * run);
   lsmeans bact temp / stderr pdiff;
proc mixed;
   class bact temp run cage;
   model bill = leg bact | temp | run;
   random cage;
   lsmeans bact temp;
/* and simplify */
proc glm;
   class bact temp run cage;
   model bill = leg bact cage(bact * temp * run);
   lsmeans bact / stderr pdiff;
proc mixed;
   class bact temp run cage;
   model bill = leg bact;
   random cage;
   lsmeans bact;
run;