options ps=64 ls=80;
data wasp;
   infile 'wasp.dat' firstobs=2;
   input caste $ TL WL HH HW TH TW G1L G2Wa G2L HL G1Wb G1Wa G1H;

proc corr;	/* raw correlations among responses */
   var TL--G1H;

proc glm;	/* Univariate ANOVA and Multivariate MANOVA */
   class caste;
   model TL--G1H = caste / ss1;
   manova h=caste / printe;

proc glm;	/* Stepwise DA (2 steps) using ANOVA and ANCOVA */
   class caste;
   model G1Wa = caste / ss1;
proc glm;
   class caste;
   model HW = G1Wa caste / ss1;
proc plot;	/* Plot of Most Significant Responses */
   plot HW*G1Wa=caste;

proc stepdisc;				/* Stepwise DA */
   class caste;
proc candisc distance anova out=can;	/* Canonical DA */
   class caste;
proc plot;
   plot G1Wa*can1=caste;
   plot HW*can1=caste;

/***********************************************************************/
proc glm data=caste;		/* analysis of covariance */
   class caste;
   model HW G1Wa = caste TL caste*TL;
proc plot;
   plot HW*TL=caste;
   plot G1Wa*TL=caste;

proc glm noprint data=wasp;	/* remove size as measured by TL */
   class caste;
   model WL--G1H = TL / ss1;
   output out = resids
      r = rWL rHH rHW rTH rTW rG1L rG2Wa rG2L rHL rG1Wb rG1Wa rG1H;
proc stepdisc;	/* Stepwise DA on residuals (size removed */
   class caste;
   var rWL--rG1H;
proc candisc distance anova out=rcan;
   class caste;
   var rWL--rG1H;
proc plot;
   plot rHW*rG1Wa=caste;
   plot rG1Wa*can1=caste;
   plot rHW*can1=caste;