本文是实验设计与分析（第6版，Montgomery著，傅珏生译) 第5章析因设计引导5.7节思考题5.7 R语言解题。主要涉及方差分析，正态假设检验，残差分析，交互作用图，等值线图。

dataframe <-data.frame(

force=c(2.70,2.78,2.83,2.86,2.45,2.49,2.85,2.80,2.60,2.72,2.86,2.87,2.75,2.86,2.94,2.88),

feed=gl(4,4,16),

speed=gl(2,2,16))

summary (dataframe)

dataframe.aov2 <- aov(force~feed*speed,data=dataframe)

summary (dataframe.aov2)

> summary (dataframe.aov2)

            Df  Sum Sq Mean Sq F value   Pr(>F)   

feed         3 0.09250 0.03083  11.859  0.00258 **

speed        1 0.14822 0.14822  57.010 6.61e-05 ***

feed:speed   3 0.04187 0.01396   5.369  0.02557 * 

Residuals    8 0.02080 0.00260                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

with(dataframe,interaction.plot(feed,speed,force,type="b",pch=19,fixed=T,xlab="feed",ylab="force"))

plot.design(force~feed*speed,data=dataframe)

fit <-lm(force~feed*speed,data=dataframe)

anova(fit)

> anova(fit)

Analysis of Variance Table

Response: force

           Df   Sum Sq  Mean Sq F value    Pr(>F)   

feed        3 0.092500 0.030833 11.8590  0.002582 **

speed       1 0.148225 0.148225 57.0096 6.605e-05 ***

feed:speed  3 0.041875 0.013958  5.3686  0.025567 * 

Residuals   8 0.020800 0.002600                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(fit)

> summary(fit)

Call:

lm(formula = force ~ feed * speed, data = dataframe)

Residuals:

     Min       1Q   Median       3Q      Max

-0.06000 -0.02625  0.00000  0.02625  0.06000

Coefficients:

               Estimate Std. Error t value Pr(>|t|)   

(Intercept)   2.740e+00  3.606e-02  75.994    1e-12 ***

feed2        -2.700e-01  5.099e-02  -5.295 0.000733 ***

feed3        -8.000e-02  5.099e-02  -1.569 0.155303   

feed4         6.500e-02  5.099e-02   1.275 0.238172   

speed2        1.050e-01  5.099e-02   2.059 0.073449 . 

feed2:speed2  2.500e-01  7.211e-02   3.467 0.008482 **

feed3:speed2  1.000e-01  7.211e-02   1.387 0.202934   

feed4:speed2 -5.912e-16  7.211e-02   0.000 1.000000   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05099 on 8 degrees of freedom

Multiple R-squared:  0.9314,    Adjusted R-squared:  0.8715

F-statistic: 15.53 on 7 and 8 DF,  p-value: 0.0004502

par(mfrow=c(2,2))

plot(fit)

par(mfrow=c(2,2))

plot(as.numeric(dataframe$feed), fit$residuals, xlab="feed", ylab="Residuals", type="p", pch=16)

plot(as.numeric(dataframe$speed), fit$residuals, xlab="speed", ylab="Residuals", pch=16)

dataframe<-data.frame(

force=c(2.70,2.78,2.83,2.86,2.45,2.49,2.85,2.80,2.60,2.72,2.86,2.87,2.75,2.86,2.94,2.88),

feed=c(0.015,0.015,0.015,0.015,0.030,0.030,0.030,0.030,0.045,0.045,0.045,0.045,0.060,0.060,0.060,0.060),

speed=c(125,125,200,200,125,125,200,200,125,125,200,200,125,125,200,200))

fit <-lm(force~feed*speed+feed*I(speed^2)+I(feed^2)*speed+I(feed^2)+I(speed^2),data=dataframe)

tmp.speed <- seq(125,200,by=.5)

tmp.feed <- seq(0.015,0.060,by=.005)

tmp <- list(feed=tmp.feed,speed=tmp.speed)

new <- expand.grid(tmp)

new$fit <- c(predict(fit,new))

require(lattice)

contourplot (fit~feed*speed ,data=new, cuts=8,region=T,col.regions=gray(7:16/16))