【R言語】バートレット検定　bartlett.testの使い方

R言語でバートレット検定を実行する関数とその実行例について解説していきます。

この記事では、多群の等分散性の検定であるバートレット検定の概要やR言語での実行例を扱っています。

2群のF検定については、下の記事で解説しています。

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

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

R言語　バートレット検定

1 ファイル 0.77 KB

ダウンロード

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

RStudioで解析したい方へ

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

Amazon

楽天市場

Yahooショッピング

ポチップ

バートレット検定

バートレット検定の概要とR言語の関数をまとめました。

概要

標本$x_{i1}, x_{i2}, \ldots, x_{in_i},\ i= 1, 2, \ldots, k$を$N(\mu_i, \sigma_i^2)$からの無作為標本とする。また、これらの標本に対応する確率変数を$X_{i1}, X_{i2}, \ldots, X_{in_i},\ i= 1, 2, \ldots, k$とする。次の「すべての母集団の分散が等しい」という仮説検定を考える。

\begin{align}&H_0:\ \sigma_1^2 = \sigma_2^2 = \cdots = \sigma_k^2\\ &H_1:\ \mathrm{not}\ H_0\end{align}

この仮説の検定統計量は次で与えられる。

\begin{align}X^2 &= \cfrac{(N-k) \log S^2+ \sum_{i = 1}^k (n_i - 1)\log S_i^2}{\left. 1 + \Bigl[ \sum_{i = 1}^k \bigl\{1 / (n_i - 1)\bigr\} - 1 / (N- k ) \Bigr\} \right/ 3(k - 1)} , \end{align}

ここに$N = \sum_{i = 1}^k n_i$、

\begin{align}S_i &= \cfrac{1}{n_i-1} \sum_{j=1 }^{n_i}(X_{ij} - \bar{X}_i)^2,\ i = 1, 2, \ldots, k, \\ \bar{X}_i &= \cfrac{}{} \sum_{j = 1}^{n_i} X_{ij},\ i = 1,2, \ldots, k,\\ S^2 &= \cfrac{1}{N - k}\sum_{i=1}^k (n_i - 1)S_i^2 .\end{align}

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

\begin{align} (\chi_{k - 1, \alpha}^2, \infty),\end{align}

ここに、$\chi_{k, \alpha}^2$は自由度$k$のカイ二乗分布の上側$\alpha$点。

関数bartlett.test

次の関数bartlett.testを使うことで、バートレット検定を行うことができます。

bartlett.test(x, g, …) bartlett.test(formula, data, subset, na.action, …)

関数bartlett.testの引数
x	数値データのベクトル、またはそれぞれの水準を表すリスト型のデータ、もしくは線形モデルのオブジェクト。
g	群に対応するfactor型のベクトル。
formula	式。次の引数dataで与えたデータセットに列名aとbがあるとき、bの水準に関するaの分散の比の検定がしたい場合a~bとする。
data	データセット。
subset	検定に用いる標本の部分集合をベクトルで指定できる。
na.action	データセットに欠損値がある場合に行うこと。

実行例

バートレット検定の実行例について見ていきます。

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

library(ggplot2)
library(gridExtra)

1 2	library(ggplot2) library(gridExtra)

では、実際にバートレット検定を行っていきます。

データセットには次のirisを使います。SpeciesはSetosa、Versicolor、Virginicaという3つの水準があるため、これら3群間の分散の検定を行いたいと思います。

> data <- iris
> head(data)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

> data <- iris

> head(data)

Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1 5.1 3.5 1.4 0.2 setosa

2 4.9 3.0 1.4 0.2 setosa

3 4.7 3.2 1.3 0.2 setosa

4 4.6 3.1 1.5 0.2 setosa

5 5.0 3.6 1.4 0.2 setosa

6 5.4 3.9 1.7 0.4 setosa

各SpeciesごとのSepal.Lengthのヒストグラムと箱ひげ図を描くには次のように関数ggplotを実行します。

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

boxPlot <- ggplot(dataset, aes(x = Sepal.Length, y = Species, fill = Species)) + 
             geom_boxplot(alpha = 0.75, colour = "gray") +
             xlim(min(dataset$Sepal.Length), max(dataset$Sepal.Length))

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

densityPlot <- ggplot(dataset, aes(x = Sepal.Length, y = ..density.., colour = Species, fill = Species)) +

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

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

xlim(min(dataset$Sepal.Length), max(dataset$Sepal.Length))

boxPlot <- ggplot(dataset, aes(x = Sepal.Length, y = Species, fill = Species)) +

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

xlim(min(dataset$Sepal.Length), max(dataset$Sepal.Length))

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

ヒストグラムと箱ひげ図

バートレット検定によって「SpeciesによってSepal.Lengthの分散が等しい」という仮説検定を行いたいと思います。

次のように関数bartlett.testのの第1引数にデータ、第2引数にラベル（factor型のベクトル）を指定することで、バートレット検定をことができます。

また、2行目のように引数dataとformulaを使って実行することもできます。

> testResult <- bartlett.test(data$Sepal.Length, data$Species)
> bartlett.test(Sepal.Length ~ Species, data = dataset)
> testResult

	Bartlett test of homogeneity of variances

data:  data$Sepal.Length and data$Species
Bartlett's K-squared = 16.006, df = 2, p-value = 0.0003345

> testResult <- bartlett.test(data$Sepal.Length, data$Species)

> bartlett.test(Sepal.Length ~ Species, data = dataset)

> testResult

Bartlett test of homogeneity of variances

data: data$Sepal.Length and data$Species

Bartlett's K-squared = 16.006, df = 2, p-value = 0.0003345

ポイント

上の結果から、p値=0.0003345<0.05より「Speciesの3水準間のSepal.Lengthの分散は全て等しくい」という仮説は棄却され、3水準間の分散は全て等しくないことが分かりました。

bartlett.testによる検定結果も以下のようにして取得することができます。検定統計量やp値などの検定結果を抽出することができます。

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

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

最後に検定結果をcsvファイルへ出力する例を紹介します。

次のように、関数data.frameを用いて検定結果（検定統計量やp値）データフレームに格納し、write.csvを用いてデータフレームを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)