> ## Ranitidine Experiment, table 9.2
> dat=read.table("table9_2.dat", h=F)
> dat
   V1    V2    V3    V4        V5     V6
1   1 -1.00 -1.00 -1.00    17.293  2.850
2   2  1.00 -1.00 -1.00    45.488  3.817
3   3 -1.00  1.00 -1.00    10.311  2.333
4   4  1.00  1.00 -1.00 11757.084  9.372
5   5 -1.00 -1.00  1.00    16.942  2.830
6   6  1.00 -1.00  1.00    25.400  3.235
7   7 -1.00  1.00  1.00 31697.199 10.364
8   8  1.00  1.00  1.00 12039.201  9.396
9   9  0.00  0.00 -1.67     7.474  2.011
10 10  0.00  0.00  1.67     6.312  1.842
11 11  0.00 -1.68  0.00    11.145  2.411
12 12  0.00  1.68  0.00     6.664  1.897
13 13 -1.68  0.00  0.00 16548.749  9.714
14 14  1.68  0.00  0.00 26351.811 10.179
15 15  0.00  0.00  0.00     9.854  2.288
16 16  0.00  0.00  0.00     9.606  2.262
17 17  0.00  0.00  0.00     8.863  2.182
18 18  0.00  0.00  0.00     8.783  2.173
19 19  0.00  0.00  0.00     8.013  2.081
20 20  0.00  0.00  0.00     8.059  2.087
> A=dat[,2]
> B=dat[,3]
> C=dat[,4]
> CEF=dat[,5]
> lnCEF=dat[,6]
> # g=lm(log(CEF) ~ (A+B+C)^2 +I(A^2) +I(B^2)+I(C^2), subset=c(1:6, 8:20))
> g=lm(lnCEF ~ (A+B+C)^2 +I(A^2) +I(B^2)+I(C^2), subset=-7)  # model (9.5) using all runs except run 7
> summary(g)  # Table 9.11 (the results are slightly different due to rounding)

Call:
lm(formula = lnCEF ~ (A + B + C)^2 + I(A^2) + I(B^2) + I(C^2), 
    subset = -7)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.541024 -0.168735 -0.003401  0.155146  1.748843 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.18540    0.57808   3.780  0.00435 ** 
A            1.11752    0.42429   2.634  0.02719 *  
B            0.79483    0.42429   1.873  0.09380 .  
C            0.01080    0.42552   0.025  0.98031    
I(A^2)       2.71301    0.37955   7.148 5.38e-05 ***
I(B^2)      -0.04794    0.37955  -0.126  0.90226    
I(C^2)      -0.13054    0.38337  -0.340  0.74129    
A:B          1.46539    0.58863   2.490  0.03445 *  
A:C         -0.19311    0.58863  -0.328  0.75037    
B:C          0.20411    0.58863   0.347  0.73675    
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 1.417 on 9 degrees of freedom
Multiple R-Squared: 0.8938,	Adjusted R-squared: 0.7877 
F-statistic:  8.42 on 9 and 9 DF,  p-value: 0.001997 

> 
> # to get the exact results as in Table 9.11, use the following coding
> A[13:14] = c((2-5.5)/(5.5-3.42), (9-5.5)/(7.58-5.5))
> B[11:12] = c((9.9-20)/(20-14), (30.1-20)/(26-20))
> C[9:10] = c((0-5)/(5-2), (10-5)/(8-5))
> g=lm(lnCEF ~ (A+B+C)^2 +I(A^2) +I(B^2)+I(C^2), subset=-7)  # model (9.5) using all runs except run 7
> summary(g)  # Table 9.11 

Call:
lm(formula = lnCEF ~ (A + B + C)^2 + I(A^2) + I(B^2) + I(C^2), 
    subset = -7)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.547573 -0.163802 -0.003029  0.152761  1.746257 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.18503    0.57852   3.777  0.00437 ** 
A            1.11691    0.42424   2.633  0.02724 *  
B            0.79263    0.42416   1.869  0.09450 .  
C            0.01007    0.42623   0.024  0.98166    
I(A^2)       2.70608    0.37882   7.143  5.4e-05 ***
I(B^2)      -0.04598    0.37858  -0.121  0.90600    
I(C^2)      -0.12952    0.38496  -0.336  0.74424    
A:B          1.46672    0.58901   2.490  0.03441 *  
A:C         -0.19178    0.58901  -0.326  0.75218    
B:C          0.20278    0.58901   0.344  0.73855    
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 1.418 on 9 degrees of freedom
Multiple R-Squared: 0.8937,	Adjusted R-squared: 0.7874 
F-statistic: 8.406 on 9 and 9 DF,  p-value: 0.002009 

> ## Fig. 9.9, contour plots with significant terms only
> g=lm(lnCEF ~ (A+B)^2 +I(A^2), subset=-7)  # using all runs except run 7
> g

Call:
lm(formula = lnCEF ~ (A + B)^2 + I(A^2), subset = -7)

Coefficients:
(Intercept)            A            B       I(A^2)          A:B  
     2.0373       1.1543       0.7552       2.7103       1.5306  

> x=seq(-2, 2, .01)
> y=seq(-2, 2, .01)
> f=function(x, y)  2.0373 + 1.1543*x + 0.7552*y + 2.7103*x^2 + 1.5306*x*y 
> z=outer(x, y, f)
> contour(x*(5.5-3.42)+5.5, y*(20-14)+20, z, xlim=c(1.34, 9.66))  # Fig. 9.9 (cannot match exactly)