餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

尖度で自由度を補正したF検定による等分散性の検証

R言語

教科書の分散比のF検定をつかった等分散性の検証は、F検定が正規分布に依存しているので利用できる場合が限られるわけですが、それは分布の裾野の厚さ、尖度が問題であって、尖度が分かればそれで自由度を補正してより広い分布に使えることが知られています。

理屈

O'Neill (2014)のp.294の S^2_n/\sigma^2 \sim \chi^2(DF_n)から、 \frac{S^2_m/\sigma^2}{S^2_n/\sigma^2} \sim \frac{\chi^2(DF_m)}{\chi^2(DF_n)}=F(DF_m,DF_n)であることが直ちに分かります。ここで、 n,mはサブサンプルの観測数で、 DF_n, DF_mは尖度による補正をした自由度です。 DF_n = 2n / (\kappa - (n-3)/(n-1))であり、 DF_mも同様です。 \sigma, \kappaは帰無仮説の上で比較する2つのサンプルにとって共通の母集団の標準偏差と尖度になります。

尖度 \kappaの推定が問題になりますが、比較する2つの標本の平均がそれぞれゼロになるように、それぞれの標本の観測値をシフトした上で、2つの標本をプールした上で尖度の不偏推定量を計算します。こうすることで、平均がシフトしているだけで分散が同じ分布の比較で、帰無仮説を棄却しないようになります。

なお、参考にした元ネタではO'Neill (2014)のResult 15を参照しろと言う話になっていたのですが、Result 15はサンプルとサブサンプルの分散を比較する話で、サンプルとサンプルの分散を比較する話ではないので、そのまま利用しませんでした。

コード

kurtosis <- function(...){
    # 尖度の不偏推定量を計算（注意：正規分布の尖度が0ではなく3になる方）
    args <- list(...)
    kappa <- numeric(length(args))
    for(i in 1:length(args)){
        a <- args[[i]]
        n <- length(a)
        m <- mean(a)
        kappa[i] <- n * (n + 1) / (n - 1) * sum((a - m)^4) / sum((a - m)^2)^2
    }
    kappa
}

on2014ss2ss <- function(a, b, alternative = NULL){
    # サブサンプルとサブサンプルの分散比の統計量を、自由度補正してF検定
    n <- c(length(a), length(b))
    s2 <- c(var(a), var(b))
    kappa <- kurtosis(c(a - mean(a), b - mean(b)))
    DF_n <- function(n) 2*n / (kappa - (n - 3)/(n - 1))
    cpf <- function(i, j) pf(s2[i]/s2[j], DF_n(n[i]), DF_n(n[j]))

    if(is.null(alternative)) alternative <- "two.sided"
    (switch(alternative, "two.sided" = function(){
        2*min(cpf(1, 2), cpf(2, 1))
    }, "less" = function(){
        cpf(1, 2)
    }, "greater" = function(){
        cpf(2, 1)
    }))()
}

使い方は簡単で、以下のように二つのサンプルのベクトルを引数に与えるだけです。

n <- 30
a <- rgamma(n, shape = 2, scale = 3)
b <- rgamma(n, shape = 2, scale = 3)
on2014ss2ss(a, b)

シミュレーション

この方法（O’Neill）がどの程度有用かシミュレーションして確かめてみました。教科書的なF検定と、非正規分布でも信頼性が高いといわれるLevene検定と比較します。

まずは偽陽性が過剰にでないか、有意水準5%で同一の正規分布、一様分布、ラプラス分布、ガンマ分布、平均値だけシフトしたラプラス分布、平均値だけシフトした一様分布から生成した2つの乱数列で検定を行なってみました。

norm laplace laplace/shift unif unif/shift gamma
F-test 0.041 0.132 0.137 0.006 0.006 0.114
O’Neill 0.041 0.058 0.053 0.043 0.052 0.046
Levene 0.043 0.046 0.043 0.043 0.044 0.044

0.05が望ましい値ですが、教科書的なF-testは裾野が厚い一様分布では検出力が過小で、裾野が薄いラプラス分布では件出力が過剰になっています。今回紹介した方法とLeveneはこのような問題はありません。

2つの乱数列を生成した分布のパラメーター（と分散）が異なる場合も試してみましょう。

norm laplace unif gamma
F-test 0.833 0.681 0.553 0.966
O’Neill 0.824 0.536 0.797 0.914
Levene 0.661 0.440 0.449 0.895

こちらは高い方が望ましいと言うことになりますが、今回紹介した方法はLevene法よりも高くなっています。なお、教科書的なF-testは高いのですが、偽陽性も高いので頼りにはなりません。

なお、シミュレーションに使ったRのコードは以下になります。

library(VGAM) # ラプラス分布用
library(car) # Levene検定用
set.seed(547)
nr <- 1000 # 母集団分布×手法ごとの試行回数
n <- c(70, 30) # 試行ごとのサンプルサイズ
names_of_methods <- c("F-test", "O’Neill", "levene")

calcPs <- function(a, b){
    pv <- numeric(length(names_of_methods))
    n <- c(length(a), length(b))

    df01 <- data.frame(y = c(a, b), g = as.factor(c(rep("a", n[1]), rep("b", n[2]))))
    r <- leveneTest(y ~ g, data=df01)
    pv[3] <- r[[3]][1]

    # Old fashion F-test
    s2 <- c(var(a), var(b))
    df <- n - 1
    qf <- s2[1]/s2[2]
    pv[1] <- pf(qf, df[1], df[2])

    # Newbie based on O'Neill (2014)
    pv[2] <- on2014ss2ss(a, b)
    pv
}

simulate <- function(dgp1, dgp2, n, nr = 1000, alpha = 0.05){
    pv <- matrix(NA, nr, length(names_of_methods), dimnames = list(NULL, names_of_methods))
    for(i in 1:nr){
        pv[i, ] <- calcPs(dgp1(n[1]), dgp2(n[2]))
    }
    apply(pv, 2, function(x) sum(x < alpha))
}

r_fp <- sapply(c(
    function(n) rnorm(n, mean = 0, sd = 1),
    function(n) runif(n, min = 0, max = 1),
    function(n) rlaplace(n, location = 0, scale = 1),
    function(n) rgamma(n, shape = 2, scale = 3)
), function(dgp, n, nr){
    simulate(dgp, dgp, n, nr)
}, n = n, nr = nr)

r_fp <- cbind(r_fp, simulate(
    function(n) rlaplace(n, location = 0, scale = 1),
    function(n) rlaplace(n, location = 10, scale = 1),
    n = n, nr = nr
), simulate(
    function(n) runif(n, min = 0, max = 1),
    function(n) runif(n, min = 20, max = 21),
    n = n, nr = nr
))

r_fp <-  r_fp / nr
colnames(r_fp) <- c("norm", "unif", "laplace", "gamma", "laplace/shift", "unif/shift")
r_fp <- r_fp[, c(1, 3, 5, 2, 6, 4)]

r_pp <- cbind(simulate(
    function(n) rnorm(n, mean = 0, sd = 1),
    function(n) rnorm(n, mean = 0, sd = 1.5),
    n = n, nr = 1000
),simulate(
    function(n) rlaplace(n, location = 0, scale = 1),
    function(n) rlaplace(n, location = 0, scale = 1.5),
    n = n, nr = 1000
), simulate(
    function(n) rgamma(n, shape = 2, scale = 3),
    function(n) rgamma(n, shape = 2, scale = 6),
    n = n, nr = 1000
), simulate(
    function(n) runif(n, min = 0, max = 1),
    function(n) runif(n, min = -0.15, max = 1.15),
    n = n, nr = 1000
))

r_pp <- r_pp / nr
colnames(r_pp) <- c("norm", "laplace", "gamma", "unif")
rownames(r_pp) <- names_of_methods
r_pp <- r_pp[, c(1, 2, 4, 3)]

print("偽陽性率")
print(r_fp)

print("真陽性率（違う分布の値は比較不可）")
print(r_pp)