二元配置分散分析     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)) {                       # 繰り返し数が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) 各水準の繰り返しが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

Made with Macintosh