目的

分散・共分散行列の同等性の検定（Box の M 検定）を行う

使用法

BoxM(x, gvar)

引数

x      データ行列（データフレーム）
gvar   データ行列の各行がどのグループに属するかを表すベクトル（factor 化して使われる）

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/BoxM.R", encoding="euc-jp")

# 分散・共分散行列の同等性の検定
BoxM <- function(    x,                                              # データ行列（データフレーム）
                        gvar)                                           # 群を表す変数
{
        method <- "分散・共分散行列の同等性の検定"
        data.name <- paste(deparse(substitute(x)), "~", deparse(substitute(gvar)))
        x <- as.data.frame(x)
        nv <- ncol(x)                                                        # 変数の個数
        if (nv < 2) stop("分散共分散行列を計算する変数は2個以上必要です")
        gvar <- as.factor(gvar)
        ni <- table(gvar)                                            # 各群のサンプルサイズ
        n <- length(gvar)                                            # サンプルサイズ
        g <- length(ni)                                                      # 群の数
        y <- split(x, gvar)                                          # 群ごとに分割したデータ行列
        Si <- lapply(y, var)                                         # 分散共分散行列
        log.det.Si <- sapply(Si, function(x) log(det(x)))            # 行列式の対数値
        S <- sapply(y, function(x) (nrow(x)-1)*var(x))                       # 変動行列
        S <- matrix(rowSums(S), nv, nv)/(n-g)                                # プールした変動行列
        M <- (n-g)*log(det(S))-sum((ni-1)*log.det.Si)                        # Box の M 統計量
        f1 <- (g-1)*nv*(nv+1)/2                                              # 第 1 自由度
        rho <- 1-(2*nv^2+3*nv-1)/(6*(nv+1)*(g-1))*(sum(1/(ni-1))-1/(n-g))
        tau <- (nv-1)*(nv+2)/(6*(g-1))*(sum(1/(ni-1)^2)-1/(n-g)^2)
        f2 <- (f1+2)/abs(tau-(1-rho)^2)                                      # 第 2 自由度
        gamma <- (rho-f1/f2)/f1
        F <- M*gamma                                                 # F 値
        p <- pf(F, f1, f2, lower.tail=FALSE)                         # P 値
        return(structure(list(statistic=c(M=M, F=F),
                parameter=c(df1=f1, df2=f2), p.value=p,
                method=method, data.name=data.name), class="htest"))
}


使用例

> x <- matrix(c( # 1~3 列がデータ，4 列目はグループ
+ 	2.9, 161.7, 120.8, 1,
+ 	2.3, 114.8,  85.2, 1,
+ 	2.0, 128.4,  92.0, 1,
+ 	3.2, 149.2,  97.3, 1,
+ 	2.7, 126.0,  81.1, 1,
+ 	4.4, 133.8, 107.6, 1,
+ 	4.1, 161.3, 114.0, 1,
+ 	2.1, 111.5,  77.3, 1,
+ 	4.8, 198.7, 172.9, 2,
+ 	3.6, 199.3, 157.9, 2,
+ 	2.0, 188.4, 152.7, 2,
+ 	4.9, 183.6, 164.2, 2,
+ 	3.9, 173.5, 172.2, 2,
+ 	4.4, 184.9, 163.2, 2
+ ), byrow=TRUE, ncol=4)
> 
> BoxM(x[,1:3], x[,4])

	分散・共分散行列の同等性の検定

data:  x[, 1:3] ~ x[, 4]
M = 9.7304, F = 1.1535, df1 = 6.000, df2 = 795.195, p-value = 0.3294

 解説ページ（2群の場合）