> ## sec 2.3 Bolt Exp't
> #
> data=read.table("http://www.stat.ucla.edu/%7Ehqxu/stat225/data/bolt.dat", h=T, skip=2)  
> data
   P.O C.W HT M.B
1   10  24 32   M
2   13  18 22   M
3   17  17 30   M
4   16  17 35   M
5   15  15 32   M
6   14  23 28   M
7   11  14 27   M
8   14  18 28   M
9   15  12 30   M
10  16  11 30   M
11  25  20 26   B
12  40  16 40   B
13  30  17 28   B
14  17  18 38   B
15  16  15 38   B
16  45  16 30   B
17  49  19 26   B
18  33  14 38   B
19  30  15 45   B
20  20  24 38   B
> 
> # make sure that we have a correct model matrix
> y=c(data$C.W, data$HT, data$P.O)
> test= as.factor(rep(data$M.B, 3))
> plating=as.factor(rep(c("CW", "HT", "PO"), c(20,20,20)))
> 
> # The baseline constraint is known as "contr.treatment" in R.
> options(contrasts=c("contr.treatment", "contr.poly"))
> g=lm(y~test*plating)
> anova(g)  # table 2.10
Analysis of Variance Table

Response: y
             Df Sum Sq Mean Sq F value    Pr(>F)    
test          1  821.4  821.40 22.4563 1.604e-05 ***
plating       2 2290.6 1145.32 31.3118 9.363e-10 ***
test:plating  2  665.1  332.55  9.0916 0.0003952 ***
Residuals    54 1975.2   36.58                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> # Table 2.11 can be read from the summary because the default contrast coding is baseline constraint
> summary(g)

Call:
lm(formula = y ~ test * plating)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.500  -2.525   0.100   2.600  18.500 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       17.400      1.913   9.098 1.74e-12 ***
testM             -0.500      2.705  -0.185 0.854030    
platingHT         17.300      2.705   6.396 3.92e-08 ***
platingPO         13.100      2.705   4.843 1.11e-05 ***
testM:platingHT   -4.800      3.825  -1.255 0.214926    
testM:platingPO  -15.900      3.825  -4.157 0.000116 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.048 on 54 degrees of freedom
Multiple R-squared:  0.6566,	Adjusted R-squared:  0.6248 
F-statistic: 20.65 on 5 and 54 DF,  p-value: 1.806e-11

> 
> ## fit the model again using zero-sum contraints
> options(contrasts=c("contr.sum", "contr.poly"))
> g=lm(y~test*plating)
> summary(g)

Call:
lm(formula = y ~ test * plating)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.500  -2.525   0.100   2.600  18.500 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     23.8333     0.7808  30.525  < 2e-16 ***
test1            3.7000     0.7808   4.739 1.60e-05 ***
plating1        -6.6833     1.1042  -6.053 1.40e-07 ***
plating2         8.2167     1.1042   7.441 7.91e-10 ***
test1:plating1  -3.4500     1.1042  -3.124  0.00286 ** 
test1:plating2  -1.0500     1.1042  -0.951  0.34588    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.048 on 54 degrees of freedom
Multiple R-squared:  0.6566,	Adjusted R-squared:  0.6248 
F-statistic: 20.65 on 5 and 54 DF,  p-value: 1.806e-11

> 
> ## more plots 
> par(mfrow=c(3,2))
> bolt=as.data.frame(cbind(B.CW=y[11:20], B.HT=y[31:40], B.PO=y[51:60], 
+ M.CW=y[1:10], M.HT=y[21:30],M.PO=y[41:50]))
> boxplot(bolt, ylab="Responses", main="Bolt Expriment")
> 
> interaction.plot(test, plating, y, ylab="Mean Response", 
+ main="Interaction Plot")
> interaction.plot(plating, test, y, ylab="Mean Response", 
+ main="Interaction Plot")
> 
> # diagnostics
> plot(g, 1)
> res.df=as.data.frame(cbind(B.CW=g$res[11:20], B.HT=g$res[31:40], B.PO=g$res[51:60], M.CW=g$res[1:10], M.HT=g$res[21:30],M.PO=g$res[41:50]))
> boxplot(res.df, ylab="Residuals", main="Residuals vs. Treatments")
> plot(g, 2)