目的

多変量に拡張された平均値の差の検定（ウィルクスのΛ）を行う
rrcov の Wilks.test は，method != "mcd" のときにはいつの場合にも「それ以外の場合」のカイ二乗近似を使うようだ

使用法

wilks(dat, g)

引数

dat	データ行列
g	群分け変数

ソース

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

# 多変量に拡張された平均値の差の検定（ウィルクスのΛ）
wilks <- function(   dat,                                    # データ行列
                        g)                                      # 群変数
{
        method <- "多変量に拡張された平均値の差の検定（ウィルクスのΛ）"
        data.name <- paste(deparse(substitute(dat)), "~", deparse(substitute(g)))
        OK <- complete.cases(dat, g)                         # 欠損値を持つケースを除く
        dat <- dat[OK,]
        g <- as.factor(g[OK])
        nv <- ncol(dat)                                              # 変数の個数
        g.case <- table(g)                                   # 各群のケース数
        n <- sum(g.case)                                     # 全ケース数
        k <- length(g.case)                                  # 群の個数

        vars <- by(dat, g, function(x) var(x)*(nrow(x)-1))   # 各群の変動・共変動行列
        w <- matrix(rowSums(mapply("+", vars)), nv)          # 群内変動・共変動行列
        t <- (var(dat)*(n-1))                                        # 全変動・共変動行列

        LAMBDA <- det(w)/det(t)
        sl <- sqrt(LAMBDA)
        if (nv == 2) {                                          # 2 変数データの場合
                f <- (1-sl)*(n-k-1)/sl/(k-1)
                p <- pf(f, df1 <- 2*(k-1), df2 <- 2*(n-k-1), lower.tail=FALSE)
                result <- list(statistic=c("F value"=f), parameter=c("df1"=df1, "df2"=df2))
        }
        else if (k == 2) {                                      # 2 群の場合
                f <- (1-LAMBDA)*(n-k-nv+1)/LAMBDA/nv
                p <- pf(f, nv, df2 <- n-k-nv+1, lower.tail=FALSE)
                result <- list(statistic=c("F value"=f), parameter=c("df1"=nv, "df2"=df2))
        }
        else if (k == 3) {                                      # 3 群の場合
                f <- (1-sl)*(n-k-nv+1)/sl/nv
                p <- pf(f, df1 <- 2*nv, df2 <- 2*(n-k-nv+1), lower.tail=FALSE)
                result <- list(statistic=c("F value"=f), parameter=c("df1"=df1, "df2"=df2))
        }
        else {                                                  # それ以外の場合
                chi.sq <- ((nv+k)/2-n+1)*log(LAMBDA)
                p <- pchisq(chi.sq, df <- nv*(k-1), lower.tail=FALSE)
                result <- list(statistic=c("chi sq. value"=chi.sq), parameter=c("df"=df))
        }
        return(structure(c(result, p.value=p, method=method, data.name=data.name), class="htest"))
}


使用例

> dat <- matrix(c(
+ 	1, 2.9, 161.7, 120.8,
+ 	1, 2.3, 114.8, 85.2,
+ 	1, 2, 128.4, 92,
+ 	1, 3.2, 149.2, 97.3,
+ 	1, 2.7, 126, 81.1,
+ 	1, 4.4, 133.8, 107.6,
+ 	1, 4.1, 161.3, 114,
+ 	1, 2.1, 111.5, 77.3,
+ 	2, 4.8, 198.7, 172.9,
+ 	2, 3.6, 199.3, 157.9,
+ 	2, 2, 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
+ ), byrow=TRUE, ncol=4)
> 
> wilks(dat[, 2:4], dat[, 1])

	多変量に拡張された平均値の差の検定（ウィルクスのΛ）

data:  dat[, 2:4] ~ dat[, 1] 
F value = 33.805, df1 = 3, df2 = 10, p-value = 1.517e-05

 解説ページ