> dat=read.table("http://www2.isye.gatech.edu/%7Ejeffwu/book/data/WaterResistance.dat", h=F, skip=1)
> dat  # 
   V1 V2 V3 V4   V5
1   4  2  2  1 53.5
2   4  2  4  1 32.5
3   4  2  1  1 46.6
4   4  2  3  1 35.4
5   5  2  4  2 44.6
6   5  2  1  2 52.2
7   5  2  3  2 45.9
8   5  2  2  2 48.3
9   1  1  3  1 40.8
10  1  1  1  1 43.0
11  1  1  2  1 51.8
12  1  1  4  1 45.5
13  2  1  2  2 60.9
14  2  1  4  2 55.3
15  2  1  3  2 51.1
16  2  1  1  2 57.4
17  6  2  1  3 32.1
18  6  2  4  3 30.1
19  6  2  2  3 34.4
20  6  2  3  3 32.2
21  3  1  1  3 52.8
22  3  1  3  3 51.7
23  3  1  4  3 55.3
24  3  1  2  3 59.2
> names(dat)=c("Whole Plot", "A","B","Repl","Resistance")  # rename the variables
> 
> A=as.factor(dat[,2])
> B=as.factor(dat[,3])
> Repl=as.factor(dat[,4])
> Resistance=dat[,5]
> 
> g=lm(Resistance~(A+B+Repl)^3)
> anova(g)
Analysis of Variance Table

Response: Resistance
          Df Sum Sq Mean Sq F value Pr(>F)
A          1 782.04  782.04               
B          3 266.01   88.67               
Repl       2 376.99  188.49               
A:B        3  62.79   20.93               
A:Repl     2 398.38  199.19               
B:Repl     6  95.65   15.94               
A:B:Repl   6  56.87    9.48               
Residuals  0   0.00                       
> 
> g=lm(Resistance~A+Repl+A:Repl+B+A:B)  
> anova(g)  # Correct SS and MS, but wrong F and p values for A and replicate,  you can adjust these manually
Analysis of Variance Table

Response: Resistance
          Df Sum Sq Mean Sq F value    Pr(>F)    
A          1 782.04  782.04 61.5303 4.596e-06 ***
Repl       2 376.99  188.49 14.8304 0.0005711 ***
B          3 266.01   88.67  6.9763 0.0056928 ** 
A:Repl     2 398.38  199.19 15.6719 0.0004503 ***
A:B        3  62.79   20.93  1.6468 0.2309105    
Residuals 12 152.52   12.71                      
---
> 
> # F and p for A
> 782.04/199.19  
[1] 3.926101
> 1-pf(782.04/199.19, 1,2) 
[1] 0.1860529
> # F and p for Replicate
> 188.49/199.19
[1] 0.9462824
> 1-pf(188.49/199.19, 2,3)
[1] 0.4801501