一元配置分散分析と多重比較 Last modified: Apr 09, 2015
目的
生データに対して，等分散性の検定，一元配置分散分析と多重比較を行う all in one 関数。
使用法
oneway.ANOVA(x, g, verbose=TRUE)
引数
x 分析対象となる観察値のベクトル
g 観察値がどの群に属するかを表すベクトル
verbose デフォルトでは結果をプリントする
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/oneway-ANOVA.R", encoding="euc-jp")
oneway.ANOVA <- function(x, g, verbose=TRUE) { # x: データベクトル，g: 群変数ベクトル
# 基礎統計量
d <- data.frame(x, g)
level <- levels(d$g)
k <- length(level)
n <- aggregate(x ~ g, d, length)[,2]
N <- sum(n)
Mean <- aggregate(x ~ g, d, mean)[,2]
Var <- aggregate(x ~ g, d, var)[,2]
SD <- aggregate(x ~ g, d, sd)[,2]
SE <- SD/sqrt(n)
result <- data.frame(水準=level, n, Mean,
SD, "Mean-SD"=Mean-SD, "Mean+SD"=Mean+SD,
SE, "Mean-SE"=Mean-SE, "Mean+SE"=Mean+SE, check.names=FALSE)
# デザインプロット
if (verbose) {
x.pos <- seq_along(level)
plot(x.pos, Mean, xlim=c(0.5, k+0.5), ylim=range(Mean-SD, Mean+SD), pch=16,
xlab="", xaxt="n", ylab="")
axis(1, at=x.pos, labels=level)
arrows(x.pos, Mean-SD, x.pos, Mean+SD, length=0.10, angle=90, code=3, lty=3)
arrows(x.pos, Mean-SE, x.pos, Mean+SE, length=0.05, angle=90, code=3, lwd=2)
}
# 等分散性の検定
r.bartlett <- bartlett.test(d$x, d$g)
r.Levene <- Levene.test(d$x, d$g)
r.eqvar <- data.frame(手法=c("Bartlett", "Levene"), 分布=c("chi square", "F"),
統計量=c(r.bartlett$statistic, r.Levene$statistic),
自由度=c(r.bartlett$parameter,
sprintf("(%i, %i)", r.Levene$parameter[1], r.Levene$parameter[2])),
"P 値"=c(r.bartlett$p.value, r.Levene$p.value))
# 分散分析（等分散性を仮定する）
r.oneway <- oneway.test(x ~ g, d, var.equal=TRUE)
gm <- mean(d$x)
St <- sum((d$x-gm)^2)
Sb <- sum(n*(Mean-gm)^2)
Sw <- St-Sb
SS <- c(Sb, Sw, St)
df <- c(k-1, N-k, N-1)
MS <- SS/df
F <- c(MS[1]/MS[2], NA, NA)
p <- c(pf(F[1], df[1], df[2], lower.tail=FALSE), NA, NA)
MS[3] <- NA
r.ANOVA <- data.frame(要因=c("級間", "級内", "全体"), 平方和=c(Sb, Sw, St),
自由度=df, 平均平方=MS, "F 値"=F, "P 値"=p, check.names=FALSE)
# 分散分析（等分散性を仮定しない）
r.Welch <- oneway.test(x ~ g, d)
r.BF <- Brown.Forsythe.test(x, g)
r.ANOVA2 <- data.frame(手法=c("Welch", "Brown-Forsythe"),
"F 値"=c(r.Welch$statistic, r.BF$statistic),
"第 1 自由度"=c(k-1, r.BF$parameter[1]),
"第 2 自由度"=c(r.Welch$parameter[2], r.BF$parameter[2]),
"P 値"=c(r.Welch$p.value, r.BF$p.value), check.names=FALSE)
# 多重比較
r.tukey <- tukey(d$x, d$g)$Tukey
lc <- t(sapply(rownames(r.tukey), function(s) as.integer(unlist(strsplit(s, ":")))))
diff <- Mean[lc[,2]]-Mean[lc[,1]]
r.GH <- tukey(d$x, d$g, method="G")$Games.Howell
r.multi <- data.frame(群1=level[lc[,1]], 群2=level[lc[,2]],
平均1=Mean[lc[,1]], 平均2=Mean[lc[,2]], 差=diff,
差の標準誤差=diff/r.tukey[,1], "t 値(Tukey)"=r.tukey[,1],
"P 値"=r.tukey[,2], "t 値(Games-Howell)"=r.GH[,1],
"自由度"=r.GH[,2], "P 値"=r.GH[,3], check.names=FALSE)
if (verbose) {
cat("\n基礎統計量\n\n")
print(result, row.names=FALSE)
cat("\n等分散性の検定\n\n")
print(r.eqvar, row.names=FALSE)
cat("\n分散分析（等分散性を仮定する）\n\n")
print(r.ANOVA, row.names=FALSE)
cat("\n分散分析（等分散性を仮定しない）\n\n")
print(r.ANOVA2, row.names=FALSE)
cat("\n多重比較\n\n")
print(r.multi, row.names=FALSE)
}
invisible(list(result, r.eqvar, r.ANOVA, r.ANOVA2, r.multi))
}
他に，以下の 3 つの関数も必要
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/Brown-Forsythe-test.R", encoding="euc-jp")
Brown.Forsythe.test <- function(x, group) { # x: データベクトル，group: 群変数ベクトル
data.name <- paste(deparse(substitute(x)), "~", deparse(substitute(group)))
OK <- complete.cases(x, group) # 欠損値を持つケースを除く
x <- x[OK]
group <- as.factor(group[OK])
group <- group[, drop=TRUE]
d <- split(x, group)
df1 <- length(d)-1
ni <- sapply(d, length)
mi <- sapply(d, mean)
ui <- sapply(d, var)
ci <- (1-ni/sum(ni))*ui
F.BF <- sum(ni*(mi-mean(x))^2)/sum(ci)
C <- ci/sum(ci)
df2 <- 1/sum(C^2/(ni-1))
p <- pf(F.BF, df1, df2, lower.tail=FALSE)
method <- "Brown-Forsythe 検定"
return(structure(list(statistic=c(F=F.BF), parameter=c("df1"=df1, "df2"=df2), "p.value"=p, method=method, data.name=data.name), class="htest"))
}
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/levene-test.R", encoding="euc-jp")
# 等分散性の検定
Levene.test <- function(x, # データベクトル
group, # 群変数ベクトル
method = c("mean", "median")) # 検定統計量の計算方法
{
data.name <- paste(deparse(substitute(x)), "~", deparse(substitute(group)))
OK <- complete.cases(x, group) # 欠損値を持つケースを除く
x <- x[OK]
fac <- as.factor(group[OK])
fac <- fac[, drop=TRUE]
n <- length(x) # 全体のデータ個数
n.i <- tapply(x, fac, length) # 各群のデータ個数
k <- length(n.i) # 群の数
method <- match.arg(method) # 引数の補完
x <- abs(x-tapply(x, fac, method)[fac]) # 測定値からそのデータが属する群の平均値（または中央値）を差し引く
sw <- sum((n.i-1)*tapply(x, fac, var)) # 群内変動
dfw <- n-k # 群内変動の自由度
dfb <- k-1 # 群間変動の自由度
f <- ((var(x)*(n-1)-sw)/dfb)/(sw/dfw) # 検定統計量
P <- pf(f, dfb, dfw, lower.tail=FALSE) # P 値
method <- paste("等分散性の検定（", method, " で調整）", sep="")
return(structure(list(statistic=c(F=f), parameter=c("df b"=dfb, "df w"=dfw), p.value=P,
method=method, data.name=data.name), class="htest"))
}
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/tukey.R", encoding="euc-jp")
# Tukey の方法による多重比較
# Games-Howell の方法も選択できるように拡張 2009/08/03
tukey <- function( data, # 観察値ベクトル
group, # 群変数ベクトル
method=c("Tukey", "Games-Howell")) # 手法の選択
{
OK <- complete.cases(data, group) # 欠損値を持つケースを除く
data <- data[OK]
group <- factor(group[OK])
n <- tapply(data, group, length) # 各群のケース数
a <- length(n) # 群の数
phi.e <- sum(n)-a # 誤差分散（群内不偏分散）の自由度
Mean <- tapply(data, group, mean) # 各群の平均値
Variance <- tapply(data, group, var) # 各群の不偏分散
result1 <- cbind(n, Mean, Variance) # 各群の統計量
rownames(result1) <- paste("Group", 1:a, sep="")
method <- match.arg(method)
if (method == "Tukey") { # Tukey の方法
v.e <- sum((n-1)*Variance)/phi.e # 誤差分散（群内不偏分散）
t <- combn(a, 2, function(ij) # 対比較
abs(diff(Mean[ij]))/sqrt(v.e*sum(1/n[ij])) )
p <- ptukey(t*sqrt(2), a, phi.e, lower.tail=FALSE) # 有意確率を計算する
Tukey <- cbind(t, p) # 対比較の結果
rownames(Tukey) <- combn(a, 2, paste, collapse=":")
return(list(result1=result1, Tukey=Tukey, phi=phi.e, v=v.e))
}
else { # Games-Howell の方法
t.df <- combn(a, 2, function(ij) { # 対比較
t <- abs(diff(Mean[ij]))/sqrt(sum(Variance[ij]/n[ij]))
df <- sum(Variance[ij]/n[ij])^2/sum((Variance[ij]/n[ij])^2/(n[ij]-1))
return(c(t, df))} )
t <- t.df[1,]
df <- t.df[2,]
p <- ptukey(t*sqrt(2), a, df, lower.tail=FALSE) # 有意確率を計算する
Games.Howell <- cbind(t, df, p) # 対比較の結果
rownames(Games.Howell) <- combn(a, 2, paste, collapse=":")
return(list(result1=result1, Games.Howell=Games.Howell))
}
}
使用例
> x <- c(
205, 206, 164, 190, 194, 203,
201, 221, 197, 185,
248, 265, 197, 220, 212, 281,
202, 276, 237, 254, 230
)
> g <- rep(paste("A", 1:4, sep=""), c(6, 4, 6, 5))
> oneway.ANOVA(x, g)
基礎統計量
水準 n Mean SD Mean-SD Mean+SD SE Mean-SE Mean+SE
A1 6 -6.333333 15.88290 -22.21624 9.549572 6.484169 -12.817502 0.1508354
A2 4 1.000000 14.96663 -13.96663 15.966630 7.483315 -6.483315 8.4833148
A3 6 37.166667 32.72563 4.44104 69.892294 13.360181 23.806485 50.5268480
A4 5 39.800000 27.58985 12.21015 67.389853 12.338557 27.461443 52.1385575
等分散性の検定
手法 分布 統計量 自由度 P.値
Bartlett chi square 3.31855 3 0.3450692
Levene F 2.38566 (3, 17) 0.1050173
分散分析（等分散性を仮定する）
要因 平方和 自由度 平均平方 F 値 P 値
級間 9284.271 3 3094.7571 5.091555 0.01072138
級内 10332.967 17 607.8216 NA NA
全体 19617.238 20 NA NA NA
分散分析（等分散性を仮定しない）
手法 F 値 第 1 自由度 第 2 自由度 P 値
Welch 4.927451 3 8.860454 0.02766496
Brown-Forsythe 5.440656 3 13.322350 0.01163382
多重比較
群1 群2 平均1 平均2 差 差の標準誤差 t 値(Tukey) P 値 t 値(Games-Howell) 自由度 P 値
A1 A2 -6.333333 1.00000 7.333333 15.91411 0.4608069 0.96656131 0.7406106 6.871719 0.87779880
A1 A3 -6.333333 37.16667 43.500000 14.23402 3.0560597 0.03276031 2.9291841 7.231682 0.07993964
A1 A4 -6.333333 39.80000 46.133333 14.92876 3.0902319 0.03058705 3.3097548 6.139828 0.05810207
A2 A3 1.000000 37.16667 36.166667 15.91411 2.2726160 0.14395624 2.3617950 7.413345 0.16830990
A2 A4 1.000000 39.80000 38.800000 16.53843 2.3460508 0.12655164 2.6887452 6.340133 0.12046051
A3 A4 37.166667 39.80000 2.633333 14.92876 0.1763933 0.99797118 0.1447992 8.990836 0.99883370
直前のページへ戻る E-mail to Shigenobu AOKI
