Multiple Regression – Muscle Study - SAS Program
options ps=54 ls=76 nodate nonumber;
data muscle2;
* infile 'C:\data\muscle2.dat';
infile '../data/muscle2.dat';
input m w h;
run;
proc plot hpct=50 vpct=50;
plot m*w m*h w*h;
run;
proc reg;
model h = m w / r influence vif;
run;
quit;
Multiple Regression – Muscle Study - SAS Output
The SAS System
Plot of m*w. A=1, B=2, etc. Plot of m*h. A=1, B=2, etc.
m | m |
80 + A 80 + A
| |
| A | A
| A B | A A A
| A A A A | BAA
70 + A C AC B 70 + A AA B AAA A A
| A C A | AA B A
| A A AA | A A A A
| AA | A A
| A B A | A A AA
60 + A 60 + A
| A A | A A
| |
| |
| |
50 + 50 +
--+------+------+------+------+-- -+------+------+-
75 100 125 150 175 2500 3000 3500
w h
Plot of w*h. A=1, B=2, etc.
175 +
|
| A A
| A A
150 +
| A A B AB A A A
w | AA A
| A AAAB B AB
125 +
| A A A
| A A A
|
100 + A A A
|
|
|
75 +
-+------+------+
2500 3000 3500
h
The REG Procedure
Model: MODEL1
Dependent Variable: h
Number of Observations Read 37
Number of Observations Used 37
Analysis of Variance
Sum of Mean
Source DF Squares Square F Value Pr > F
Model 2 713520 356760 16.48 <.0001
Error 34 736024 21648
Corrected Total 36 1449544
Root MSE 147.13176 R-Square 0.4922
Dependent Mean 3015.70270 Adj R-Sq 0.4624
Coeff Var 4.87885
Parameter Estimates
Parameter Standard Variance
Variable DF Estimate Error t Value Pr > |t| Inflation
Intercept 1 977.42536 376.05311 2.60 0.0137 0
m 1 17.77770 4.94281 3.60 0.0010 1.00965
w 1 6.24364 1.52213 4.10 0.0002 1.00965
The REG Procedure
Model: MODEL1
Dependent Variable: h
Output Statistics
Dependent Predicted Std Error Std Error Student
Obs Variable Value Mean Predict Residual Residual Residual
1 3398 3311 57.0719 86.9116 135.6 0.641
2 2988 2957 41.0637 30.6256 141.3 0.217
3 3048 3105 29.2528 -57.0960 144.2 -0.396
4 2781 2899 59.1551 -117.8746 134.7 -0.875
5 2912 3194 41.7618 -281.9845 141.1 -1.999
6 3135 3003 25.3345 131.6324 144.9 0.908
7 3261 3096 28.8805 164.7929 144.3 1.142
8 3030 3060 55.4891 -30.4228 136.3 -0.223
9 3139 3000 25.1862 139.1880 145.0 0.960
10 2996 3030 32.4587 -34.4296 143.5 -0.240
11 3248 3030 27.5737 217.9659 144.5 1.508
12 3117 3096 28.8805 20.7929 144.3 0.144
13 2891 3000 37.1307 -109.2075 142.4 -0.767
14 2667 2794 45.7559 -126.8196 139.8 -0.907
15 3403 3174 37.9975 228.5710 142.1 1.608
16 2999 2875 41.5245 123.8750 141.2 0.878
17 3318 3082 26.8921 235.8577 144.7 1.631
18 2989 3153 66.1218 -164.4391 131.4 -1.251
19 2936 2875 41.5245 60.8750 141.2 0.431
20 3020 2988 30.7577 32.0796 143.9 0.223
21 2812 2857 46.0121 -45.4427 139.8 -0.325
22 2962 2871 64.2723 91.2119 132.4 0.689
23 3236 3010 24.9516 225.7644 145.0 1.557
24 3214 3032 26.5824 182.4312 144.7 1.261
25 3389 3246 51.3027 143.3173 137.9 1.039
26 2908 2968 25.7989 -59.5691 144.9 -0.411
27 3063 3022 26.6143 41.3019 144.7 0.285
28 2956 3051 29.0770 -95.1243 144.2 -0.660
29 3023 3174 36.6978 -150.6044 142.5 -1.057
30 3001 3017 25.3530 -16.3467 144.9 -0.113
31 2841 3085 31.9054 -243.7159 143.6 -1.697
32 3117 3224 50.6348 -107.3494 138.1 -0.777
33 2733 2912 30.5245 -178.6623 143.9 -1.241
34 2808 3059 44.2339 -251.2173 140.3 -1.790
35 2813 2787 58.7554 25.7671 134.9 0.191
36 2615 2606 75.4471 9.0997 126.3 0.0720
37 2814 2936 40.3077 -121.7552 141.5 -0.860
The REG Procedure
Model: MODEL1
Dependent Variable: h
Output Statistics
Cook's Hat Diag Cov
Obs -2-1 0 1 2 D RStudent H Ratio DFFITS
1 | |* | 0.024 0.6352 0.1505 1.2413 0.2673
2 | | | 0.001 0.2137 0.0779 1.1812 0.0621
3 | | | 0.002 -0.3910 0.0395 1.1230 -0.0793
4 | *| | 0.049 -0.8719 0.1616 1.2184 -0.3829
5 | ***| | 0.117 -2.0961 0.0806 0.8176 -0.6205
6 | |* | 0.008 0.9058 0.0296 1.0470 0.1583
7 | |** | 0.017 1.1476 0.0385 1.0115 0.2297
8 | | | 0.003 -0.2201 0.1422 1.2694 -0.0896
9 | |* | 0.009 0.9590 0.0293 1.0375 0.1666
10 | | | 0.001 -0.2366 0.0487 1.1438 -0.0535
11 | |*** | 0.028 1.5382 0.0351 0.9209 0.2935
12 | | | 0.000 0.1420 0.0385 1.1354 0.0284
13 | *| | 0.013 -0.7623 0.0637 1.1085 -0.1988
14 | *| | 0.029 -0.9045 0.0967 1.1250 -0.2960
15 | |*** | 0.062 1.6482 0.0667 0.9243 0.4406
16 | |* | 0.022 0.8746 0.0797 1.1094 0.2573
17 | |*** | 0.031 1.6731 0.0334 0.8863 0.3110
18 | **| | 0.132 -1.2619 0.2020 1.1898 -0.6348
19 | | | 0.005 0.4261 0.0797 1.1689 0.1253
20 | | | 0.001 0.2198 0.0437 1.1387 0.0470
21 | | | 0.004 -0.3208 0.0978 1.2010 -0.1056
22 | |* | 0.037 0.6837 0.1908 1.2958 0.3320
23 | |*** | 0.024 1.5917 0.0288 0.9020 0.2739
24 | |** | 0.018 1.2721 0.0326 0.9793 0.2337
25 | |** | 0.050 1.0406 0.1216 1.1301 0.3871
26 | | | 0.002 -0.4062 0.0307 1.1116 -0.0723
27 | | | 0.001 0.2815 0.0327 1.1226 0.0518
28 | *| | 0.006 -0.6540 0.0391 1.0950 -0.1318
29 | **| | 0.025 -1.0589 0.0622 1.0550 -0.2727
30 | | | 0.000 -0.1111 0.0297 1.1259 -0.0194
31 | ***| | 0.047 -1.7473 0.0470 0.8801 -0.3881
32 | *| | 0.027 -0.7725 0.1184 1.1757 -0.2831
33 | **| | 0.023 -1.2516 0.0430 0.9944 -0.2654
34 | ***| | 0.106 -1.8532 0.0904 0.8934 -0.5842
35 | | | 0.002 0.1883 0.1595 1.2970 0.0820
36 | | | 0.001 0.0710 0.2629 1.4832 0.0424
37 | *| | 0.020 -0.8571 0.0751 1.1069 -0.2441
The SAS System
The REG Procedure
Model: MODEL1
Dependent Variable: h
Output Statistics
------DFBETAS------
Obs Intercept m w
1 -0.2278 0.1692 0.1558
2 0.0014 0.0276 -0.0444
3 0.0310 -0.0152 -0.0402
4 -0.2390 0.3347 -0.1323
5 0.4795 -0.4506 -0.1846
6 0.0078 0.0268 -0.0412
7 -0.0743 0.0250 0.1200
8 -0.0164 0.0542 -0.0647
9 0.0138 0.0218 -0.0429
10 -0.0100 0.0250 -0.0278
11 -0.0533 0.1245 -0.0778
12 -0.0092 0.0031 0.0148
13 -0.0705 0.1261 -0.0947
14 -0.2541 0.1760 0.1612
15 -0.3210 0.2886 0.1502
16 0.1994 -0.2072 -0.0077
17 -0.1160 0.1023 0.0791
18 0.3784 -0.5684 0.2161
19 0.0971 -0.1010 -0.0037
20 0.0184 -0.0265 0.0144
21 -0.0839 0.0894 0.0006
22 0.0411 0.1355 -0.2881
23 0.0055 0.0452 -0.0539
24 -0.0396 0.0883 -0.0483
25 -0.2223 0.0657 0.3270
26 -0.0295 0.0217 0.0105
27 0.0081 -0.0155 0.0165
28 0.0412 -0.0702 0.0273
29 0.1869 -0.1413 -0.1342
30 0.0009 -0.0047 0.0039
31 0.0523 0.0694 -0.2490
32 0.1368 -0.0155 -0.2456
33 -0.1503 0.0698 0.1386
34 0.2288 -0.4359 0.2634
35 0.0384 0.0042 -0.0747
36 0.0394 -0.0266 -0.0273
37 -0.0321 -0.0845 0.1835
Sum of Residuals 0
Sum of Squared Residuals 736024
Predicted Residual SS (PRESS) 852966
Multiple Regression – Muscle Study - R Program
pdf("muscle2.pdf")
muscle2 <- read.table("C:\\data\\muscle2.dat", header=F,
col.names=c("M","W","H"))
pairs(muscle2)
muscle2.reg <- lm(H ~ M + W, data=muscle2)
summary(muscle2.reg)
aov(muscle2.reg)
par(mfrow=c(2,2))
plot(muscle2.reg, which=1:4)
muscle2.rstandard <- rstandard(muscle2.reg)
muscle2.rstudent <- rstudent(muscle2.reg)
muscle2.inf <- influence.measures(muscle2.reg)
muscle2.rstandard
muscle2.rstudent
muscle2.inf
## You must have downloaded DAAG to get Variance Inflation Factors
library(DAAG)
muscle2.vif <- vif(muscle2.reg)
muscle2.vif
dev.off()
------
Multiple Regression – Muscle Study - R Output
> summary(muscle2.reg)
Call:
lm(formula = H ~ M + W, data = muscle2)
Residuals:
Min 1Q Median 3Q Max
-282.0 -109.2 9.1 123.9 235.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 977.425 376.053 2.599 0.013723 *
M 17.778 4.943 3.597 0.001011 **
W 6.244 1.522 4.102 0.000242 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 147.1 on 34 degrees of freedom
Multiple R-squared: 0.4922, Adjusted R-squared: 0.4624
F-statistic: 16.48 on 2 and 34 DF, p-value: 9.914e-06
> aov(muscle2.reg)
Call:
aov(formula = muscle2.reg)
Terms:
M W Residuals
Sum of Squares 349284.5 364235.6 736023.6
Deg. of Freedom 1 1 34
Residual standard error: 147.1318
Estimated effects may be unbalanced
> par(mfrow=c(2,2))
> plot(muscle2.reg, which=1:4)
> muscle2.rstandard <- rstandard(muscle2.reg)
> muscle2.rstudent <- rstudent(muscle2.reg)
> muscle2.inf <- influence.measures(muscle2.reg)
> muscle2.rstandard
1 2 3 4 5 6
0.64088516 0.21676409 -0.39596519 -0.87498523 -1.99874900 0.90822206
7 8 9 10 11 12
1.14225776 -0.22325885 0.96018174 -0.23991633 1.50815450 0.14412538
13 14 15 16 17 18
-0.76707120 -0.90691533 1.60806293 0.87760899 1.63050403 -1.25108759
19 20 21 22 23 24
0.43127694 0.22295904 -0.32516622 0.68916608 1.55698960 1.26066295
25 26 27 28 29 30
1.03930157 -0.41124039 0.28542219 -0.65953210 -1.05700926 -0.11278928
31 32 33 34 35 36
-1.69682257 -0.77708127 -1.24130914 -1.79025245 0.19102184 0.07203953
37
-0.86044390
> muscle2.rstudent
1 2 3 4 5 6
0.63523866 0.21370031 -0.39100130 -0.87189416 -2.09613088 0.90582162
7 8 9 10 11 12
1.14756891 -0.22011253 0.95904825 -0.23656214 1.53815029 0.14203346
13 14 15 16 17 18
-0.76233173 -0.90448601 1.64815278 0.87456915 1.67308982 -1.26194146
19 20 21 22 23 24
0.42605428 0.21981651 -0.32084795 0.68374812 1.59171611 1.27207161
25 26 27 28 29 30
1.04056596 -0.40615900 0.28153096 -0.65395742 -1.05889237 -0.11113903
31 32 33 34 35 36
-1.74730240 -0.77245866 -1.25160803 -1.85323788 0.18829280 0.07097764
37
-0.85707886
> muscle2.inf
Influence measures of
lm(formula = H ~ M + W, data = muscle2) :
dfb.1_ dfb.M dfb.W dffit cov.r cook.d hat inf
1 -0.227811 0.16922 0.155829 0.2673 1.241 0.024249 0.1505
2 0.001412 0.02762 -0.044410 0.0621 1.181 0.001323 0.0779
3 0.030986 -0.01523 -0.040243 -0.0793 1.123 0.002151 0.0395
4 -0.239043 0.33475 -0.132321 -0.3829 1.218 0.049207 0.1616
5 0.479499 -0.45062 -0.184612 -0.6205 0.818 0.116686 0.0806
6 0.007809 0.02676 -0.041171 0.1583 1.047 0.008401 0.0296
7 -0.074315 0.02501 0.119971 0.2297 1.012 0.017429 0.0385
8 -0.016397 0.05422 -0.064739 -0.0896 1.269 0.002755 0.1422 *
9 0.013783 0.02179 -0.042943 0.1666 1.038 0.009277 0.0293
10 -0.010028 0.02497 -0.027809 -0.0535 1.144 0.000982 0.0487
11 -0.053259 0.12449 -0.077822 0.2935 0.921 0.027598 0.0351
12 -0.009198 0.00309 0.014849 0.0284 1.135 0.000277 0.0385
13 -0.070523 0.12609 -0.094729 -0.1988 1.108 0.013341 0.0637
14 -0.254132 0.17604 0.161166 -0.2960 1.125 0.029354 0.0967
15 -0.321034 0.28863 0.150237 0.4406 0.924 0.061597 0.0667
16 0.199413 -0.20723 -0.007693 0.2573 1.109 0.022219 0.0797
17 -0.116027 0.10230 0.079072 0.3110 0.886 0.030628 0.0334
18 0.378448 -0.56839 0.216104 -0.6348 1.190 0.132041 0.2020
19 0.097146 -0.10096 -0.003747 0.1253 1.169 0.005366 0.0797
20 0.018432 -0.02651 0.014356 0.0470 1.139 0.000757 0.0437
21 -0.083948 0.08937 0.000582 -0.1056 1.201 0.003820 0.0978
22 0.041063 0.13548 -0.288112 0.3320 1.296 0.037335 0.1908 *
23 0.005476 0.04521 -0.053933 0.2739 0.902 0.023928 0.0288
24 -0.039593 0.08835 -0.048285 0.2337 0.979 0.017876 0.0326
25 -0.222333 0.06569 0.326994 0.3871 1.130 0.049834 0.1216
26 -0.029455 0.02171 0.010525 -0.0723 1.112 0.001788 0.0307
27 0.008061 -0.01549 0.016496 0.0518 1.123 0.000919 0.0327
28 0.041216 -0.07022 0.027315 -0.1318 1.095 0.005893 0.0391
29 0.186894 -0.14127 -0.134178 -0.2727 1.055 0.024706 0.0622
30 0.000896 -0.00467 0.003916 -0.0194 1.126 0.000130 0.0297
31 0.052339 0.06941 -0.249024 -0.3881 0.880 0.047357 0.0470
32 0.136778 -0.01551 -0.245551 -0.2831 1.176 0.027042 0.1184
33 -0.150290 0.06976 0.138596 -0.2654 0.994 0.023101 0.0430
34 0.228751 -0.43587 0.263445 -0.5842 0.893 0.106157 0.0904
35 0.038447 0.00418 -0.074678 0.0820 1.297 0.002308 0.1595 *
36 0.039372 -0.02665 -0.027292 0.0424 1.483 0.000617 0.2629 *
37 -0.032102 -0.08445 0.183504 -0.2441 1.107 0.020025 0.0751
> ## You must have downloaded DAAG to get Variance Inflation Factors
> library(DAAG)
> muscle2.vif <- vif(muscle2.reg)
> muscle2.vif
M W
1.0096 1.0096
1-Way ANOVA – CRD – Amoeba Study – SAS Program
options ps=54 ls=76;
data amoeba;
*infile 'C:\data\entozamoeba.dat';
infile '../../data/entozamoeba.dat';
input trt yield;
run;
proc glm;
class trt;
model yield = trt;
means trt / bon tukey hovtest;
run;
proc npar1way wilcoxon;
class trt;
var yield;
run;
quit;
1-Way ANOVA – CRD – Amoeba Study – SAS Output
The GLM Procedure
Class Level Information
Class Levels Values
trt 5 1 2 3 4 5
Number of Observations Read 50
Number of Observations Used 50
The GLM Procedure
Dependent Variable: yield
Sum of
Source DF Squares Mean Square F Value Pr > F
Model 4 19665.68000 4916.42000 5.04 0.0019
Error 45 43858.40000 974.63111
Corrected Total 49 63524.08000
R-Square Coeff Var Root MSE yield Mean
0.309578 14.81543 31.21908 210.7200
Source DF Type I SS Mean Square F Value Pr > F
trt 4 19665.68000 4916.42000 5.04 0.0019
Source DF Type III SS Mean Square F Value Pr > F
trt 4 19665.68000 4916.42000 5.04 0.0019
The GLM Procedure
Levene's Test for Homogeneity of yield Variance
ANOVA of Squared Deviations from Group Means
Sum of Mean
Source DF Squares Square F Value Pr > F
trt 4 10086852 2521713 1.73 0.1604
Error 45 65661767 1459150
Tukey's Studentized Range (HSD) Test for yield
NOTE: This test controls the Type I experimentwise error rate, but it
generally has a higher Type II error rate than REGWQ.
Alpha 0.05
Error Degrees of Freedom 45
Error Mean Square 974.6311
Critical Value of Studentized Range 4.01842
Minimum Significant Difference 39.671
Means with the same letter are not significantly different.
Tukey Grouping Mean N trt
A 246.40 10 1
A
B A 216.50 10 4
B
B 203.50 10 2
B
B 196.90 10 3
B
B 190.30 10 5
Bonferroni (Dunn) t Tests for yield
NOTE: This test controls the Type I experimentwise error rate, but it
generally has a higher Type II error rate than REGWQ.
Alpha 0.05
Error Degrees of Freedom 45
Error Mean Square 974.6311
Critical Value of t 2.95208
Minimum Significant Difference 41.216
Means with the same letter are not significantly different.
Bon Grouping Mean N trt
A 246.40 10 1
A
B A 216.50 10 4
B
B 203.50 10 2
B
B 196.90 10 3
B
B 190.30 10 5
The NPAR1WAY Procedure
Wilcoxon Scores (Rank Sums) for Variable yield
Classified by Variable trt
Sum of Expected Std Dev Mean
trt N Scores Under H0 Under H0 Score
------
1 10 386.50 255.0 41.221156 38.650
2 10 230.50 255.0 41.221156 23.050
3 10 211.50 255.0 41.221156 21.150
4 10 274.00 255.0 41.221156 27.400
5 10 172.50 255.0 41.221156 17.250
Average scores were used for ties.
Kruskal-Wallis Test
Chi-Square 12.6894
DF 4
Pr > Chi-Square 0.0129
------
1-Way ANOVA – CRD – Amoeba Study – R Program
pdf("C:\\Rmisc\\graphs\\amoeba.pdf")
amoeba1 <- read.fwf("C:\\data\\amoeba.txt", width=c(8,8), col.names=c("Trt", "Yield"))
fTrt <- factor(amoeba1$Trt, levels=1:5)
levels(fTrt) <- c("None", "H", "F", "HF", "FH")
amoeba <- data.frame(amoeba1, fTrt)
attach(amoeba)
kruskal.test(Yield~fTrt)
amoeba1.aov <- aov(Yield~fTrt)
summary(amoeba1.aov)
# Obtain Tukey's Comparisons among levels of treatment
TukeyHSD(amoeba1.aov, "fTrt")
# Obtain Bonferroni's Comparisons among levels of treatment (does not use MSE from ANOVA and not pooled df)
pairwise.t.test(Yield, fTrt, p.adj="bonf")
# Histogram of Residuals and Plot Residuals versus Time order
hist(residuals(amoeba1.aov), breaks=seq(-75,75,15))
plot(residuals(amoeba1.aov), type="o")
# Bartlett's test for equal variances
bartlett.test(Yield~fTrt)
dev.off()
1-Way ANOVA – CRD – Amoeba Study – R Output
> kruskal.test(Yield~fTrt)
Kruskal-Wallis rank sum test
data: Yield by fTrt
Kruskal-Wallis chi-squared = 12.6894, df = 4, p-value = 0.01290
amoeba1.aov <- aov(Yield~fTrt)
> summary(amoeba1.aov)
Df Sum Sq Mean Sq F value Pr(>F)
fTrt 4 19666 4916.4 5.0444 0.001911 **
Residuals 45 43858 974.6
---
> # Obtain Tukey's Comparisons among levels of treatment
> TukeyHSD(amoeba1.aov, "fTrt")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Yield ~ fTrt)
$fTrt
diff lwr upr p adj
H-None -42.9 -82.57118 -3.228817 0.0281550
F-None -49.5 -89.17118 -9.828817 0.0078618
HF-None -29.9 -69.57118 9.771183 0.2208842
FH-None -56.1 -95.77118 -16.428817 0.0019658
F-H -6.6 -46.27118 33.071183 0.9894377
HF-H 13.0 -26.67118 52.671183 0.8832844
FH-H -13.2 -52.87118 26.471183 0.8774716
HF-F 19.6 -20.07118 59.271183 0.6284290
FH-F -6.6 -46.27118 33.071183 0.9894377
FH-HF -26.2 -65.87118 13.471183 0.3444453
> # Obtain Bonferroni's Comparisons among levels of treatment (does not use MSE from ANOVA and not pooled df)
> pairwise.t.test(Yield, fTrt, p.adj="bonf")
Pairwise comparisons using t tests with pooled SD
data: Yield and fTrt
None H F HF
H 0.0360 - - -
F 0.0093 1.0000 - -
HF 0.3768 1.0000 1.0000 -
FH 0.0022 1.0000 1.0000 0.6707
P value adjustment method: bonferroni
> # Bartlett's test for equal variances
> bartlett.test(Yield~fTrt)
Bartlett test of homogeneity of variances
data: Yield by fTrt
Bartlett's K-squared = 8.2024, df = 4, p-value = 0.08444
2-Way ANOVA – RBD – Caffeine Study – SAS Program
options ps=54 ls=76;
goptions reset=all device=pdf gsfname=output gsfmode=replace;
data caffeine;
infile '../../data/caffeine.dat';
* infile 'C:\caffeine.dat';
input subject dose endurance;
run;
filename output 'caffeine1.pdf';
proc gplot;
plot endurance*subject=dose;
run;
proc glm;
class subject dose;
model endurance = dose subject;
means dose / bon tukey;
output out=caffout r=caffres p=caffpred;
run;
filename output 'caffeine2.pdf';
proc gchart;
vbar caffres;
run;
/* The statistic for Friedman's Test is Statistic 2:
"Row Mean Scores Differ". In the tables statement, use the ordering:
block*treatment*response */
proc freq;
tables subject*dose*endurance / cmh2 scores=rank noprint;
run;
quit;
2-Way ANOVA – RBD – Caffeine Study – SAS Output
The GLM Procedure
Class Level Information
Class Levels Values
subject 9 1 2 3 4 5 6 7 8 9
dose 4 0 5 9 13
Number of Observations Read 36
Number of Observations Used 36
The GLM Procedure
Dependent Variable: endurance
Sum of
Source DF Squares Mean Square F Value Pr > F
Model 11 6491.115722 590.101429 11.23 <.0001
Error 24 1261.657078 52.569045
Corrected Total 35 7752.772800
R-Square Coeff Var Root MSE endurance Mean
0.837264 13.12616 7.250451 55.23667
Source DF Type I SS Mean Square F Value Pr > F
dose 3 933.121622 311.040541 5.92 0.0036
subject 8 5557.994100 694.749263 13.22 <.0001
Source DF Type III SS Mean Square F Value Pr > F
dose 3 933.121622 311.040541 5.92 0.0036
subject 8 5557.994100 694.749263 13.22 <.0001
Tukey's Studentized Range (HSD) Test for endurance
NOTE: This test controls the Type I experimentwise error rate, but it
generally has a higher Type II error rate than REGWQ.
Alpha 0.05
Error Degrees of Freedom 24
Error Mean Square 52.56904
Critical Value of Studentized Range 3.90126
Minimum Significant Difference 9.4286
Means with the same letter are not significantly different.
Tukey Grouping Mean N dose
A 58.681 9 9
A
A 58.149 9 13
A
A 57.677 9 5
B 46.440 9 0
The GLM Procedure
Bonferroni (Dunn) t Tests for endurance
NOTE: This test controls the Type I experimentwise error rate, but it
generally has a higher Type II error rate than REGWQ.
Alpha 0.05
Error Degrees of Freedom 24
Error Mean Square 52.56904
Critical Value of t 2.87509
Minimum Significant Difference 9.8268
Means with the same letter are not significantly different.
Bon Grouping Mean N dose
A 58.681 9 9
A
A 58.149 9 13
A
A 57.677 9 5
B 46.440 9 0
The FREQ Procedure
Summary Statistics for dose by endurance
Controlling for subject
Cochran-Mantel-Haenszel Statistics (Based on Rank Scores)
Statistic Alternative Hypothesis DF Value Prob
------
1 Nonzero Correlation 1 10.4533 0.0012
2 Row Mean Scores Differ 3 14.2000 0.0026
Total Sample Size = 36
2-Way ANOVA – RBD – Caffeine Study – R Program
pdf("C:\\Rmisc\\graphs\\caffeine1.pdf")
caffeine <- read.fwf("C:\\data\\caffeine.txt", width=c(1,5,8),
col.names=c("subject_c", "dose_c", "endtime"))
attach(caffeine)
dose_c <- factor(dose_c)
subject_c <- factor(subject_c)
caffeine1.rbd <- aov(endtime ~ subject_c + dose_c)
summary(caffeine1.rbd)
TukeyHSD(caffeine1.rbd, "dose_c")
interaction.plot(subject_c,dose_c,endtime)
caffeine.res <- residuals(caffeine1.rbd)
hist(caffeine.res)
friedman.test(endtime ~ dose_c | subject_c)
dev.off()
2-Way ANOVA – RBD – Caffeine Study – R Output
> caffeine1.rbd <- aov(endtime ~ subject_c + dose_c)
> summary(caffeine1.rbd)
Df Sum Sq Mean Sq F value Pr(>F)
subject_c 8 5558.0 694.75 13.2159 4.174e-07 ***
dose_c 3 933.1 311.04 5.9168 0.003591 **
Residuals 24 1261.7 52.57
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> TukeyHSD(caffeine1.rbd, "dose_c")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = endtime ~ subject_c + dose_c)
$dose_c
diff lwr upr p adj
5-0 11.2366667 1.808030 20.665303 0.0153292
9-0 12.2411111 2.812474 21.669748 0.0076616
13-0 11.7088889 2.280252 21.137526 0.0110929
9-5 1.0044444 -8.424192 10.433081 0.9909369
13-5 0.4722222 -8.956414 9.900859 0.9990313
13-9 -0.5322222 -9.960859 8.896414 0.9986162
> interaction.plot(subject_c,dose_c,endtime)
> caffeine.res <- residuals(caffeine1.rbd)
> hist(caffeine.res)
> friedman.test(endtime ~ dose_c | subject_c)
Friedman rank sum test
data: endtime and dose_c and subject_c
Friedman chi-squared = 14.2, df = 3, p-value = 0.002645
2-Factor ANOVA – Lacrosse Study – SAS Program
options nodate nonumber ps=54 ls=76;
data one;
infile 'C:\data\lacrosse.dat';
input brand side gadd;
run;
proc means mean;
class brand;
var gadd;
run;
proc means mean;
class side;
var gadd;
run;
proc means mean std;
class brand side;
var gadd;
run;
proc glm;
class brand side;
model gadd = brand|side;
run;
quit;
2-Factor ANOVA – Lacrosse Study – SAS Output
The MEANS Procedure
Analysis Variable : gadd
N
brand Obs Mean
------
1 20 1070.30
2 20 1069.85
3 20 1116.65
4 20 1359.65
------
side Obs Mean
------
1 40 1090.88
2 40 1217.35
------
brand side Obs Mean Std Dev
------
1 1 10 1166.10 152.3998147
2 10 974.4998000 72.9999420
2 1 10 1117.60 216.1999812
2 10 1022.10 105.0999403
3 1 10 857.0001000 151.4998416
2 10 1376.30 211.4000969
4 1 10 1222.80 123.1000267
2 10 1496.50 183.9998678
------
The GLM Procedure
Class Level Information
Class Levels Values
brand 4 1 2 3 4
side 2 1 2
Number of Observations Read 80
Number of Observations Used 80
Dependent Variable: gadd
Sum of
Source DF Squares Mean Square F Value Pr > F
Model 7 3107552.083 443936.012 17.50 <.0001
Error 72 1826954.058 25374.362
Corrected Total 79 4934506.141
R-Square Coeff Var Root MSE gadd Mean
0.629759 13.80224 159.2933 1154.112
Source DF Type I SS Mean Square F Value Pr > F
brand 3 1155477.556 385159.185 15.18 <.0001
side 1 319918.386 319918.386 12.61 0.0007
brand*side 3 1632156.141 544052.047 21.44 <.0001
Source DF Type III SS Mean Square F Value Pr > F
brand 3 1155477.556 385159.185 15.18 <.0001
side 1 319918.386 319918.386 12.61 0.0007
brand*side 3 1632156.141 544052.047 21.44 <.0001
2-Factor ANOVA – Lacrosse Study – R Program
pdf("lacrosse.pdf")
lacrosse1 <- read.fwf("C:\\data\\lacrosse.dat", width=c(8,8,14),
col.names=c("brand", "side", "gadd"))
attach(lacrosse1)
brand <- factor(brand, levels=1:4,
labels=c("SHC", "SCHAF", "SCHUL", "BUL"))
side <- factor(side, levels=1:2,
labels=c("Front", "Back"))
tapply(gadd, brand, mean) # marginal mean for brand
tapply(gadd, side, mean) # marginal mean for side
tapply(gadd, list(brand,side), mean) # cell means
tapply(gadd, list(brand,side), sd) # cell SDs
# Fit Model with Interaction (A*B says to include main effects and
# interactions)
lacrosse1.aov <- aov(gadd ~ brand*side)
summary(lacrosse1.aov)
# Plot Means of Front and Back versus Brand
interaction.plot(brand,side,gadd)
# Plot residuals versus predicted values
yhat <- predict(lacrosse1.aov)
e <- resid(lacrosse1.aov)
plot(yhat,e)
dev.off()
------
2-Factor ANOVA – Lacrosse Study – R Output
> tapply(gadd, brand, mean) # marginal mean for brand
SHC SCHAF SCHUL BUL
1070.300 1069.850 1116.650 1359.650
> tapply(gadd, side, mean) # marginal mean for side
Front Back
1090.875 1217.350
> tapply(gadd, list(brand,side), mean) # cell means
Front Back
SHC 1166.0999 974.4998
SCHAF 1117.5999 1022.1000
SCHUL 857.0001 1376.3000
BUL 1222.8001 1496.5001
> tapply(gadd, list(brand,side), sd) # cell SDs
Front Back
SHC 152.3998 72.99994
SCHAF 216.2000 105.09994
SCHUL 151.4998 211.40010
BUL 123.1000 183.99987
> # Fit Model with Interaction (A*B says to include main effects and interactions)
> lacrosse1.aov <- aov(gadd ~ brand*side)
> summary(lacrosse1.aov)
Df Sum Sq Mean Sq F value Pr(>F)
brand 3 1155478 385159 15.179 9.459e-08 ***
side 1 319918 319918 12.608 0.0006816 ***
brand:side 3 1632156 544052 21.441 4.988e-10 ***
Residuals 72 1826954 25374
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> # Plot Means of Front and Back versus Brand
> interaction.plot(brand,side,gadd)
> # Plot residuals versus predicted values
> yhat <- predict(lacrosse1.aov)
> e <- resid(lacrosse1.aov)
> plot(yhat,e)
2-Factor Nested ANOVA – Swamp Study – SAS Program
options ps=54 ls=76;
data swamp;
infile 'C:\data\swamp1.dat';
input size swampid watlev;
run;
proc means mean;
class size;
var watlev;
run;
proc means mean std;
class swampid;
var watlev;
run;
proc glm;
class size swampid;
model watlev = size swampid(size);
test h=size e=swampid(size);
run;
quit;
2-Factor Nested ANOVA – Swamp Study – SAS Output
The MEANS Procedure
Analysis Variable : watlev
N
size Obs Mean
------
1 81 52.1666667
2 81 106.4344444
3 81 152.7665432
------
swampid Obs Mean Std Dev
------
1 27 48.0003704 15.1197564
2 27 36.2996296 20.8405574
3 27 72.2000000 20.5796685
4 27 110.2018519 8.1599468
5 27 103.0011111 11.4297881
6 27 106.1003704 15.3802679
7 27 128.8000000 15.3303610
8 27 137.1992593 10.2896905
9 27 192.3003704 14.8099627
------
The GLM Procedure
Class Level Information
Class Levels Values
size 3 1 2 3
swampid 9 1 2 3 4 5 6 7 8 9
Number of Observations Read 243
Number of Observations Used 243
Dependent Variable: watlev
Sum of
Source DF Squares Mean Square F Value Pr > F
Model 8 493781.3279 61722.6660 267.01 <.0001
Error 234 54092.1874 231.1632
Corrected Total 242 547873.5154
R-Square Coeff Var Root MSE watlev Mean
0.901269 14.64897 15.20405 103.7892
Source DF Type I SS Mean Square F Value Pr > F
size 2 410723.7365 205361.8683 888.38 <.0001
swampid(size) 6 83057.5914 13842.9319 59.88 <.0001
Source DF Type III SS Mean Square F Value Pr > F
size 2 410723.7365 205361.8683 888.38 <.0001
swampid(size) 6 83057.5914 13842.9319 59.88 <.0001
Tests of Hypotheses Using the Type III
MS for swampid(size) as an Error Term
Source DF Type III SS Mean Square F Value Pr > F
size 2 410723.7365 205361.8683 14.84 0.0048
2-Factor Nested ANOVA – Swamp Study – R Program
pdf("swamp.pdf")
swamp <- read.fwf("C:\\data\\swamp1.dat", width=c(8,8,8),
col.names=c("size", "swampid", "watlev"))
attach(swamp)
size <- factor(size, levels=1:3,labels=c("Small", "Medium", "Large"))
swampid <- factor(swampid)
tapply(watlev, size, mean)
tapply(watlev, swampid, mean)
tapply(watlev, swampid, sd)
swamp.aov1 <- aov(watlev ~ size + size/swampid)
# This provides ANOVA, not appropriate F-test for fsize
summary(swamp.aov1)
swamp.aov2 <- aov(watlev ~ size + Error(swampid))
# This provides appropriate F-test for fsize
summary(swamp.aov2)
dev.off()
2-Factor Nested ANOVA – Swamp Study – R Output
> tapply(watlev, size, mean)
Small Medium Large
52.16667 106.43444 152.76654
> tapply(watlev, swampid, mean)
1 2 3 4 5 6 7 8
48.00037 36.29963 72.20000 110.20185 103.00111 106.10037 128.80000 137.19926
9
192.30037
> tapply(watlev, swampid, sd)
1 2 3 4 5 6 7 8
15.119756 20.840557 20.579668 8.159947 11.429788 15.380268 15.330361 10.289690
9
14.809963
> swamp.aov1 <- aov(watlev ~ size + size/swampid)
> # This provides ANOVA, not appropriate F-test for fsize
> summary(swamp.aov1)
Df Sum Sq Mean Sq F value Pr(>F)
size 2 410724 205362 888.385 < 2.2e-16 ***
size:swampid 6 83058 13843 59.884 < 2.2e-16 ***
Residuals 234 54092 231
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> swamp.aov2 <- aov(watlev ~ size + Error(swampid))
> # This provides appropriate F-test for fsize
> summary(swamp.aov2)
Error: swampid
Df Sum Sq Mean Sq F value Pr(>F)
size 2 410724 205362 14.835 0.004759 **
Residuals 6 83058 13843
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 234 54092 231.16
Split-Plot ANOVA – Wool Friction Study – SAS Program
options ps=54 ls=76;
data woolfrict;
infile '../../data/woolfrict1.dat';
input run trt revs shrink;
run;
proc glm;
class run trt revs;
model shrink = run trt run*trt revs trt*revs;
test h=run e=run*trt;
test h=trt e=run*trt;
contrast 'Revs Linear' Revs -3 -2 -1 0 1 2 3;
contrast 'Revs Quadratic' Revs 5 0 -3 -4 -3 0 5;
contrast 'Revs Cubic' Revs -1 1 1 0 -1 -1 1;
contrast 'Revs Quartic' Revs 3 -7 1 6 1 -7 3;
contrast 'Revs Quintic' Revs -1 4 -5 0 5 -4 1;
contrast 'Revs Sextic' Revs 1 -6 15 -20 15 -6 1;
contrast 'Trt*Revs Linear'
trt*revs -3 -2 -1 0 1 2 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 2 1 0 -1 -2 -3,
trt*revs 0 0 0 0 0 0 0 -3 -2 -1 0 1 2 3 0 0 0 0 0 0 0 3 2 1 0 -1 -2 -3,
trt*revs 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -3 -2 -1 0 1 2 3 3 2 1 0 -1 -2 -3;
contrast 'Trt*Revs Quadratic'
trt*revs 5 0 -3 -4 -3 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 0 3 4 3 0 -5,
trt*revs 0 0 0 0 0 0 0 5 0 -3 -4 -3 0 5 0 0 0 0 0 0 0 -5 0 3 4 3 0 -5,
trt*revs 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 -3 -4 -3 0 5 -5 0 3 4 3 0 -5;
contrast 'Trt*Revs Cubic'
trt*revs -1 1 1 0 -1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 -1 0 1 1 -1,
trt*revs 0 0 0 0 0 0 0 -1 1 1 0 -1 -1 1 0 0 0 0 0 0 0 1 -1 -1 0 1 1 -1,
trt*revs 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 1 0 -1 -1 1 1 -1 -1 0 1 1 -1;
contrast 'Trt*Revs Quartic'
trt*revs 3 -7 1 6 1 -7 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -3 7 -1 -6 -1 7 -3,
trt*revs 0 0 0 0 0 0 0 3 -7 1 6 1 -7 3 0 0 0 0 0 0 0 -3 7 -1 -6 -1 7 -3,
trt*revs 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 -7 1 6 1 -7 3 -3 7 -1 -6 -1 7 -3;
contrast 'Trt*Revs Quintic'
trt*revs -1 4 -5 0 5 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 5 0 -5 4 -1,
trt*revs 0 0 0 0 0 0 0 -1 4 -5 0 5 -4 1 0 0 0 0 0 0 0 1 -4 5 0 -5 4 -1,
trt*revs 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -5 0 5 -4 1 1 -4 5 0 -5 4 -1;
contrast 'Trt*Revs Sextic'
trt*revs 1 -6 15 -20 15 -6 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 6 -15 20 -15 6 -1,
trt*revs 0 0 0 0 0 0 0 1 -6 15 -20 15 -6 1 0 0 0 0 0 0 0 -1 6 -15 20 -15 6 -1,
trt*revs 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -6 15 -20 15 -6 1 -1 6 -15 20 -15 6 -1;
run;
quit;
Split-Plot ANOVA – Wool Friction Study – SAS Output
The GLM Procedure
Class Level Information
Class Levels Values
run 4 1 2 3 4
trt 4 1 2 3 4
revs 7 200 400 600 800 1000 1200 1400
Number of Observations Read 112
Number of Observations Used 112
Dependent Variable: shrink
Sum of
Source DF Squares Mean Square F Value Pr > F
Model 39 14572.74062 373.66002 271.32 <.0001
Error 72 99.15929 1.37721
Corrected Total 111 14671.89991
R-Square Coeff Var Root MSE shrink Mean
0.993242 5.050423 1.173547 23.23661
Source DF Type I SS Mean Square F Value Pr > F
run 3 124.28598 41.42866 30.08 <.0001
trt 3 3012.52812 1004.17604 729.14 <.0001
run*trt 9 114.63723 12.73747 9.25 <.0001
revs 6 11051.77929 1841.96321 1337.46 <.0001
trt*revs 18 269.51000 14.97278 10.87 <.0001
Source DF Type III SS Mean Square F Value Pr > F
run 3 124.28598 41.42866 30.08 <.0001
trt 3 3012.52813 1004.17604 729.14 <.0001
run*trt 9 114.63723 12.73747 9.25 <.0001
revs 6 11051.77929 1841.96321 1337.46 <.0001
trt*revs 18 269.51000 14.97278 10.87 <.0001
Contrast DF Contrast SS Mean Square F Value Pr > F
Revs Linear 1 10846.82893 10846.82893 7875.93 <.0001
Revs Quadratic 1 198.87574 198.87574 144.40 <.0001
Revs Cubic 1 2.87042 2.87042 2.08 0.1532
Revs Quartic 1 1.43196 1.43196 1.04 0.3113
Revs Quintic 1 0.12190 0.12190 0.09 0.7669
Revs Sextic 1 1.65033 1.65033 1.20 0.2773
Trt*Revs Linear 3 154.36768 51.45589 37.36 <.0001
Trt*Revs Quadratic 3 104.77063 34.92354 25.36 <.0001
Trt*Revs Cubic 3 7.41375 2.47125 1.79 0.1559
Trt*Revs Quartic 3 0.76898 0.25633 0.19 0.9055
Trt*Revs Quintic 3 0.86732 0.28911 0.21 0.8892
Trt*Revs Sextic 3 1.32165 0.44055 0.32 0.8110
Tests of Hypotheses Using the Type III MS for run*trt as an Error Term
Source DF Type III SS Mean Square F Value Pr > F
run 3 124.285982 41.428661 3.25 0.0739
trt 3 3012.528125 1004.176042 78.84 <.0001
Split-Plot ANOVA – Wool Friction Study – R Program
pdf("woolfrict.pdf")
wf1 <- read.fwf("C:\\data\\woolfrict.dat",
width=c(8,8,8,8),
col.names=c("runnum", "trt", "revs", "shrink"))
attach(wf1)
trt <- factor(trt, levels=1:4,
labels=c("Untreated", "AC15sec", "AC4min", "AC15min"))
runnum <- factor(runnum, levels=1:4)
revs <- factor(revs,levels=seq(200,1400,200),ordered=TRUE)
# Fit full ANOVA with all terms, F-test for trt is inappropriate
woolfrict.sp1 <- aov(shrink ~ trt + runnum + trt:runnum + revs + revs:trt)
summary(woolfrict.sp1)
# This model uses correct error terms for trt and revs and interaction
woolfrict.sp2 <- aov(shrink ~ trt*revs + Error(runnum/trt))
summary(woolfrict.sp2)
# Partitions the Revs SS into orthogonal polynomials
summary(woolfrict.sp2,split=list(revs=list(linear=1, quadratic=2,
cubic=3, quartic=4, quintic=5, sextic=6)))
dev.off()
Split-Plot ANOVA – Wool Friction Study – R Output
> woolfrict.sp1 <- aov(shrink ~ trt + runnum + trt:runnum + revs + revs:trt)
> summary(woolfrict.sp1)
Df Sum Sq Mean Sq F value Pr(>F)
trt 3 3012.5 1004.18 729.1367 < 2.2e-16 ***
runnum 3 124.3 41.43 30.0815 1.024e-12 ***
revs 6 11051.8 1841.96 1337.4577 < 2.2e-16 ***
trt:runnum 9 114.6 12.74 9.2487 3.546e-09 ***
trt:revs 18 269.5 14.97 10.8718 4.620e-14 ***