【R言語】F検定/分散の比の検定　関数var.testの使い方

R言語でF検定を実行する方法を解説していきます。

2群間の等分散性の検定であるF検定を実行する関数とその実行例について見ていきます。

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

: 【統計学】等分散性のF検定・分散の比の検定
2つの母集団の分散の比ついての仮説を検定する際に用いられるF検定を解説する。 F検定の検定統計量の導出や検定統計量がF分布に従うことの証明をしていく。また、母集団の分散の比を検定しt検定を行う際の検 ...
続きを見る

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

R言語　F検定

1 ファイル 1.36 KB

ダウンロード

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

RStudioで解析したい方へ

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

Amazon

＼楽天ポイント5倍セール！／

楽天市場

Yahooショッピング

ポチップ

F検定

F検定の概要とR言語の関数を以下にまとめました。

概要

F検定

$x_{11}, \ldots, x_{1n_1}$は$N(\mu, \sigma_1^2)$からの独立同一な標本であるとし、$x_{21}, \ldots, x_{2n_2}$は$N(\mu, \sigma_2^2)$からの独立同一な標本であるとする。このとき、次の「2つの母集団の分散$\sigma_1^2$と$\sigma_2^2$は等しいか」の仮説を検定する。

\begin{align}&H_0: \sigma_1^2 = \sigma_2^2\\&H_1: \sigma_1^2 \neq \sigma_2^2\end{align}

検定統計量として次を用いる。

\begin{align} F = \cfrac{s_1^2}{s_2^2}\sim F_{n_2-1}^{n_1-1},\end{align}

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

\begin{align}[0, F_{n_2-1, 1-\alpha/2}^{n_1-1}) \cup (F_{n_2-1, \alpha/2}^{n_1-1}, \infty),\end{align}

ここに、$s_1^2$と$s_2^2$は次で定義される不偏標本分散の確率変数である。

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

関数var.test

次の関数var.testを用いることで、F検定を実行することができます。

var.test(x, y, ratio = 1, alternative = c("two.sided", "less", "greater"), conf.level = 0.95, …)

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

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

実行例

Ｆ検定の実行例を紹介します。

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

install.packages("ggplot2")
install.packages("gridExtra")
library(ggplot2)
library(gridExtra)

install.packages("ggplot2")

install.packages("gridExtra")

library(ggplot2)

library(gridExtra)

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

> 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

データirisには3つの水準をもつSpeciesという列があります。

Speciesの水準に関するSepal.Lengthの経験分布と箱ひげ図は次のようになります。

dataset_setosa_versic <- dataset[dataset$Species == "setosa" | dataset$Species == "versicolor", ]
  
densityPlot <- ggplot(dataset_setosa_versic, 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_setosa_versic$Sepal.Length), max(dataset_setosa_versic$Sepal.Length))

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

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

dataset_setosa_versic <- dataset[dataset$Species == "setosa" | dataset$Species == "versicolor", ]

densityPlot <- ggplot(dataset_setosa_versic, 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_setosa_versic$Sepal.Length), max(dataset_setosa_versic$Sepal.Length))

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

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

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

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

SetosaとVersicolorの箱ひげ図のひげの長さが異なるため、例として「SetosaとVersicolorか」という仮説検定をしたいと思います。

等分散性のF検定は次のように関数var.testの引数dataにデータセット、formulaに変数名~グループ名を渡すことで実行することができます。

また、3行目以降のように第1引数と第2引数に2つの群のデータを渡す方法でも実行できます。

testResult <- var.test(Sepal.Length ~ Species, data = dataset_setosa_versic)

sepalLength_setosa <- data$Sepal.Length[data$Species == "setosa"]
sepalLength_versicolor <- data$Sepal.Length[data$Species == "versicolor"] 
var.test(sepalLength_setosa, sepalLength_versicolor)

testResult <- var.test(Sepal.Length ~ Species, data = dataset_setosa_versic)

sepalLength_setosa <- data$Sepal.Length[data$Species == "setosa"]

sepalLength_versicolor <- data$Sepal.Length[data$Species == "versicolor"]

var.test(sepalLength_setosa, sepalLength_versicolor)

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

> testResult

	F test to compare two variances

data:  Sepal.Length by Species
F = 0.46634, num df = 49, denom df = 49, p-value = 0.008657
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.2646385 0.8217841
sample estimates:
ratio of variances 
         0.4663429

> testResult

F test to compare two variances

data: Sepal.Length by Species

F = 0.46634, num df = 49, denom df = 49, p-value = 0.008657

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

95 percent confidence interval:

0.2646385 0.8217841

sample estimates:

ratio of variances

0.4663429

ポイント

p値=0.008657<0.05より、「SetosaとVersicolorの群のSepal.Lengthの分散は等しいとは言えない」という結果が得られ、SetosaとVersicolor間で分散が異なることが分かりました。

上の実行結果に書いてある検定結果は次のようにtestResultの後ろに$を付けることで取得することが可能です。

stat <- unname(testResult$statistic) #統計量
df <- unname(testResult$parameter)   #自由度
pValue <- testResult$p.value         #p値
ci <- testResult$conf.int            #母分散の比の信頼区間
estimate <- testResult$estimate      #母分散の比の推定値
nullValue <- testResult$null.value   #帰無仮説の下での母分散の比の値

stat <- unname(testResult$statistic) #統計量

df <- unname(testResult$parameter) #自由度

pValue <- testResult$p.value #p値

ci <- testResult$conf.int #母分散の比の信頼区間

estimate <- testResult$estimate #母分散の比の推定値

nullValue <- testResult$null.value #帰無仮説の下での母分散の比の値

csvファイル等に保存したいときは、以下のように関数data.frameを用いこれらの検定結果の値をデータフレームに格納し、write.tableやwrite.csvなどでcsvファイルへ出力することができます。

resultTable <- data.frame(stat = stat, num.d.f. = df[1], denom.d.f. = df[2], p.value = pValue,
                          estimate = estimate, c.i.lower = ci[1], c.i.upper = ci[2],
                          null.value = nullValue, row.names = NULL)

write.csv(resultTable, "F検定.csv", row.names = FALSE)

resultTable <- data.frame(stat = stat, num.d.f. = df[1], denom.d.f. = df[2], p.value = pValue,

estimate = estimate, c.i.lower = ci[1], c.i.upper = ci[2],

null.value = nullValue, row.names = NULL)

write.csv(resultTable, "F検定.csv", row.names = FALSE)

作業ディレクトリに次の"F検定.csv"が保存されます。

F検定結果

次に片側検定の例について見ていきます。

データセットに次の正規乱数のデータxとyを用います。

#片側検定
x <- rnorm(40, sd = 2)
y <- rnorm(60)

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

#片側検定

x <- rnorm(40, sd = 2)

y <- rnorm(60)

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

以下、xとyのヒストグラムと箱ひげ図です。xの方が標準偏差が大きいのが分かります。

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

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

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

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

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

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

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

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

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

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

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

xとyのヒストグラムと箱ひげ図

「xの母分散はyの母分散よりも大きい」という仮説検定を実行するには、次のように関数var.testの引数alternativeを"greater"にします。

testResult <- var.test(data ~ group, data = dataset, alternative = "greater")

1	testResult <- var.test(data ~ group, data = dataset, alternative = "greater")

testResultを参照すると片側検定の結果が出力されます。

> testResult

	F test to compare two variances

data:  data by group
F = 4.2596, num df = 39, denom df = 59, p-value = 2.967e-07
alternative hypothesis: true ratio of variances is greater than 1
95 percent confidence interval:
 2.658963      Inf
sample estimates:
ratio of variances 
          4.259621

> testResult

F test to compare two variances

data: data by group

F = 4.2596, num df = 39, denom df = 59, p-value = 2.967e-07

alternative hypothesis: true ratio of variances is greater than 1

95 percent confidence interval:

2.658963 Inf

sample estimates:

ratio of variances

4.259621

ポイント

p値=2.967e-07<0.05より、「xとyの母分散は等しいとは言えない」という結果が得られ、xの母分散はyの母分散よりも大きいことが分かりました。

両側検定と同様に、上で計算した検定結果を次のようにして取得することが可能です。

stat <- unname(testResult$statistic) #統計量
df <- unname(testResult$parameter)   #自由度
pValue <- testResult$p.value         #p値
ci <- testResult$conf.int            #母分散の比の信頼区間
estimate <- testResult$estimate      #母分散の比の推定値
nullValue <- testResult$null.value   #帰無仮説の下での母分散の比の値

resultTable <- data.frame(stat = stat, num.d.f. = df[1], denom.d.f. = df[2], p.value = pValue,
                          estimate = estimate, c.i.lower = ci[1], c.i.upper = ci[2],
                          null.value = nullValue, row.names = NULL)