东莞网站建设智搜宝网站建设北京海淀

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

dataframe<-data.frame(

strength=c(9.60,9.69,8.43,9.98,11.28,10.10,11.01,10.44,9.00,9.57,9.03,9.80),

Temperature=gl(3,4,12),

pressure=gl(4,1,12))

summary (dataframe)

dataframe.aov2 <- aov(strength~pressure+Temperature,data=dataframe)

summary (dataframe.aov2)

> summary (dataframe.aov2)

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

pressure     3  0.581  0.1936   0.539 0.6727 

Temperature  2  4.658  2.3288   6.487 0.0316 *

Residuals    6  2.154  0.3590                

---

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

with(dataframe,interaction.plot(Temperature,pressure,strength,type="b",pch=19,fixed=T,xlab="Temperature (°F)",ylab="strength"))

plot.design(strength~pressure+Temperature,data=dataframe)

fit <-lm(strength~pressure+Temperature,data=dataframe)

anova(fit)

> anova(fit)

Analysis of Variance Table

Response: strength

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

pressure     3 0.5807 0.19356  0.5392 0.67270 

Temperature  2 4.6576 2.32882  6.4873 0.03162 *

Residuals    6 2.1539 0.35898                 

---

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

summary(fit)

> summary(fit)

Call:

lm(formula = strength ~ pressure + Temperature, data = dataframe)

Residuals:

    Min      1Q  Median      3Q     Max

-0.6575 -0.4902  0.1233  0.3067  0.6400

Coefficients:

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

(Intercept)    9.5575     0.4237  22.559 4.97e-07 ***

pressure2     -0.1733     0.4892  -0.354   0.7352   

pressure3     -0.4700     0.4892  -0.961   0.3738   

pressure4      0.1133     0.4892   0.232   0.8245   

Temperature2   1.2825     0.4237   3.027   0.0232 * 

Temperature3  -0.0750     0.4237  -0.177   0.8653   

---

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

Residual standard error: 0.5991 on 6 degrees of freedom

Multiple R-squared:  0.7086,    Adjusted R-squared:  0.4658

F-statistic: 2.918 on 5 and 6 DF,  p-value: 0.1124

par(mfrow=c(2,2))

plot(fit)

par(mfrow=c(2,2))

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

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