1 OPTIONS NONOTES NOSTIMER NOSOURCE NOSYNTAXCHECK;
61
62 /********************* scab1.sas (2018 verion) **********************/
63 title 'Scab Disease Data based on Cochran and Cox (1958)';
64
65 data potato;
66 infile '/folders/myfolders/441s18/Lecture/ScabDisease.data.txt';
67 length condition $ 10.; /* Default length of character values is 8. */
68 input condition $ infection;
69 label infection = 'Ave percent surface covered';
70 /* Make dummy variables for condition.
71 * Names of the dummy variables will correspond to the names of
72 the treatments.
73 * There are never missing values for a randomly assigned
74 explanatory variable, so this code is safe.
75 */
76 if condition = 'Control' then Control=1; else Control=0;
77 if condition = 'Fall300' then Fall300=1; else Fall300=0;
78 if condition = 'Fall600' then Fall600=1; else Fall600=0;
79 if condition = 'Fall1200' then Fall1200=1; else Fall1200=0;
80 if condition = 'Spring300' then Spring300=1; else Spring300=0;
81 if condition = 'Spring600' then Spring600=1; else Spring600=0;
82 if condition = 'Spring1200' then Spring1200=1; else Spring1200=0;
83
84 /* Check is commented out
85 proc freq data = potato;
86 tables (Control--Spring1200) * condition / norow nocol nopercent;
87 run;
88 */
89
NOTE: The infile '/folders/myfolders/441s18/Lecture/ScabDisease.data.txt' is:
Filename=/folders/myfolders/441s18/Lecture/ScabDisease.data.txt,
Owner Name=root,Group Name=vboxsf,
Access Permission=-rwxrwx---,
Last Modified=11Feb2018:15:28:04,
File Size (bytes)=578
NOTE: 33 records were read from the infile '/folders/myfolders/441s18/Lecture/ScabDisease.data.txt'.
The minimum record length was 0.
The maximum record length was 16.
NOTE: SAS went to a new line when INPUT statement reached past the end of a line.
NOTE: The data set WORK.POTATO has 32 observations and 9 variables.
NOTE: DATA statement used (Total process time):
real time 0.00 seconds
cpu time 0.00 seconds
90 proc freq data = potato;
91 tables condition;
92 run;
NOTE: There were 32 observations read from the data set WORK.POTATO.
NOTE: PROCEDURE FREQ used (Total process time):
real time 0.04 seconds
cpu time 0.04 seconds
93
94 /* In the following proc means, the order=data option means that the
95 classification levels will be displayed in the order that they appear
96 in the data file:
97 Control Spring300 Spring600 Spring1200 Fall300 Fall600 Fall1200
98 The default is alphabetical, and inconvenient:
99 Control Fall1200 Fall300 Fall600 Spring1200 Spring300 Spring600
100 */
101
102 proc means data=potato order=data;
103 class condition;
104 var infection;
105 run;
NOTE: There were 32 observations read from the data set WORK.POTATO.
NOTE: PROCEDURE MEANS used (Total process time):
real time 0.03 seconds
cpu time 0.02 seconds
106
107 proc glm data=potato order=data;
108 title2 'Basic proc glm';
109 class condition;
110 model infection=condition;
111
NOTE: PROCEDURE GLM used (Total process time):
real time 0.36 seconds
cpu time 0.19 seconds
112 proc reg data=potato plots=none;
113 title2 'Custom tests with proc reg';
114 model infection = Control Fall300 Fall600 Fall1200
115 Spring300 Spring600 Spring1200 / noint;
116 /* 1. Is amount of scab disease affected by experimental treatment? */
117 treatment: test Control = Fall300=Fall600=Fall1200
118 = Spring300=Spring600=Spring1200;
119
120 /* 2. Is the amount of scab disease for each treatment different from the
121 amount in the control condition? */
122 Control_vs_Spring300: test Control = Spring300;
123 Control_vs_Spring600: test Control = Spring600;
124 Control_vs_Spring1200: test Control = Spring1200;
125 Control_vs_Fall300: test Control = Fall300;
126 Control_vs_Fall600: test Control = Fall600;
127 Control_vs_Fall1200: test Control = Fall1200;
128
129 /* 3. The other 6 choose 2 = 15 comparisons are easier in proc glm */
130
131 /* 4. Is the average expected amount of scab disease in the 3 Fall
132 conditions different from the expected amount in the control
133 condition? */
134 Control_vs_Fall: test 3*Control = Fall300+Fall600+Fall1200;
135
136 /* 5. Is the average expected amount of scab disease in the 3 Spring
137 conditions different from the amount in the control condition? */
138 Control_vs_Spring: test 3*Control = Spring300+Spring600+Spring1200;
139
140 /* 6. Is the average expected amount of scab disease in the 3 Spring
141 conditions different from the average expected amount in the
142 three Fall conditions? */
143 Spring_vs_Fall: test Spring300+Spring600+Spring1200 =
144 Fall300+Fall600+Fall1200;
145
146 /* 7. Is amount of scab disease affected by the amount of sulfur when the
147 sulfur is applied in the Spring? This one test. */
148 SpringSulphur: test Spring300=Spring600=Spring1200;
149
150 /* 8. Is amount of scab disease affected by the amount of sulfur when the
151 sulfur is applied in the Fall? This one test. */
152 FallSulphur: test Fall300=Fall600=Fall1200;
153
154 /* 9. For each amount of sulphur, is the expected amount of scab disease
155 different depending on whether the treatment is applied in Fall or
156 the Spring? This is 3 tests. */
157 Fall300_vs_Spring300: test Fall300=Spring300;
158 Fall600_vs_Spring600: test Fall600=Spring600;
159 Fall1200_vs_Spring1200: test Fall1200=Spring1200;
160 run;
161
NOTE: PROCEDURE REG used (Total process time):
real time 0.26 seconds
cpu time 0.24 seconds
162 proc glm data = potato order=data;
163 title2 'Using proc glm -- NOT basic';
164 class condition;
165 model infection = condition / clparm;
166 /* clparm gives CIs for contrasts down in the estimate statements. */
167 /* Pairwise multiple comparisons */
168 lsmeans condition / pdiff tdiff; /* Unadjusted: Questions 2 and 3 */
169 lsmeans condition / pdiff tdiff adjust = tukey;
170 lsmeans condition / pdiff tdiff adjust = bon;
171 lsmeans condition / pdiff tdiff adjust = scheffe;
172 /* Test some custom contrasts. Order of treatments is:
173 Control Spring300 Spring600 Spring1200 Fall300 Fall600 Fall1200 */
174 contrast '1. Overall test' condition 1 -1 0 0 0 0 0,
175 condition 1 0 -1 0 0 0 0,
176 condition 1 0 0 -1 0 0 0,
177 condition 1 0 0 0 -1 0 0,
178 condition 1 0 0 0 0 -1 0,
179 condition 1 0 0 0 0 0 -1;
180 contrast '4. Av. of Fall vs. Control' condition 3 0 0 0 -1 -1 -1;
181 contrast '5. Av. of Spring vs. Control' condition 3 -1 -1 -1 0 0 0;
182 contrast '6. Fall vs. Spring' condition 0 1 1 1 -1 -1 -1;
183
184 contrast '7. Spring Amount' condition 0 1 -1 0 0 0 0,
185 condition 0 0 1 -1 0 0 0;
186
187 contrast '8. Fall Amount' condition 0 0 0 0 1 -1 0,
188 condition 0 0 0 0 0 1 -1;
189
190 contrast '9a. Spr vs. Fall for 300 lbs.' condition 0 1 0 0 -1 0 0;
191 contrast '9b. Spr vs. Fall for 600 lbs.' condition 0 0 1 0 0 -1 0;
192 contrast '9c. Spr vs. Fall, 1200 lbs' condition 0 0 0 1 0 0 -1;
193
194 /* Estimate and CI for Control vs. Fall 1200 (F = t^2) */
195 estimate 'Control vs. Fall 1200'
196 condition 1 0 0 0 0 0 -1;
197 estimate 'Control vs. Mean of Fall'
198 condition 3 0 0 0 -1 -1 -1 / divisor = 3;
199 run;
200
201 /* Scheffe critical value for a test of s contrasts is critval * (p-1)/s.
202 For p=7 means and a single contrast, it's critval * (7-1)/1 */
NOTE: PROCEDURE GLM used (Total process time):
real time 1.92 seconds
cpu time 1.09 seconds
203 proc iml;
NOTE: IML Ready
204 title2 'Critical value for all possible 1-df Scheffe F-tests';
205 numdf = 6;
205 ! /* p-1 = Numerator degrees of freedom for initial test */
206 dendf = 25;
206 ! /* n-p = Denominator degrees of freedom for initial test */
207 alpha = 0.05;
208 critval = finv(1-alpha,numdf,dendf);
208 ! print critval;
209 ScheffeCritval = critval*numdf;
209 ! print ScheffeCritval;
210 run;
NOTE: Module MAIN is undefined in IML; cannot be RUN.
211
NOTE: Exiting IML.
NOTE: PROCEDURE IML used (Total process time):
real time 0.02 seconds
cpu time 0.02 seconds
212 proc iml;
NOTE: IML Ready
213 title2 'Table of Scheffe critical values for COLLECTIONS of contrasts';
214 numdf = 6;
214 ! /* Numerator degrees of freedom for initial test */
215 dendf = 25;
215 ! /* Denominator degrees of freedom for initial test */
216 alpha = 0.05;
217 critval = finv(1-alpha,numdf,dendf);
218 /* Make empty matrix */
219 zero = {0 0};
220 S_table = repeat(zero,numdf,1);
221 /* Syntax is repeat(matrix, nrowrep, ncolrep) */
222 /* Label the columns */
223 namz = {"Number of Contrasts in followup F-test"
224 " Scheffe Critical Value"};
225 mattrib S_table colname=namz;
226 do i = 1 to numdf;
227 s_table(|i,1|) = i;
228 s_table(|i,2|) = numdf/i * critval;
229 end;
230 reset noname;
230 ! /* Makes output look nicer in this case */
231 print "Initial test has" numdf " and " dendf "degrees of freedom."
232 "Using significance level alpha = " alpha;
233 print s_table;
234
235 /******** Example of subsetting vs. contrasts ********/
236
NOTE: Exiting IML.
NOTE: PROCEDURE IML used (Total process time):
real time 0.02 seconds
cpu time 0.03 seconds
237 proc glm data=potato order=data;
238 title2 'Test differences among treatments excluding control';
239 title3 'Using contrasts of all means (traditional)';
240 class condition;
241 model infection = condition;
242 contrast 'Differences excluding control' condition 0 1 -1 0 0 0 0,
243 condition 0 0 1 -1 0 0 0,
244 condition 0 0 0 1 -1 0 0,
245 condition 0 0 0 0 1 -1 0,
246 condition 0 0 0 0 0 1 -1;
247 run;
NOTE: PROCEDURE GLM used (Total process time):
real time 0.32 seconds
cpu time 0.18 seconds
248 proc glm data=potato order=data;
249 title2 'Test differences among treatments excluding control';
250 title3 'Subset the data by excluding control';
251 where condition ne 'Control'; /* Case sensitive */
252 class condition;
253 model infection = condition;
254 run;
255
256
257
258
259
260
261 OPTIONS NONOTES NOSTIMER NOSOURCE NOSYNTAXCHECK;
274