/********************* scab1.sas (2018 verion) **********************/ title 'Scab Disease Data based on Cochran and Cox (1958)'; data potato; infile '/folders/myfolders/441s18/Lecture/ScabDisease.data.txt'; length condition $ 10.; /* Default length of character values is 8. */ input condition $ infection; label infection = 'Ave percent surface covered'; /* Make dummy variables for condition. * Names of the dummy variables will correspond to the names of the treatments. * There are never missing values for a randomly assigned explanatory variable, so this code is safe. */ if condition = 'Control' then Control=1; else Control=0; if condition = 'Fall300' then Fall300=1; else Fall300=0; if condition = 'Fall600' then Fall600=1; else Fall600=0; if condition = 'Fall1200' then Fall1200=1; else Fall1200=0; if condition = 'Spring300' then Spring300=1; else Spring300=0; if condition = 'Spring600' then Spring600=1; else Spring600=0; if condition = 'Spring1200' then Spring1200=1; else Spring1200=0; /* Check is commented out proc freq data = potato; tables (Control--Spring1200) * condition / norow nocol nopercent; run; */ proc freq data = potato; tables condition; run; /* In the following proc means, the order=data option means that the classification levels will be displayed in the order that they appear in the data file: Control Spring300 Spring600 Spring1200 Fall300 Fall600 Fall1200 The default is alphabetical, and inconvenient: Control Fall1200 Fall300 Fall600 Spring1200 Spring300 Spring600 */ proc means data=potato order=data; class condition; var infection; run; proc glm data=potato order=data; title2 'Basic proc glm'; class condition; model infection=condition; proc reg data=potato plots=none; title2 'Custom tests with proc reg'; model infection = Control Fall300 Fall600 Fall1200 Spring300 Spring600 Spring1200 / noint; /* 1. Is amount of scab disease affected by experimental treatment? */ treatment: test Control = Fall300=Fall600=Fall1200 = Spring300=Spring600=Spring1200; /* 2. Is the amount of scab disease for each treatment different from the amount in the control condition? */ Control_vs_Spring300: test Control = Spring300; Control_vs_Spring600: test Control = Spring600; Control_vs_Spring1200: test Control = Spring1200; Control_vs_Fall300: test Control = Fall300; Control_vs_Fall600: test Control = Fall600; Control_vs_Fall1200: test Control = Fall1200; /* 3. The other 6 choose 2 = 15 comparisons are easier in proc glm */ /* 4. Is the average expected amount of scab disease in the 3 Fall conditions different from the expected amount in the control condition? */ Control_vs_Fall: test 3*Control = Fall300+Fall600+Fall1200; /* 5. Is the average expected amount of scab disease in the 3 Spring conditions different from the amount in the control condition? */ Control_vs_Spring: test 3*Control = Spring300+Spring600+Spring1200; /* 6. Is the average expected amount of scab disease in the 3 Spring conditions different from the average expected amount in the three Fall conditions? */ Spring_vs_Fall: test Spring300+Spring600+Spring1200 = Fall300+Fall600+Fall1200; /* 7. Is amount of scab disease affected by the amount of sulfur when the sulfur is applied in the Spring? This one test. */ SpringSulphur: test Spring300=Spring600=Spring1200; /* 8. Is amount of scab disease affected by the amount of sulfur when the sulfur is applied in the Fall? This one test. */ FallSulphur: test Fall300=Fall600=Fall1200; /* 9. For each amount of sulphur, is the expected amount of scab disease different depending on whether the treatment is applied in Fall or the Spring? This is 3 tests. */ Fall300_vs_Spring300: test Fall300=Spring300; Fall600_vs_Spring600: test Fall600=Spring600; Fall1200_vs_Spring1200: test Fall1200=Spring1200; run; proc glm data = potato order=data; title2 'Using proc glm -- NOT basic'; class condition; model infection = condition / clparm; /* clparm gives CIs for contrasts down in the estimate statements. */ /* Pairwise multiple comparisons */ lsmeans condition / pdiff tdiff; /* Unadjusted: Questions 2 and 3 */ lsmeans condition / pdiff tdiff adjust = tukey; lsmeans condition / pdiff tdiff adjust = bon; lsmeans condition / pdiff tdiff adjust = scheffe; /* Test some custom contrasts. Order of treatments is: Control Spring300 Spring600 Spring1200 Fall300 Fall600 Fall1200 */ contrast '1. Overall test' 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 '4. Av. of Fall vs. Control' condition 3 0 0 0 -1 -1 -1; contrast '5. Av. of Spring vs. Control' condition 3 -1 -1 -1 0 0 0; contrast '6. Fall vs. Spring' condition 0 1 1 1 -1 -1 -1; contrast '7. Spring Amount' condition 0 1 -1 0 0 0 0, condition 0 0 1 -1 0 0 0; contrast '8. Fall Amount' condition 0 0 0 0 1 -1 0, condition 0 0 0 0 0 1 -1; contrast '9a. Spr vs. Fall for 300 lbs.' condition 0 1 0 0 -1 0 0; contrast '9b. Spr vs. Fall for 600 lbs.' condition 0 0 1 0 0 -1 0; contrast '9c. 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; run; /* 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; run; 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); /* Make empty matrix */ zero = {0 0}; S_table = repeat(zero,numdf,1); /* Syntax is repeat(matrix, nrowrep, ncolrep) */ /* 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 data=potato order=data; 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; run; proc glm data=potato order=data; 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;