> ## sec 2.9 Starch Exp't
> #
> dat<-read.table("http://www2.isye.gatech.edu/%7Ejeffwu/book/data/starch.dat", h=T)
> dat
   starch strength thickness
1      CA    791.7       7.7
2      CA    610.0       6.3
3      CA    710.0       8.6
...
12     CA    990.0      12.5
13     CA    862.7      11.7
14     CO    731.0       8.0
15     CO    710.0       7.3
...
32     CO    592.5       7.2
33     PO    983.3      13.0
34     PO    958.8      13.3
...
48     PO   1050.0      14.1
49     PO    973.3      13.7
> g<- lm( strength ~ thickness + as.factor(starch), data = dat)
> summary(g) # table 2.28

Call:
lm(formula = strength ~ thickness + as.factor(starch), data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-203.63  -99.45  -57.84   56.72  637.61 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           158.26     179.78   0.880 0.383360    
thickness              62.50      17.06   3.664 0.000653 ***
as.factor(starch)CO   -83.67      86.10  -0.972 0.336351    
as.factor(starch)PO    70.36      67.78   1.038 0.304795    
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 164.7 on 45 degrees of freedom
Multiple R-Squared: 0.6815,	Adjusted R-squared: 0.6602 
F-statistic: 32.09 on 3 and 45 DF,  p-value: 3.001e-11 

> anova(g) # table 2.29
Analysis of Variance Table

Response: strength
                  Df  Sum Sq Mean Sq F value    Pr(>F)    
thickness          1 2553357 2553357 94.1859 1.318e-12 ***
as.factor(starch)  2   56725   28362  1.0462    0.3597    
Residuals         45 1219940   27110                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
> 
> g<- lm( strength ~ as.factor(starch), data = dat)
> anova(g) # table 2.30
Analysis of Variance Table

Response: strength
                  Df  Sum Sq Mean Sq F value    Pr(>F)    
as.factor(starch)  2 2246204 1123102  32.619 1.512e-09 ***
Residuals         46 1583818   34431                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
>

# some plots
pdf("starch2.pdf", w=8, h=4.5)
par(mfrow=c(1,2))
plot(strength ~ starch , data=dat)
plot(strength ~ thickness, data=dat, col=rep(c(1,2,3), c(13,19,17)) ) 
dev.off()