Two Way layout
We will work with Torque data
First let's read data in:
> data<-read.table("/afs/stat.lsa.umich.edu/users/student/kutsyy/.public_html/classes/data/bolt.dat",header=T)
> data[1:3,]
torque plating locknut
1 10
P&O mandrel
2 13
P&O mandrel
3 17
P&O mandrel
> attach(data)
> torque
1 2 3 4 5
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26
10 13 17 16 15 14 11 14 15 16 25 40
30 17 16 45 49 33 30 20 24 18 17 17 15 23
27 28 29 30 31 32 33 34 35 36 37 38
39 40 41 42 43 44 45 46 47 48 49 50 51 52
14 18 12 11 20 16 17 18 15 16 19 14
15 24 32 22 30 35 32 28 27 28 30 30 26 40
53 54 55 56 57 58 59 60
28 38 38 30 26 38 45 38
Note that I used attach comand.
Now let's fir linear model.
> fit<-lm(torque~plating*locknut)
> summary(fit)
Call: lm(formula = torque ~ plating * locknut)
Residuals:
Min
1Q Median 3Q Max
-14.5 -2.525 0.1
2.6 18.5
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept)
17.4000 1.9125 9.0979 0.0000
platingHT 17.3000 2.7047 6.3962
0.0000
platingP&O
13.1000 2.7047 4.8434 0.0000
locknut -0.5000 2.7047 -0.1849
0.8540
platingHTlocknut -4.8000
3.8251 -1.2549 0.2149
platingP&Olocknut -15.9000
3.8251 -4.1568 0.0001
Residual standard error: 6.048 on 54 degrees
of freedom
Multiple R-Squared: 0.6566
F-statistic: 20.65 on 5 and 54 degrees of
freedom, the p-value is 1.806e-11
Correlation of Coefficients:
(Intercept) platingHT platingP&O locknut platingHTlocknut
platingHT -0.7071
platingP&O
-0.7071 0.5000
locknut -0.7071 0.5000
0.5000
platingHTlocknut 0.5000
-0.7071 -0.3536 -0.7071
platingP&Olocknut 0.5000
-0.3536 -0.7071 -0.7071 0.5000
> model.matrix(fit)
(Intercept) platingHT platingP&O
locknut platingHTlocknut platingP&Olocknut
1
1 0
1 1
0
1
2
1 0
1 1
0
1
3
1 0
1 1
0
1
4
1 0
1 1
0
1
5
1 0
1 1
0
1
6
1 0
1 1
0
1
7
1 0
1 1
0
1
8
1 0
1 1
0
1
9
1 0
1 1
0
1
10
1 0
1 1
0
1
11
1 0
1 0
0
0
12
1 0
1 0
0
0
13
1 0
1 0
0
0
14
1 0
1 0
0
0
15
1 0
1 0
0
0
16
1 0
1 0
0
0
17
1 0
1 0
0
0
18
1 0
1 0
0
0
19
1 0
1 0
0
0
20
1 0
1 0
0
0
21
1 0
0 1
0
0
22
1 0
0 1
0
0
23
1 0
0 1
0
0
24
1 0
0 1
0
0
25
1 0
0 1
0
0
26
1 0
0 1
0
0
27
1 0
0 1
0
0
28
1 0
0 1
0
0
29
1 0
0 1
0
0
30
1 0
0 1
0
0
31
1 0
0 0
0
0
32
1 0
0 0
0
0
33
1 0
0 0
0
0
34
1 0
0 0
0
0
35
1 0
0 0
0
0
36
1 0
0 0
0
0
37
1 0
0 0
0
0
38
1 0
0 0
0
0
39
1 0
0 0
0
0
40
1 0
0 0
0
0
41
1 1
0 1
1
0
42
1 1
0 1
1
0
43
1 1
0 1
1
0
44
1 1
0 1
1
0
45
1 1
0 1
1
0
46
1 1
0 1
1
0
47
1 1
0 1
1
0
(Intercept) platingHT platingP&O
locknut platingHTlocknut platingP&Olocknut
48
1 1
0 1
1
0
49
1 1
0 1
1
0
50
1 1
0 1
1
0
51
1 1
0 0
0
0
52
1 1
0 0
0
0
53
1 1
0 0
0
0
54
1 1
0 0
0
0
55
1 1
0 0
0
0
56
1 1
0 0
0
0
57
1 1
0 0
0
0
58
1 1
0 0
0
0
59
1 1
0 0
0
0
60
1 1
0 0
0
0
Is your model matrix looks the same? If not,
what is wrong wioth your computer (or my)? How can you fix it?
Now let's get ANOVA table:
> anova(fit)
Analysis of Variance Table
Response: torque
Terms added sequentially (first to last)
Df Sum of Sq Mean Sq F Value
Pr(F)
plating
2 2290.633 1145.317 31.31182 0.0000000009
locknut
1 821.400 821.400 22.45626 0.0000160406
plating:locknut 2 665.100
332.550 9.09159 0.0003952398
Residuals
54 1975.200 36.578
Now let's take a look at plots:
> motif()
> par(mfrow=c(2,3))
> plot(fit)
Do you know what is each plot is for?
SAS
Now let's do the same in SAS.
data bolt;
infile "/afs/stat.lsa.umich.edu/users/student/kutsyy/.public_html/classes/data/bolt1.dat";
input torque 9-10 plating $ 17-19 locknut
$ 26-32;
run;
proc print data=bolt;
run;
proc glm data=bolt;
class plating locknut;
model torque= plating
locknut plating*locknut;
output out=outANOVA predicted=P
residual=R;
means torque/ bon tukey
scheffe;
run;
proc plot;
plot R*P;
run;
|