R言語で分散が既知のときの母分散の比の信頼区間を算出する関数やその実行例について解説します。
実行例では、データセットから母分散の比の信頼区間を計算するだけでなく、母分散の比の信頼区間の意味を図で説明します。
この記事で扱うプログラミングコードは以下からダウンロードできます。
母分散の比の信頼区間差については、次の記事を参照してください。
母分散の比の信頼区間
分散の信頼区間とR言語の関数についてまとめました。
母分散の比の信頼区間の概要について列挙しています。
母分散の比の信頼区間
正規分布\(N(\mu_1, \sigma_1^2)\)から大きさ\(n_1\)の無作為標本\(x_{11}, x_{12}, \ldots, x_{1n_1}\)、\(N(\mu_2, \sigma_2^2)\)から大きさ\(n_2\)の無作為標本\(x_{21}, x_{22}, \ldots, x_{2n_2}\)のが得られたとする。母分散の比\(\sigma_1^2 / \sigma_2^2\)の\(100(1- \alpha)\)%信頼区間を紹介する。以降、標本平均を\(\bar{x}_i = n_i^{-1} \sum_{j = 1}^{n_i} x_{ij}\)、不偏標本分散を\(s_i^2 = (n_i -1)^{-1} \sum_{j=1}^{n_i}(x_{ij} - \bar{x}_i)^2\)、第一自由度\(n_1\)、第二自由度\(n_2\)のF分布の上側\(\alpha\)点を\(F_{n_1, n_2, \alpha}\)とする。
母分散の比の信頼区間
母分散の比\(\sigma_1^2 / \sigma_2^2\)の\(100(1 - \alpha)\)%信頼区間は次で与えられる。
関数var.test
次の関数var.testを用いることで母分散の比の検定を行うことができます。 以下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 | データセットに欠損値がある場合に行うこと。 |
実行例
母分散の比の信頼区間の実行例について紹介します。
グラフの作成の際にパッケージEnvStats、ggplot2、gridExtraをを用います。あらかじめインストールしておいてください。
1 2 3 4 5 6 | install.packages("ggplot2") install.packages("gridExtra") install.packages("EnvStats") library(ggplot2) library(gridExtra) library(EnvStats) |
次の木に関するデータセットirisを例に母分散の比の信頼区間の計算を行っていきます。
1 2 3 4 5 6 7 8 9 | > dataset <- iris > head(dataset) 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 |
setosaとversicolor種のSepal.Lengthのヒストグラムは以下の通りです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | dataset_setosa_versic <- dataset[dataset$Species == "setosa" | dataset$Species == "versicolor", ] dataset_setosa_versic$Species <- factor(dataset_setosa_versic$Species) boxPlot <- ggplot(dataset_setosa_versic, aes(x = 1, y = Sepal.Length, color = Species, fill = Species)) + geom_boxplot() + scale_y_continuous() + coord_flip() + geom_point(position = position_dodge(width = 0.9), alpha = 0.3, size = 1) histPlot <- ggplot(dataset_setosa_versic, aes(x = Sepal.Length, color = Species, fill = Species)) + geom_histogram(bins = 20) + scale_y_continuous() grid.arrange(histPlot, boxPlot) |
setosaとversicolorのSepal.Lengthのヒストグラムと箱ひげ図
SpeciesがsetosaであるSepal.Lengthの母分散の比の信頼区間の計算を行いたいと思います。
setosaに限定する理由は、すべてのSpeciesを含めてしまうと正規性が成り立たなくなるためです。
ヒストグラムを描けば分かりますが、全てを含めた場合Sepal.Lengthのヒストグラムは多峰性を持ちます。
SpeciesがsetosaであるSepal.Lengthの95%信頼区間を算出するには、次のように関数varTestの引数conf.levelを0.95を与えて後ろに$conf.intを付けます。
3行目以降のように2群のデータをベクトルとしてvar.testの第一、第二引数に渡す方法でも計算することができます。
1 2 3 4 5 6 | conf_level <- 0.95 ci <- var.test(Sepal.Length ~ Species, data = dataset_setosa_versic , conf.level = conf_level)$conf.int sepalLen_setosa <- dataset$Sepal.Length[dataset$Species == "setosa"] sepalLen_versicolor <- dataset$Sepal.Length[dataset$Species == "versicolor"] var.test(sepalLen_setosa, sepalLen_versicolor, conf.level = conf_level)$conf.int |
ciを参照すると信頼区間の下限と上限、また信頼水準が格納されているのが分かります。
1 2 3 4 | > ci [1] 0.2646385 0.8217841 attr(,"conf.level") [1] 0.95 |
信頼区間の上限と下限を取得するには、次のようにciの要素をインデックスで取得すれば大丈夫です。
1 2 | ci_lower <- ci[1] #信頼区間の下限 ci_upper <- ci[2] #信頼区間の上限 |
計算した信頼区間を他のファイルに保存したいときは、次のように関数data.frameでデータフレームに計算した値を格納しwrite.csvなどで他のファイルへ保存します。
1 2 3 4 | resultTable <- data.frame(sample.var.setosa = var(dataset$Sepal.Length[dataset$Species == "setosa"]), sample.var.versicolor = var(dataset$Sepal.Length[dataset$Species == "versicolor"]), c.i.lower = ci[1], c.i.upper = ci[2], row.names = NULL) write.csv(resultTable, "母分散の比の信頼区間.csv", row.names = FALSE) |
母分散の比の信頼区間
上で算出した信頼区間をグラフにすると下の画像のようになります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | ci_setosa <- varTest(dataset$Sepal.Length[dataset$Species == "setosa"], conf.level = conf_level)$conf.int ci_versicolor <- varTest(dataset$Sepal.Length[dataset$Species == "versicolor"], conf.level = conf_level)$conf.int dataset_ci <- data.frame(x = c(rep(tapply(dataset_setosa_versic$Sepal.Length, dataset_setosa_versic$Species, mean), 2) + c(tapply(dataset_setosa_versic$Sepal.Length, dataset_setosa_versic$Species,sd), -tapply(dataset_setosa_versic$Sepal.Length, dataset_setosa_versic$Species,sd)), rep(tapply(dataset_setosa_versic$Sepal.Length, dataset_setosa_versic$Species, mean), each = 2) + c(sqrt(c(ci_setosa, ci_versicolor)), -sqrt(c(ci_setosa, ci_versicolor)))), name = factor(c(rep("SD", 4), paste0(100 * conf_level, paste0("% c.i.", rep(c("lower", "upper"), 4)))), levels = c("SD", paste0(100 * conf_level, "% c.i.", c("lower", "upper")))), group = factor(c(rep(levels(dataset_setosa_versic$Species), 2), rep(levels(dataset_setosa_versic$Species), each = 2, times = 2)))) ggplot(dataset_setosa_versic, aes(x = Species, y = Sepal.Length, color = Species)) + geom_dotplot(binaxis = "y", stackdir = "center") + stat_summary(aes(x = Species, y = Sepal.Length), fun.min = function(x) mean(x) - sd(x), fun.max = function(x) mean(x) + sd(x), geom = 'errorbar', width = 0.15, color = "#00BA38", size = 1) + stat_summary(aes(x = as.numeric(Species) - 0.05, y = Sepal.Length), fun.min = function(x) mean(x) - sqrt(varTest(x)$conf.int[1]), fun.max = function(x) mean(x) + sqrt(varTest(x)$conf.int[1]), geom = 'errorbar', width = 0.1, color = "#619CFF", size = 1) + stat_summary(aes(x = as.numeric(Species) - 0.1, y = Sepal.Length), fun.min = function(x) mean(x) - sqrt(varTest(x)$conf.int[2]), fun.max = function(x) mean(x) + sqrt(varTest(x)$conf.int[2]), geom = 'errorbar', width = 0.1, color = "#F8766D", size = 1) + geom_text(data = dataset_ci, size = 3.5, aes(x = as.numeric(group) - 0.35, y = x, color = name, label = paste0(c(rep("SD", 4), rep(c("c.i.lower", "c.i.upper"), 4)), ": ", round(x, 2)))) + labs(color = "") + scale_colour_manual(breaks = c("SD", "95% c.i.lower", "95% c.i.upper"), values = c("#00BA38", "#619CFF", "#F8766D", "#000000", "#000000")) |
2群の母分散の信頼区間
- 緑線はsetosa及びversicolor種のSepal.Lengthの標本平均±標準偏差
- 赤線はsetosa及びversicolor種のSepal.Lengthの標本平均±標準偏差の信頼区間の上限
- 青線はsetosa及びversicolor種のSepal.Lengthの標本平均±標準偏差の信頼区間の下限
を表しています。
setosa種のSepal.Lengthの母分散の信頼区間の下限をversicolor種の信頼区間の上限で割ったものがだいたい0.20、
setosa種のSepal.Lengthの母分散の信頼区間の上限をversicolor種の信頼区間の下限で割ったものがだいたい1.04
であるように、上で求めた母分散の比の信頼区間は、グラフのsetosaの赤線と青線とversicolorの青線と赤線の比率に近い値になることが確認できます。（上のグラフの信頼区間は母分散の比ではなく母分散に基づくので、計算した値とは異なります）
まとめ
R言語で母分散の比の信頼区間を計算する関数やその実行例をみていきました。
関数var.testを用いることで母分散の比の信頼区間を計算することができます。
