二元配置分散分析 Last modified: Dec 13, 2004
目的
二元配置分散分析を行う
使用法
twoway.anova(x, a, b)
引数
x 測定値
a 要因 A を表す自然数（1から始まる数）
b 要因 B を表す自然数（1から始まる数）
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/twoway_anova.R", encoding="euc-jp")
# 二元配置分散分析
twoway.anova <- function( x, # データベクトル
a, # 要因 A の factor ベクトル
b) # 要因 B の factor ベクトル
{
name.a <- deparse(substitute(a))
name.b <- deparse(substitute(b))
name.ab <- paste(name.a, name.b, sep="*")
OK <- complete.cases(x, a, b)
x <- x[OK]
a <- a[OK]
b <- b[OK]
tbl <- table(a, b) # 各水準の組合せの繰り返し数
n <- length(x) # 全データ数
dft <- n-1 # 自由度
na <- nrow(tbl) # 要因 A の水準数
dfa <- na-1 # 自由度
nb <- ncol(tbl) # 要因 B の水準数
dfb <- nb-1 # 自由度
gm <- mean(x) # 全データの平均値
MSt <- var(x) # 不偏分散
sst <- dft*MSt # 全変動は不偏分散に自由度を掛けたもの
ma <- tapply(x, a, mean) # 要因 A の水準ごとの平均値
mb <- tapply(x, b, mean) # 要因 B の水準ごとの平均値
hirei <- function(tbl) # 繰り返し数が周辺度数に比例するなら TRUE を返す関数
{
for (i in 2:nrow(tbl)) {
temp <- tbl[i,]/tbl[1,] # 繰り返し数のクロス集計表の，2 行目以降について，1行目のベクトルとの比を計算
if (any(temp != temp[1])) { # どれか一つでも比が違えば，比例していない
return(FALSE)
}
}
return(TRUE)
}
if (any(tbl == 0)) {
result <- "いくつかの水準の組合せにおいて，データがないセルがあります"
}
else if (all(tbl == 1)) { # 繰り返し数が１の場合
ssa <- nb*sum((ma-gm)^2) # 要因 A による変動
ssb <- na*sum((mb-gm)^2) # 要因 B による変動
sse <- sst-ssa-ssb # 誤差変動
dfe <- dft-dfa-dfb # 自由度
ss <- c(ssa, ssb, sse, sst) # 平方和
df <- c(dfa, dfb, dfe, dft) # 自由度
ms <- ss/df # 平均平方
f <- ms/ms[3] # F 値
p <- pf(f, df, df[3], lower.tail=FALSE) # P 値
f[3:4] <- p[3:4] <- NA # 計算しないセル
result <- cbind(ss, df, ms, f, p)
colnames(result) <- c("SS", "d.f.", "MS", "F value", "P value")
rownames(result) <- c(name.a, name.b, "e", "T")
}
else if(all(tbl >= 2) && hirei(tbl)) { # 繰り返し数が周辺度数に比例する場合
mab <- matrix(0, na, nb) # 水準の組み合わせごとの平均値を計算する
for (i in 1:na) {
for (j in 1:nb) {
mab[i, j] <- mean(x[a == i & b == j])
}
}
sse <- 0 # 残差平方和
for (i in 1:na) {
for (j in 1:nb) {
sse <- sse+sum((x[a == i & b == j]-mab[i,j])^2)
}
}
ssa <- sum(rowSums(tbl)*(ma-gm)^2) # 要因 A による平方和
ssb <- sum(colSums(tbl)*(mb-gm)^2) # 要因 B による平方和
ss <- c(ssa, ssb, sst-ssa-ssb-sse, sse, sst)
dfab <- dfa*dfb
df <- c(dfa, dfb, dfab, dft-dfa-dfb-dfab, dft)
ms <- ss/df # 平均平方
f <- ms/ms[4] # F 値
p <- pf(f, df, df[4], lower.tail=FALSE) # P 値
f[4:5] <- p[4:5] <- NA # 計算の不要なセル
result <- cbind(ss, df, ms, f, p) # 母数モデル
f[1:2] <- ms[1:2]/ms[3] # F 値
p[1:2] <- pf(f[1:2], df[1:2], df[3], lower.tail=FALSE) # P 値
result2 <- cbind(ss, df, ms, f, p) # 変量モデル
colnames(result) <- colnames(result2) <- c("SS", "d.f.", "MS", "F value", "P value")
rownames(result) <- rownames(result2) <- c(name.a, name.b, name.ab, "e", "T")
result <- list(Model1=result, Model2=result2)
}
else { # 一般の場合
mreg <- function(dat) # 回帰分析を行い，回帰変動，誤差変動，全変動を返す関数
{
nc <- ncol(dat)
ans <- lm(dat[,nc] ~ dat[, -nc])
St <- var(dat[, nc])*(nrow(dat)-1)
Se <- sum(ans$residuals^2)
return(c(St-Se, Se, St)) # 回帰変動，誤差変動，全変動
}
d1 <- model.matrix(~factor(a))[,-1] # 要因 A に関するダミー変数
d2 <- model.matrix(~factor(b))[,-1] # 要因 B に関するダミー変数
d3 <- model.matrix(~factor(a)*factor(b))[,-1] # 要因 A，B の交互作用も含むダミー変数
r.a.b.ab <- mreg(cbind(d3, x)) # フルモデル
r.a.b <- mreg(cbind(d1, d2, x)) # 要因 A，B の主効果のみを含むモデル
r.a <- mreg(cbind(d1, x)) # 要因 A を含むモデル
r.b <- mreg(cbind(d2, x)) # 要因 B を含むモデル
ss <- c(r.a.b[1]-r.b[1], r.a.b[1]-r.a[1], r.a.b.ab[1]-r.a.b[1], r.a.b.ab[2], r.a.b.ab[3]) # 変動の構成
df <- c(dfa, dfb, dfa*dfb, n-na*nb, n-1)
ms <- ss/df # 平均平方
f <- ms/ms[4] # F 値
p <- pf(f, df, df[4], lower.tail=FALSE) # P 値
f[4:5] <- p[4:5] <- NA # 計算不要のセル
result <- cbind(ss, df, ms, f, p)
colnames(result) <- c("SS", "d.f.", "MS", "F value", "P value")
rownames(result) <- c(name.a, name.b, name.ab, "e", "T")
}
return(result)
}
使用例
(1) 各水準の繰り返しが１の場合
> a <- rep(1:3, each=4)
> b <- rep(1:4, 3)
> x <- c(9, 17, 12, 16, 1, 21, 16, 11, 7, 19, 6, 9)
> data.frame(x, a, b)
x a b
1 9 1 1
2 17 1 2
3 12 1 3
4 16 1 4
5 1 2 1
6 21 2 2
7 16 2 3
8 11 2 4
9 7 3 1
10 19 3 2
11 6 3 3
12 9 3 4
> twoway.anova(x, a, b)
SS d.f. MS F value P value # 母数モデル
a 21.50000 2 10.75000 0.6592845 0.55102999
b 268.66667 3 89.55556 5.4923339 0.03719245
e 97.83333 6 16.30556 NA NA # NA はこの場合，値を表示する必要のない箇所であることを示しているだけ
T 388.00000 11 35.27273 NA NA
aov 関数を使う場合
# 上の場合だと，引き続いて
> df1 <- data.frame(f1=as.factor(a), f2=as.factor(b),x=x)
> summary(aov(x ~ f1+f2, df1))
Df Sum Sq Mean Sq F value Pr(>F)
f1 2 21.500 10.750 0.6593 0.55103
f2 3 268.667 89.556 5.4923 0.03719 *
Residuals 6 97.833 16.306
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
(2) 各水準の繰り返しが周辺度数に比例する場合（全て等しい場合を含む）
> x <- c(
24.8, 23.9, 24.1, 28.8, 22.6, 28.0, 26.4, 27.4, 29.4, 30.0, 28.7,
29.2, 25.0, 26.6, 27.9, 28.5, 27.1, 25.2, 27.9, 29.2, 26.7, 25.0,
29.9, 29.4, 27.5, 32.5, 29.5, 26.3, 28.2, 31.8, 30.2, 31.7, 29.2,
30.3, 30.9, 28.4, 29.8, 26.7, 30.7, 28.5, 26.5, 26.7, 31.7, 27.2,
25.5, 29.9, 27.3, 28.8, 28.5, 25.7, 28.7, 31.3, 29.4, 29.8, 30.3,
29.6, 31.7, 33.6, 32.0, 34.3)
> a <- rep(rep(1:4, each=3), 5)
> b <- rep(1:5, each=12)
> twoway.anova(x, a, b)
$Model1
SS d.f. MS F value P value # 母数モデル
a 51.39733 3 17.132444 4.937067 0.0052039202
b 106.28733 4 26.571833 7.657221 0.0001114609
a*b 52.85267 12 4.404389 1.269215 0.2738926131
e 138.80667 40 3.470167 NA NA # NA はこの場合，値を表示する必要のない箇所であることを示しているだけ
T 349.34400 59 5.921085 NA NA
$Model2
SS d.f. MS F value P value # 変量モデル
a 51.39733 3 17.132444 3.889857 0.037390836
b 106.28733 4 26.571833 6.033035 0.006719395
a*b 52.85267 12 4.404389 1.269215 0.273892613
e 138.80667 40 3.470167 NA NA
T 349.34400 59 5.921085 NA NA
aov 関数を使う場合（母数モデル）
# 上の場合は，引き続いて
> df2 <- data.frame(f1=as.factor(a), f2=as.factor(b),x=x)
> summary(aov(x ~ f1+f2+f1:f2, df2))
Df Sum Sq Mean Sq F value Pr(>F)
f1 3 51.397 17.132 4.9371 0.0052039 **
f2 4 106.287 26.572 7.6572 0.0001115 ***
f1:f2 12 52.853 4.404 1.2692 0.2738926
Residuals 40 138.807 3.470
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> dat3 <- matrix(c(
1, 1, 7.0, 1, 1, 5.7,
1, 2, 4.6, 1, 2, 5.6,
1, 3, 7.0, 1, 3, 7.4,
1, 3, 6.6, 1, 4, 7.3,
1, 4, 6.3, 1, 4, 7.8,
1, 5, 5.3, 1, 5, 4.7,
2, 1, 5.8, 2, 1, 8.9,
2, 2, 10.0, 2, 2, 7.1,
2, 3, 6.5, 2, 3, 9.8,
2, 3, 5.0, 2, 4, 9.7,
2, 4, 6.4, 2, 4, 9.9,
2, 5, 7.8, 2, 5, 5.4,
3, 1, 8.1, 3, 1, 6.1, 3, 1, 10.2, 3, 1, 6.9,
3, 2, 5.2, 3, 2, 5.7, 3, 2, 11.8, 3, 2, 6.1,
3, 3, 10.4, 3, 3, 5.6, 3, 3, 9.4, 3, 3, 5.9, 3, 3, 12.1, 3, 3, 5.7,
3, 4, 5.3, 3, 4, 9.7, 3, 4, 5.5, 3, 4, 10.2, 3, 4, 13.8, 3, 4, 7.0,
3, 5, 11.1, 3, 5, 4.2, 3, 5, 7.1, 3, 5, 5
), byrow=TRUE, ncol=3)
> a <- dat3[,1]
> b <- dat3[,2]
> x <- dat3[,3]
> twoway.anova(x, a, b)
$Model1
SS d.f. MS F value P value # 母数モデル
a 20.796875 2 10.3984375 1.7589914 0.1879828
b 19.536458 4 4.8841146 0.8261929 0.5180374
a*b 7.848958 8 0.9811198 0.1659654 0.9940031
e 195.082500 33 5.9115909 NA NA # NA はこの場合，値を表示する必要のない箇所であることを示しているだけ
T 243.264792 47 5.1758466 NA NA
$Model2
SS d.f. MS F value P value # 変量モデル
a 20.796875 2 10.3984375 10.5985401 0.005636406
b 19.536458 4 4.8841146 4.9781022 0.026003536
a*b 7.848958 8 0.9811198 0.1659654 0.994003145
e 195.082500 33 5.9115909 NA NA
T 243.264792 47 5.1758466 NA NA
aov 関数を用いる場合（母数モデル）
# 上に引き続いて
> df3 <- data.frame(f1=as.factor(a), f2=as.factor(b),x=x)
> summary(aov(x ~ f1+f2+f1:f2, df3))
Df Sum Sq Mean Sq F value Pr(>F)
f1 2 20.797 10.398 1.7590 0.188
f2 4 19.536 4.884 0.8262 0.518
f1:f2 8 7.849 0.981 0.1660 0.994
Residuals 33 195.082 5.912
(3) 一般の場合（要因が直交していない場合）
> a <- c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> b <- c(1, 1, 2, 2, 1, 1, 2, 2, 2, 2)
> x <- c(17, 16, 25, 22, 18, 26, 34, 30, 34, 30)
> twoway.anova(x, a, b)
SS d.f. MS F value P value # 母数モデル
a 121.440476 1 121.440476 13.7479784 0.009995282
b 177.190476 1 177.190476 20.0592992 0.004197947
a*b 5.142857 1 5.142857 0.5822102 0.474368685
e 53.000000 6 8.833333 NA NA # NA はこの場合，値を表示する必要のない箇所であることを示しているだけ
T 415.600000 9 46.177778 NA NA
aov 関数を使う場合は，要因を指定する順序で結果が変わる
# 上に引き続いて
> df4 <- data.frame(f1=as.factor(a), f2=as.factor(b),x=x)
> summary(aov(x ~ f2+f1 * f2, df4))
Df Sum Sq Mean Sq F value Pr(>F)
f2 1 236.017 236.017 26.7189 0.002076 **
f1 1 121.440 121.440 13.7480 0.009995 ** # twoway.anova(x, a, b) の要因 a の行と同じ結果
f2:f1 1 5.143 5.143 0.5822 0.474369
Residuals 6 53.000 8.833
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> summary(aov(x ~ f1+f1 * f2, df4))
Df Sum Sq Mean Sq F value Pr(>F)
f1 1 180.267 180.267 20.4075 0.004028 **
f2 1 177.190 177.190 20.0593 0.004198 ** # twoway.anova(x, a, b) の要因 b の行と同じ結果
f1:f2 1 5.143 5.143 0.5822 0.474369
Residuals 6 53.000 8.833
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
以上の結果より，「aov 関数は Type I の平方和を使っている」ことがわかる。
要因の指定順序によらない結果を得るためには，Type II の平方和を使う。car パッケージの Anova 関数により，以下のように行う。
> library(car)
> Anova(aov(x ~ f1+f1 * f2, df4))
Anova Table (Type II tests)
Response: x
Sum Sq Df F value Pr(>F)
f1 121.440 1 13.7480 0.009995
f2 177.190 1 20.0593 0.004198
f1:f2 5.143 1 0.5822 0.474369
Residuals 53.000 6
> Anova(aov(x ~ f2+f1 * f2, df4))
Anova Table (Type II tests)
Response: x
Sum Sq Df F value Pr(>F)
f2 177.190 1 20.0593 0.004198
f1 121.440 1 13.7480 0.009995
f2:f1 5.143 1 0.5822 0.474369
Residuals 53.000 6
解説ページ
直前のページへ戻る E-mail to Shigenobu AOKI
