本文是實驗設計與分析（第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))