本文是实行计划与分析（第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))


免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。