install.packages(“DoE.base”)
library(DoE.base)
library(lattice)
library(latticeExtra)
# =============================================================================(
plot.summaly <- function(a, df, df.mean, …){
x.lab <- df.mean[[a]]$label
# xy plot
plot.default(df[[a]], df$x, xaxt = “n”, ann = FALSE, …)
par(new = TRUE)
plot.default(df.mean[[a]]$label, df.mean[[a]]$x,
ylab = “x”, xlab = a,
xaxt = “n”,
type = “l”, col = “red”,
…)
axis(side = 1, at = x.lab, labels = levels(x.lab))
# boxplot
boxplot(formula(sprintf(“x~%s”, a)), data = df, …)
}
plot.Interaction <- function(a, x.df, …){
f00 <- function(a, x.list, y, add.flg = FALSE, …){
f01 <- function(x, y, a, data.list, …){
list.name <- names(data.list)
legend.lab <- list.name
pchs = rep(1, length(legend.lab))
ltys = rep(“b”, length(legend.lab))
for(i in list.name){
order.z <- order(data.list[[i]][[x]])
df <- data.list[[i]][order.z,]
x.lab <- df[[x]]
main.lab <- paste(a, x, sep = ” – “)
if(rev(list.name)[1] == i){
plot.default(df[[x]], df[[y]], col = df[[a]],
ylab = y, xlab = x, main = main.lab,
xaxt = “n”,
…)
} else {
plot.default(df[[x]], df[[y]], col = df[[a]],
xaxt = “n”, ann = FALSE, …)
par(new = add.flg)
}
if(!add.flg){
axis(side = 1, at = x.lab, labels = levels(x.lab))
legend(“topleft”, legend = legend.lab, col = (1:length(legend.lab)),
pch = pchs, #lty = ltys,
bg = “transparent”)
}
}
if(add.flg){
# if(levels(x.lab)[1] == “D1”){
# browser()
# }
x.lab <- df[[x]]
axis(side = 1, at = x.lab, labels = levels(x.lab))
legend(“topleft”, legend = legend.lab, col = (1:length(legend.lab)),
pch = pchs, #lty = ltys,
bg = “transparent”)
}
TRUE
}
x.list2 <- x.list[[a]]
tg.x <- names(x.list2[[1]])
z <- grep(paste(c(a, y), collapse = “|”), tg.x)
lapply(tg.x[-z], f01, y, a, x.list2, …)
}
s.df <- lapply(a, function(a, x){split(x, x[[a]])}, x.df)
names(s.df) <- a
lapply(a, f00, s.df, “x”, ylim = c(0, 0.5), type = “b”, …)
}
cal.LSD <- function(tg.f, data.aov, a = 0.05, tg.x = “x”){
f.lsd <- function(n.comb, n, v.value, t.value){
((1/n[n.comb][1]) + (1/n[n.comb][2]) * v.value)^(1/2) * t.value
}
f.delta <- function(n.comb, x){
rtn <- data.frame(x[n.comb][1] – x[n.comb][2])
names(rtn) <- paste(names(x)[n.comb][1], names(x)[n.comb][2], sep = ” “)
row.names(rtn) <- NULL
rtn
}
aov.sum <- summary(data.aov)[[1]]
tg.data <- data.aov$model[[tg.f]]
if(is.null(tg.data)){
print(sprintf(“error : 対象の要素がNULLです。(tg.f : %s)”, tg.f))
return(NA)
}
lv <- levels(tg.data)
lv.comb <- combn(length(lv),2)
z <- length(row.names(aov.sum)[[1]])
n <-table(tg.data)
e.tValue <- abs(qt(a/2, aov.sum[[“Df”]][z]))
x.mean <- tapply(data.aov$model[[tg.x]], data.aov$model[[tg.f]], mean)
x.lsd <- apply(lv.comb, 2, f.lsd, n, aov.sum[[“Mean Sq”]][z], e.tValue)
names(x.lsd) <- names(x.mean)
x.delta.0 <- apply(lv.comb, 2,
function(n.comb, x){x[n.comb][1] – x[n.comb][2]},
x.mean)
names(x.delta.0) <- apply(lv.comb, 2,
function(n.comb, x){paste(x[n.comb][1], x[n.comb][2], sep = ” – “)},
names(x.mean))
x.delta <- as.data.frame(apply(lv.comb, 2, f.delta, x.mean))
# as.data.frame後にnamesに含まれる“-“が“.”に置き換わるのでgsub()再置換する。
names(x.delta) <- gsub(“[.]”, ” – “, names(x.delta))
# ————————————————————————-
print(“(Residuals)”)
print(sprintf(“factor : %s”, paste(lv, collapse = “, “)))
#print(sprintf(“factor : %s”, lv))
print(” – mean -“)
print(x.mean)
print(sprintf(“t-Value : %s”, round(e.tValue, 3)))
print(” – L.S.D -“)
print(x.lsd)
print(” – delta – “)
print(x.delta)
#print(x.delta.0)
rtn <- list()
rtn$factor <- lv
rtn$n <- n
rtn$mean <- x.mean
rtn$lsd <- x.lsd
rtn$delta <- x.delta
rtn
}
# =============================================================================
old.wd <- getwd()
setwd(“/Documents/R”)
# L9(3^4)の直行表を読み込む
show.oas(nruns=9, parents.only=FALSE)
oaTableL9 <- oa.design(ID=L9.3.4, randomize=FALSE)
print(oaTableL9)
x <- read.csv(“L9_data.csv”,
header = TRUE,
stringsAsFactors = FALSE)
x.data.info <- read.csv(“L9_data_factor_levels.csv”,
header = TRUE,
stringsAsFactors = FALSE)
# C因子について実験データの割付を直行表にあわせる。
data.frame(“直行表割付水準“ = as.integer(oaTableL9.org$C),
“実験データの水準” = c(1, 2, 3, 2, 3, 1, 3, 1, 2))
# 2と3を入れ替える。
x.data.info$C <- x.data.info$C[c(1,3,2)]
x.data.info$C
x.data <- cbind(oaTableL9, x = x$x)
# x.data.2<- data.frame(
# A = c(sprintf(“A%s”, as.character(x.data$A))),
# B = c(sprintf(“B%s”, as.character(x.data$B))),
# C = c(sprintf(“C%s”, as.character(x.data$C))),
# D = c(sprintf(“D%s”, as.character(x.data$D))),
# x = x.data$x)
x.data.header <- names(x.data)
z <- (x.data.header != “x”)
x.data.2 <- data.frame(
as.data.frame(
lapply(x.data.header[z],
function(a, x){paste(a, x[[a]], sep = “”)},
x.data)),
x.data$x)
x.data.2.mean <- lapply(x.data.header[z],
function(a, x){x <- tapply(x$x, x[[a]], mean); data.frame(x.mean = c(x), label = names(x))},
x.data.2)
names(x.data.2.mean) <- x.data.header[z]
print(x.data.2.mean)
#
old.par <- par(mfcol=c(2,length(x.data.header[z])))
lapply(x.data.header[z],
plot.summaly,
x.data.2, x.data.2.mean, ylim = c(0, 0.5))
par(old.par)
# 交互作用確認(Interaction)
#
# —————————————————————————
# old.par <- par(mfcol=c(3,3))
# plot.Interaction(“A”, x.data.2, add.flg = TRUE)
# par(old.par)
# for nextを利用してplotするとparによる複数レイアウトができないので
# layout()を使ってplotを配置する。
# layout(matrix(seq(3), rep(1, 3)), 1, 3)
p.mat <- matrix(seq(12), ncol = 3, nrow = 4, byrow = TRUE)
layout(p.mat,
widths = rep.int(1, ncol(p.mat)),
heights = rep.int(1, nrow(p.mat)))
plot.Interaction(x.data.header[-5], x.data.2, add.flg = TRUE)
# —————————————————————————
# 分散分析
result.aov <- aov(x ~ A + B + C, data = x.data)
summary(result.aov)
tapply(x.data$x, x.data$A, mean)
result.2.aov <- aov(x ~ A + B + C, data = x.data.2)
summary(result.2.aov)
# 誤差分散よりも値が小さいBの要素をプーリング
result.2.aov.2 <- aov(x ~ A + C, data = x.data.2)
summary(result.2.aov.2)
# 最小有意差:LSD(Least Significant Difference)を求める
#cal.LSD(“A”, result.2.aov, “A1”, “A2”)
result.lsd <- lapply(x.data.header[c(1,3)], cal.LSD, result.2.aov.2)
# lapply(x.data.header[z][4], cal.LSD, result.2.aov)
# debug(cal.LSD)
names(result.lsd) <- x.data.header[c(1,3)]
result.lm <- lm(x ~ A + B + C, data = x.data)
summary.lm(result.lm)
result.lm.2 <- lm(x ~ A + B + C, data = x.data.2)
summary.lm(result.lm.2)
Teedamari 