【R言語】アンサリ・ブラッドリー検定　ansari.testの使い方

R言語でアンサリ・ブラッドリー検定を実行する方法を解説していきます。

この記事では、非正規母集団からの標本に対する分散の比の検定であるアンサリ・ブラッドリー検定のための関数とその実行例を紹介します。

F検定と同様に、アンサリ・ブラッドリー検定の関数が標準で実装されており、簡単に実行することができます。

また、F検定については次の記事で解説しています。

: 【R言語】分散の比の検定・F検定　関数var.testの使い方
R言語でF検定を実行する方法を解説していきます。 2群間の等分散性の検定であるF検定を実行する関数とその実行例について見ていきます。 F検定の詳細については、下の記事で解説しています。また、この記事 ...
続きを見る

また、この記事で紹介するプログラミングコードは以下のzipファイル中のスクリプトに保存しています。

R言語　アンサリ・ブラッドリー検定

1 ファイル 1.09 KB

ダウンロード

改訂2版 RユーザのためのRStudio[実践]入門〜tidyverseによるモダンな分析フローの世界

RStudioで解析したい方へ

¥3,278 （2022/06/02 22:52時点 | Amazon調べ）

Amazon

楽天市場

Yahooショッピング

ポチップ

アンサリ・ブラッドリー検定

分散の検定に用いる様々な関数を紹介します。ここでは関数の用法とその引数の説明をします。関数の使用例については、この後の実行例で行っています。

概要

$x_1, x_2, \ldots, x_m$を第一母集団$\Pi_1$、$y_1, y_2, \ldots, y_n$を第二母集団$\Pi_2$からの無作為標本とする。また、標本に対応する確率変数をそれぞれ$X_1, X_2, \ldots, X_m$、$Y_1, Y_2, \ldots, Y_n$とし、2つの母集団の中央値は等しいとする。それぞれの母集団の尺度パラメータを$\eta_1$、$\eta_2$とする。このとき、次の「2つの母集団の分散が等しいか」という仮説検定を考える。

\begin{align}&H_0:\ \eta_1 / \eta_2 = 1 \\ &H_1:\ \eta_1 / \eta_2 \neq 1\end{align}

$R_i,\ i = 1, 2, \ldots, m$を、次のように$X_1, X_2, \ldots, X_m$、$Y_1, Y_2, \ldots, Y_n$を順番に並べ中央の順位に向かって増加するような$X_1, X_2, \ldots, X_m$の順位とする。

\begin{align}R_i \in \left\{ \begin{array}{cc}\{1, 2, \ldots,(m+n)/2, (m+n)/2, \ldots, 2, 1\}, & \mathrm{if}\ m + n\ \mathrm{is}\ \mathrm{even},\\ \{1, 2, \ldots,(m+n - 1)/2, (m+n + 1)/2, (m+n - 1)/2, \ldots, 2, 1\}, & \mathrm{if}\ m + n\ \mathrm{is}\ \mathrm{odd}.\end{array} \right. \end{align}

$m, n \to \infty$の下で検定統計量は次で与えられる。

\begin{align}T = \cfrac{A - \mathrm{E}[A]}{\sqrt{\mathrm{Var}[A]}} \sim N(0, 1), \end{align}

ここに

\begin{align}A &= \sum_{i=1}^m R_i ,\\ \mathrm{E}[A] &=\left\{\begin{array}{cc} \cfrac{m(m +n + 2) }{4}, & \mathrm{if}\ m + n \ \mathrm{is}\ \mathrm{even}, \\ \cfrac{m(m +n + 1)^2 }{4(m + n)}, & \mathrm{if}\ m + n \ \mathrm{is}\ \mathrm{odd} , \end{array} \right. \\ \mathrm{Var}[A] &= \left\{\begin{array}{cc} \cfrac{mn(m + n + 2)(m + n -2)}{48(m + n - 1)}, & \mathrm{if}\ m + n \ \mathrm{is}\ \mathrm{even},\\\cfrac{mn(m +n +1)\bigl\{(m + n )^2 + 3\bigr\}}{48(m + n )^2} , & \mathrm{if}\ m + n \ \mathrm{is}\ \mathrm{odd}.\end{array} \right.\end{align}

また、有意水準$\alpha$の棄却域は次で与えられる。

\begin{align}(-\infty , -z_{\alpha/ 2}) \cup (z_{\alpha / 2}, \infty),\end{align}

ここに$z_{\alpha}$は標準正規分布の上側$\alpha$点である。

関数ansari.test

次の関数ansari.testを使うことで、アンサリ・ブラッドリー検定を行うことができます。

ansari.test(x, y, alternative = c("two.sided", "less", "greater"), exact = NULL, conf.int = FALSE, conf.level = 0.95, …)

ansari.test(formula, data, subset, na.action, …)

関数ansari.testの引数
x	片方の群のnumric型のベクトルのデータ。
y	もう一方の群のnumric型のベクトルのデータ。
alternative	帰無仮説で与える分散の比の値。デフォルトで1。alternative対立仮説の設定。"two-sided"で両側検定、"greater"で右片側検定、"less"で左片側検定を指定できる。
exact	正確なp値を計算するかどうかをlogical型で指定できる。
conf.int	信頼区間を計算するかどうかをlogical型で指定できる。
conf.level	信頼区間の信頼度。
formula	式。次の引数dataで与えたデータセットに列名aとbがあるとき、bの水準に関するaの分散の比の検定がしたい場合a~bとする。
data	データセット。
subset	検定に用いる標本の部分集合をベクトルで指定できる。
na.action	データセットに欠損値がある場合に行うこと。

実行例

関数ansari.testを使ってアンサリ・ブラッドリー検定を行っていきます。

データセットをプロットする際に以下のパッケージが必要なため、事前にインストールしたライブラリに追加しときます。

library(ggplot2)
library(gridExtra)

1 2	library(ggplot2) library(gridExtra)

まず、データセットとして次を用います。

ramsay <- c(111, 107, 100, 99, 102, 106, 109, 108, 104, 99,
            101, 96, 97, 102, 107, 113, 116, 113, 110, 98)
jung.parekh <- c(107, 108, 106, 98, 105, 103, 110, 105, 104,
                 100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99)

ramsay <- c(111, 107, 100, 99, 102, 106, 109, 108, 104, 99,

101, 96, 97, 102, 107, 113, 116, 113, 110, 98)

jung.parekh <- c(107, 108, 106, 98, 105, 103, 110, 105, 104,

100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99)

ansari.testのRDocumentationにもあるように　赤血球に関するデータとなります。

解析を行う準備としてこれらのデータを次のようにデータセットに変換しておきます。

> dataset <- data.frame(serum.iron = c(ramsay, jung.parekh),
+                       group = rep(c("ramsay", "jung.parekh"), each = length(ramsay)))
> head(dataset)
  serum.iron  group
1        111 ramsay
2        107 ramsay
3        100 ramsay
4         99 ramsay
5        102 ramsay
6        106 ramsay

> dataset <- data.frame(serum.iron = c(ramsay, jung.parekh),

+ group = rep(c("ramsay", "jung.parekh"), each = length(ramsay)))

> head(dataset)

serum.iron group

1 111 ramsay

2 107 ramsay

3 100 ramsay

4 99 ramsay

5 102 ramsay

6 106 ramsay

次のヒストグラムから分かるように、正規性を持たないことが予想され、F検定などのパラメトリック手法が適用できないことが分かります。

densityPlot <- ggplot(dataset, aes(x = serum.iron , y = ..density.., colour = group, fill = group)) +
  geom_histogram(position = "identity", bins = 25, alpha = 0.75) +
  geom_density(stat = "density", position = "identity", alpha = 0.75) +
  xlim(min(dataset$serum.iron ), max(dataset$serum.iron ))

boxPlot <- ggplot(dataset, aes(x = serum.iron , y = group, fill = group)) +
             geom_boxplot(alpha = 0.75, colour = "gray")

grid.arrange(densityPlot, boxPlot) #経験分布と箱ひげ図を同時にプロット

densityPlot <- ggplot(dataset, aes(x = serum.iron , y = ..density.., colour = group, fill = group)) +

geom_histogram(position = "identity", bins = 25, alpha = 0.75) +

geom_density(stat = "density", position = "identity", alpha = 0.75) +

xlim(min(dataset$serum.iron ), max(dataset$serum.iron ))

boxPlot <- ggplot(dataset, aes(x = serum.iron , y = group, fill = group)) +

geom_boxplot(alpha = 0.75, colour = "gray")

grid.arrange(densityPlot, boxPlot) #経験分布と箱ひげ図を同時にプロット

ヒストグラムと箱ひげ図

このように正規性を持たないまたは、2群のデータの中央値がだいたい同じであるデータ対してはアンサリ・ブラッドリー検定のようなノンパラメトリック検定を用いるのが適しています。

次のように関数ansari.testの引数dataにデータセット、formulaに式（データの変数 ~ 群の変数）を記述することで、アンサリ・ブラッドリー検定を実行することができます。

3行目以降のように、2つのグループのデータをあらかじめ与える方法でも実行することができます。

testResult <- ansari.test(serum.iron  ~ group, data = dataset)

serum.iron_ramsay <- dataset$serum.iron [dataset$group == "ramsay"]
serum.iron_jung.parekh <- dataset$serum.iron [dataset$group == "jung.parekh"]
testResult <- ansari.test(serum.iron_ramsay, serum.iron_jung.parekh)

testResult <- ansari.test(serum.iron ~ group, data = dataset)

serum.iron_ramsay <- dataset$serum.iron [dataset$group == "ramsay"]

serum.iron_jung.parekh <- dataset$serum.iron [dataset$group == "jung.parekh"]

testResult <- ansari.test(serum.iron_ramsay, serum.iron_jung.parekh)

testResultを参照すると次のように検定の結果がコンソールに表示されます。

> testResult

	Ansari-Bradley test

data:  serum.iron by group
AB = 234.5, p-value = 0.1815
alternative hypothesis: true ratio of scales is not equal to 1

> testResult

Ansari-Bradley test

data: serum.iron by group

AB = 234.5, p-value = 0.1815

alternative hypothesis: true ratio of scales is not equal to 1

ポイント

p値=0.1815より「ramsayとjung.parekhの分散のは異なるとは言えない」ことが言えます。すなわち、2群間の分散が等しいことが分かりました。

検定結果を次のようにして取得することが可能です。

stat <- unname(testResult$statistic) #統計量
pValue <- testResult$p.value         #p値

1 2	stat <- unname(testResult$statistic) #統計量 pValue <- testResult$p.value #p値

最後に、これらの検定結果をcsvファイルへ出力する例を以下に紹介します。

次のようにしてresultTableに格納された検定結果をcsvファイルに出力することができます。

resultTable <- data.frame(stat = stat, p.value = pValue, row.names = NULL)
write.csv(resultTable, "アンサリ・ブラッドリー検定.csv", row.names = FALSE)

1 2	resultTable <- data.frame(stat = stat, p.value = pValue, row.names = NULL) write.csv(resultTable, "アンサリ・ブラッドリー検定.csv", row.names = FALSE)

アンサリ・ブラッドリー検定結果

また、片側検定についても解説します。

次のt分布の乱数を使って片側検定を実行する例をみていきます。

x <- rt(40, 20)
y <- rt(40, 4)

dataset <- data.frame(sample = c(x, y), group = rep(c("1st", "2nd"), length(x)))

x <- rt(40, 20)

y <- rt(40, 4)

dataset <- data.frame(sample = c(x, y), group = rep(c("1st", "2nd"), length(x)))

乱数によって変わりますが次のような自由度の異なるt分布からのデータから成ります。

densityPlot <- ggplot(dataset, aes(x = sample , y = ..density.., colour = group, fill = group)) +
  geom_histogram(position = "identity", bins = 25, alpha = 0.75) +
  geom_density(stat = "density", position = "identity", alpha = 0.75) +
  xlim(min(dataset$sample), max(dataset$sample))

boxPlot <- ggplot(dataset, aes(x = sample , y = group, fill = group)) +
  geom_boxplot(alpha = 0.75, colour = "gray")

grid.arrange(densityPlot, boxPlot)

densityPlot <- ggplot(dataset, aes(x = sample , y = ..density.., colour = group, fill = group)) +

geom_histogram(position = "identity", bins = 25, alpha = 0.75) +

geom_density(stat = "density", position = "identity", alpha = 0.75) +

xlim(min(dataset$sample), max(dataset$sample))

boxPlot <- ggplot(dataset, aes(x = sample , y = group, fill = group)) +

geom_boxplot(alpha = 0.75, colour = "gray")

grid.arrange(densityPlot, boxPlot)

ヒストグラムと箱ひげ図

アンサリ・ブラッドリー検定により、第一母集団（1st）の群の分散がと第二母集団（2nd）の群の分散よりも大きいかどうか検定するには、次のように関数ansari.testの引数alternativeを"greater"にします。

> testResult <- ansari.test(x, y, alternative = "greater")
> testResult

	Ansari-Bradley test

data:  x and y
AB = 914, p-value = 0.9656
alternative hypothesis: true ratio of scales is greater than 1

> testResult <- ansari.test(x, y, alternative = "greater")

> testResult

Ansari-Bradley test

data: x and y

AB = 914, p-value = 0.9656

alternative hypothesis: true ratio of scales is greater than 1

ポイント

p-value = 0.9656より「第一母集団の群の分散は第二母集団よりも大きいこと」がいえました。（第一母集団の群の分散は第二母集団よりも大きくないとは言えない）

両側検定の時と同様に、検定の結果を抽出または保存したいときは次のようにします。

stat <- unname(testResult$statistic) #統計量
pValue <- testResult$p.value         #p値

resultTable <- data.frame(stat = stat, p.value = pValue, row.names = NULL)