/********************* scab1.sas ***************************/
title 'Scab Disease Data';

data potato;
     infile '/folders/myfolders/ScabDisease.data.txt';
     length condition $ 10.; /* Default length of character values is 8. */
     input condition $  infection;
     label infection = 'Ave percent surface covered';
proc freq;
     tables condition;
 
proc means;
     class condition;
     var infection;

proc glm;
     class condition;
     model infection = condition / clparm; 
     /* clparm gives CIs for contrasts down in the estimate statements. */
     /* Pairwise multiple comparisons */
     lsmeans condition / pdiff tdiff adjust = tukey; 
     lsmeans condition / pdiff tdiff adjust = bon;
     lsmeans condition / pdiff tdiff adjust = scheffe;
     /* Test some custom contrasts. Beware of alphabetical order: 
        Control Fall1200 Fall300 Fall600 Spring1200 Spring300 Spring600*/
     contrast 'Overall test for comparison' condition 1 -1  0  0  0  0  0,
                                            condition 1  0 -1  0  0  0  0,
                                            condition 1  0  0 -1  0  0  0,
                                            condition 1  0  0  0 -1  0  0,
                                            condition 1  0  0  0  0 -1  0,
                                            condition 1  0  0  0  0  0 -1;

     contrast 'Av. of Fall vs. Control'   condition 3 -1 -1 -1  0  0  0;
     contrast 'Av. of Spring vs. Control' condition 3  0  0  0 -1 -1 -1;
     contrast 'Fall vs. Spring'           condition 0  1  1  1 -1 -1 -1;
     
     contrast 'Spring Amount'             condition 0  0  0  0  1 -1  0,
                                          condition 0  0  0  0  0  1 -1;
     
     contrast 'Fall Amount'               condition 0  1 -1  0  0  0  0,
                                          condition 0  0  1 -1  0  0  0;
                                          
     contrast 'Spr vs. Fall for 300 lbs.' condition 0  0  1  0  0 -1  0; 
     contrast 'Spr vs. Fall for 600 lbs.' condition 0  0  0  1  0  0 -1;
     contrast 'Spr vs. Fall, 1200 lbs'    condition 0  1  0  0 -1  0  0;

     /* Estimate and CI for Control vs. Fall 1200 (F = t^2) */
     estimate 'Control vs. Fall 1200'     condition 1 -1  0  0  0  0  0;
     estimate 'Control vs. Mean of Fall'  
              condition 3 -1 -1 -1  0  0  0 / divisor = 3; 

/* Scheffe critical value for a test of s contrasts is critval * (p-1)/s.
   For p=7 means and a single contrast, it's critval *  (7-1)/1 */
proc iml;
     title2 'Critical value for all possible 1-df Scheffe F-tests';
     numdf = 6;   /* p-1 = Numerator degrees of freedom for initial test */
     dendf = 25;  /* n-p = Denominator degrees of freedom for initial test */
     alpha = 0.05;
     critval = finv(1-alpha,numdf,dendf);  print critval;
     ScheffeCritval = critval*numdf; print ScheffeCritval;
     
proc iml;
     title2 'Table of Scheffe critical values for COLLECTIONS of contrasts';
     numdf = 6;   /* Numerator degrees of freedom for initial test */
     dendf = 25;  /* Denominator degrees of freedom for initial test */
     alpha = 0.05;
     critval = finv(1-alpha,numdf,dendf);
     zero = {0 0}; S_table = repeat(zero,numdf,1);  /* Make empty matrix */
     /* Label the columns */
     namz = {"Number of Contrasts in followup F-test"
             "    Scheffe Critical Value"};
     mattrib S_table colname=namz;
     do i = 1 to numdf;
        s_table(|i,1|) = i;
        s_table(|i,2|) = numdf/i *  critval;
     end;
     reset noname; /* Makes output look nicer in this case */
     print "Initial test has"  numdf " and " dendf "degrees of freedom."
           "Using significance level alpha = " alpha;
     print s_table;
     
     /* Example of subsetting vs. contrasts */
     
proc glm;
     title2 'Test differences among treatments excluding control';
     title3 'Using contrasts of all means (traditional)';
     class condition;
     model infection = condition;
     contrast 'Differences excluding control' condition 0  1 -1  0  0  0  0,
                                              condition 0  0  1 -1  0  0  0,
                                              condition 0  0  0  1 -1  0  0,
                                              condition 0  0  0  0  1 -1  0,
                                              condition 0  0  0  0  0  1 -1;
proc glm;
     title2 'Test differences among treatments excluding control';
     title3 'Subset the data by excluding control';
     where condition ne 'Control'; /* Case sensitive */
     class condition;
     model infection = condition;