############################################################################## # dae.R # Montgomery's Design and Analysis of Experimentsの例を実行するためのスクリプト # 作成日時: <2011-11-30> ############################################################################## ############################################################################## # 第5章 (pp. 162--206) # ページは，Design and Analysis of Experimentsの第7版のもの ############################################################################## #バッテリーデータの読み込み battery <- read.table("battery.txt", header=TRUE) battery$Material <- as.factor(battery$Material) battery$Temperature <- as.factor(battery$Temperature) summary(battery) # 分散分析 battery.aov <- aov(Life ~ Material * Temperature, data=battery) # 交互作用プロット with(battery, interaction.plot(Temperature, Material, Life,type="b", pch=19, fixed=T, xlab="温度（F）", ylab="平均寿命")) # 効果プロット plot.design(Life ~ Material * Temperature, data=battery) # 交互作用が有意の場合の平均の対比を求める正しい方法 # 温度70Fで3つの平均を計算する mm <- with(subset(battery, Temperature==70), aggregate(Life, list(M=Material), mean)) # 次に，t分布のパーセント点と誤差（ANOVAからプールした標準偏差SDを用いて） # の積を求める． val.crit <- qtukey(.95, 3, 27) * sqrt(unlist(summary(battery.aov))[["Mean Sq4"]]/4) # 最後に，観測された平均の差と限界値とを比較する diff.mm <- c(d.3.1=mm$x[3]-mm$x[1], d.3.2=mm$x[3]-mm$x[2], d.2.1=mm$x[2]-mm$x[1]) names(which(diff.mm > val.crit)) # 交互作用をプールしたモデルの作成 summary(battery.aov2 <- aov(Life ~ Material + Temperature, data=battery)) # 化学製品不純物データ impurity <- read.table("impurity.txt", header=TRUE) impurity$Temperature <- as.factor(impurity$Temperature) impurity$Pressure <- as.factor(impurity$Pressure) summary(aov(N ~ Temperature * Pressure, impurity)) # 残差プロット library(alr3) residualPlots(lm(N ~ Temperature + Pressure, impurity)) # ボトル充填データ bottling <- read.table("bottling.txt", header=TRUE, colClasses=c("numeric", rep("factor", 3))) summary(bottling) opar <- par(mfrow=c(2, 2), cex=.8) boxplot(Deviation ~ ., data=bottling, las=2, cex.axis=.8, ylab="偏差") abline(h=0, lty=2) par(las=1) mm <- with(bottling, tapply(Deviation, Carbonation, mean)) ss <- with(bottling, tapply(Deviation, Carbonation, sd)) bp <- barplot(mm,xlab="炭酸化", ylab="偏差", ylim=c(-2,9)) arrows(bp, mm - ss / sqrt(4), bp, mm + ss / sqrt(4), code=3, angle=90, length=.1) with(bottling, interaction.plot(Carbonation, Pressure, Deviation, type="b")) with(bottling, interaction.plot(Carbonation, Speed, Deviation, type="b")) par(opar) # 3元配置分散分析 summary(bottling.aov <- aov(Deviation ~ . ^3, bottling)) # 交互作用のプーリング bottling.aov2 <- aov(Deviation ~ ., bottling) anova(bottling.aov2, bottling.aov) # バッテリー寿命データ分析で，2次項の追加 battery$Temperature.num <- as.numeric(as.character(battery$Temperature)) battery.aov3 <- aov(Life ~ Material + Temperature.num + I(Temperature.num^2) + Material:Temperature.num+Material:I(Temperature.num^2), data=battery) summary(battery.aov3) # 2次項を含むモデルの推定結果のグラフ化 new <- data.frame(Temperature.num=rep(seq(15, 125, by=5), 3), Material=gl(3, 23)) new$fit <- predict(battery.aov3, new) opar <- par(las=1) # 最初に予測値をプロットする with(new, interaction.plot(Temperature.num, Material,fit, legend=FALSE, xlab="温度", ylab="寿命", ylim=c(20, 190))) txt.leg <- paste("原材料の種類", 1:3) text(5:7,new$fit[new$Temperature.num==c(45,55,65)] - c(3, 3, -20), txt.leg, pos=1) # 次に観測値を points(rep(c(1, 15, 23), each=12), battery$Life, pch=19) par(opar) # 工具寿命データ tool <- read.table("toollife.txt", header=TRUE) tool.lm <- lm(Life ~ Angle * Speed + I(Angle^2) * I(Speed^2) + Angle:I(Speed^2) + I(Angle^2):Speed, tool) tmp.angle <- seq(15, 25, by=.1) tmp.speed <- seq(125, 175, by=.5) tmp <- list(Angle=tmp.angle, Speed=tmp.speed) new <- expand.grid(tmp) new$fit <- c(predict(tool.lm,new)) # 等高線プロット require(lattice) contourplot(fit ~ Angle * Speed,data=new, cuts=8, region=T, col.regions=gray(7:16/16)) # レーダー強度データ intensity <- read.table("intensity.txt", header=TRUE, colClasses=c("numeric", rep("factor", 3))) require(lattice) xyplot(Intensity ~ Ground|Operator, data=intensity, groups=Filter, panel=function(x, y,...){ subs <- list(...)$subscripts panel.xyplot(x, y, pch=c(1,19),...) panel.superpose(x, y, panel.groups="panel.lmline", lty=c(1,2),...) }, key=list(text=list(lab=as.character(1:2)), lines=list(lty=1:2, col=1), corner=c(1, .95), title="フィルタのタイプ", cex.title=.8), col=1) # Error()の利用 intensity.aov <- aov(Intensity ~ Ground * Filter + Error(Operator), intensity) summary(intensity.aov) ############################################################################## # 第6章 (pp. 207--272) # ページは，Design and Analysis of Experimentsの第7版のもの ############################################################################## # 収量データ yield <- read.table("yield.txt", header=T) attach(yield) rm(yield) yield.sums <- aggregate(yield, list(reactant=reactant, catalyst=catalyst), sum) yield.sums # 回帰モデル reactant.num <- reactant levels(reactant.num) <- c(25, 15) reactant.num <- as.numeric(as.character(reactant.num)) catalyst.num <- catalyst levels(catalyst.num) <- c(2, 1) catalyst.num <- as.numeric(as.character(catalyst.num)) yield.lm <- lm(yield ~ reactant.num + catalyst.num) yield.lm # 線型モデルLMの係数 # 応答曲面と等高線図の作成 library(scatterplot3d) s3d <- scatterplot3d(reactant.num, catalyst.num, yield, type="n", angle=135, scale.y=1, xlab="反応剤", ylab="触媒", zlab="収量") s3d$plane3d(yield.lm,lty.box="solid", col="darkgray") tmp <- list(reactant.num=seq(15, 25, by=.5), catalyst.num=seq(1, 2, by=.1)) new.data <- expand.grid(tmp) new.data$fit <- predict(yield.lm, new.data) contourplot(fit ~ reactant.num + catalyst.num,new.data, xlab="反応剤", ylab="触媒") # エッチングデータの読み込み plasma <- read.table("plasma.txt", header=TRUE) plasma # データの変換 plasma.df <- data.frame(etch=c(plasma$R1, plasma$R2), rbind(plasma[,2:4], plasma[,2:4])) plasma.df[,2:4] <- lapply(plasma.df[,2:4], factor) # 分散分析 plasma.df.aov <- aov(etch ~ A * B * C, data=plasma.df) summary(plasma.df.aov)