/********************* scab1.sas ***************************/
title 'Scab Disease Data';
title2 'Read Data from an Excel Spreadsheet';

proc import datafile="/home/u1407221/441s24/SAS03/ScabDisease.xlsx" 
            out=scab dbms=xlsx replace;
            getnames=yes;
/* Input data file is ScabDisease.xlsx
   Ouput data table is called scab.
   dbms=xlsx The input file is an Excel spreadsheet.
             Necessary to read an Excel spreadsheet directly under unix/linux
             Works in PC environment too except for Excel 4.0 spreadsheets
             If there are multiple sheets, use sheet="sheet1" or something.
   replace	 If the data table already exists, replace it.  Use this!
   getnames=yes   Use column names as variable names. Beware of embedded,
                  leading and trailing blanks.    
   By default, blanks in a spreadsheet are missing.  */

proc print; /* By default, use the most recently created data set. */

data potato;
     set scab;
     label infection = 'Ave percent surface covered';
     
proc freq; /* Now the most recently created data set is potato. */
     tables condition;

proc means;
     class condition;
     var infection;

proc boxplot; 
     plot infection*condition; /* Deliberately unsorted */

proc glm plots=none;
     title3 'All pairwise comparisons';
     class condition;
     model infection = condition;
     /* Pairwise multiple comparisons */ 
     lsmeans condition / pdiff tdiff; /* No correction */    
     lsmeans condition / pdiff tdiff adjust = tukey; 
     lsmeans condition / pdiff tdiff adjust = bon;
     lsmeans condition / pdiff tdiff adjust = scheffe;

proc glm order=data; /* Use order in data file, not alphabetical */
     title3 'Custom Tests';
     class condition;
     model infection = condition / clparm; 
     /* clparm gives CIs for contrasts down in the estimate statements. */
     /* Test the custom contrasts. */
     contrast '1: 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 '2a: Spring300 vs. Control'  condition 1 -1  0  0  0  0  0;
     contrast '2b: Spring600 vs. Control'  condition 1  0 -1  0  0  0  0;
     contrast '2c: Spring1200 vs. Control' condition 1  0  0 -1  0  0  0;
     contrast '2d: Fall300 vs. Control'    condition 1  0  0  0 -1  0  0;
     contrast '2e: Fall600 vs. Control'    condition 1  0  0  0  0 -1  0;
     contrast '2f: Fall1200 vs. Control'   condition 1  0  0  0  0  0 -1;
     contrast '3:  Spring vs. Control'     condition 3 -1 -1 -1  0  0  0;
     contrast '4:  Fall vs. Control'       condition 3  0  0  0 -1 -1 -1;
     contrast '5:  Spring vs. Fall'        condition 0  1  1  1 -1 -1 -1;
     
     contrast '6: Spring Amount'          condition 0  1 -1  0  0  0  0,
                                          condition 0  0  1 -1  0  0  0;

     contrast '7: Fall Amount'            condition 0  0  0  0  1 -1  0,
                                          condition 0  0  0  0  0  1 -1;
     
     contrast '8a: Spr vs. Fall for 300 lbs.' condition 0  1  0  0 -1  0  0; 
     contrast '8b: Spr vs. Fall for 600 lbs.' condition 0  0  1  0  0 -1  0;
     contrast '8c: Spr vs. Fall, 1200 lbs'    condition 0  0  0  1  0  0 -1;

     /* Estimate and CI for Control vs. Fall 1200 (F = t^2) */
     estimate 'Control vs. Fall 1200'     condition 1  0  0  0  0  0 -1;
     estimate 'Control vs. Mean of Fall'  
              condition 3  0  0  0 -1 -1 -1 / divisor = 3; 
     estimate 'Spring minus Fall'
              condition 0 1 1 1 -1 -1 -1 / 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 plots=none;
     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; 
     
run;